Solution Procedure for Transport Modeling in Effluent Recharge Based on Operator-Splitting Techniques

Abstract

The coupling of groundwater movement and reactive transport during groundwater recharge with wastewater leads to a complicated mathematical model, involving terms to describe convection-dispersion, adsorption/desorption and/or biodegradation, and so forth. It has been found very difficult to solve such a coupled model either analytically or numerically. The present study adopts operator-splitting techniques to decompose the coupled model into two submodels with different intrinsic characteristics. By applying an upwind finite difference scheme to the finite volume integral of the convection flux term, an implicit solution procedure is derived to solve the convection-dominant equation. The dispersion term is discretized in a standard central-difference scheme while the dispersion-dominant equation is solved using either the preconditioned Jacobi conjugate gradient (PJCG) method or Thomas method based on local-one-dimensional scheme. The solution method proposed in this study is applied to the demonstration project of groundwater recharge with secondary effluent at Gaobeidian sewage treatment plant (STP) successfully.

1. Introduction

Usually, the groundwater system can self-maintain a hydraulic steady-state and chemical equilibrium before some external stimulations, such as rains, irritation, or wastewater recharge, are introduced. However, the groundwater system will experience a long time unsteady process after such stimulations and reach new equilibrium state again. In the unsteady duration caused by external stimulations, especially by wastewater recharging, the subsurface system undergoes complicated physical/chemical/biological processes of solute contaminants: convection-dispersion, adsorption/desorption, dissolution/precipitation, geo-chemical reactions, and/or biodegradations. Each of theses processes can be described by suitable mathematical model. The coupling of the above processes is called the reactive transport of contaminants. The mathematical models describing the coupled phenomena of some or all processes mentioned above have been proved extremely complicated. Furthermore, to obtain solutions of such models has been proved a tough mission. Many investigations have been carried out in recent years to investigate the phenomena on hydrophysical processes with or without reactive transport in variety of porous media ([1–16]). Because several processes get involved in the groundwater recharge with wastewater, the mathematical models describing the coupled system show strongly nonlinear and can only be solved by numerical means.

Numerical solution methods for the suite of reactive transport equations may be divided into two categories: one-step approaches and two-step approaches. Meanwhile, the standard Galerkin finite element method (FEM) or its modified version is adopted for the discretization of the flow and the solute transport equations by most of the existing numerical models. In the one-step approaches, the whole equation system describing the coupled model of groundwater flow and reactive transport is solved simultaneously, while in the two-step approaches the flow/transport and reaction portions are treated separately. The one-step methods have been proved inefficient and seldom been used to solve practical problems. Most of investigations focus on the two-step or multistep methods. Few of them employed the two-step methods based on sequential iteration approach (SIA) to solve reactive transport model. D. Schäfer et al. [5] developed a 3-D model transport, biochemistry, and chemistry (TBC) to deal with reactive transport in a saturated groundwater flow. A sequential solution procedure is used where the advection-dispersion transport equations are solved independently for each species or aqueous component by adding source/sink term representing exchanges with other phases and chemical/biochemical reactions from preceding time-step. The chemical/biochemical reaction equations are solved with the concentration changes from the transport as explicit source/sink term. Tebes-Stevens et al. [8] developed a reactive transport model FEREACT to examine the coupled effects of 2-D steady-state groundwater flow, equilibrium aqueous speciation reactions, and kinetically controlled interphase reactions. An improved two-step solution algorithm of SIA is used to incorporate the effects of geochemical and microbial reaction processes in governing equations for solute transport in subsurface.

In many other existing investigations, two-step approaches based on operator-splitting techniques are adopted to solve the reactive transport model. Prior to these researches, the operator-splitting techniques were applied to solve the convection-dispersion equation and presented its advantages. Zaidel and Levi [17] proposed a second-order accurate explicit finite difference scheme based on this method to obtain the numerical solution of the convection equation. Russo et al. [18] used a standard central-difference scheme to solve the dispersion equation. In recent years, the research activities have been concentrated on the numerical solution approaches of reactive transport models by using the operator-splitting techniques. Barry et al. [19] developed an alternative operator-splitting approach for solving chemical reaction groundwater transport models. Cheng and Yeh [4] presented a 3-D numerical model DHYDROGEOCHEM of subsurface flow, heat transfer, and reactive chemical transport in which a complete suite of geochemical reactions was taken into account. Strong coupling was used for steady-state simulations while weak coupling employed for transient simulations to save computation time. The flow equation was discretized using Galerkin FEM while both the reactive chemical transport and heat transfer equations were solved by either the hybrid Lagrangian-Eulerian or the conventional Eulerian FEM. The Picard method was used to handle nonlinearity in both reactive chemical transport and unsaturated subsurface flow, while the Newton-Raphson method was employed to calculate chemical equilibrium. Barry et al. [11] compared operator-splitting methods for solving coupled chemical nonequilibrium reaction/groundwater transport models. Cirpka et al. [9, 10] and Yeh et al. [16] used both the SIA two-step method and the operator-splitting two-step method for comparisons.

Although many aspects of such reactive transport processes have been investigated experimentally and numerically, most of the existing studies are tailed to a specific problem and our understanding on the coupled process is far from satisfactory. Practical problems arising from engineering requirements of groundwater recharge with wastewater and its possible application in field scale motivate us to gain insight into mechanism/kinetics of fundamental reaction processes and interactions among them to carry out further investigations on phenomena of reactive transport, and to develop numerical models for solving the coupled problems efficiently. In the present study, a new solution procedure based on operator-splitting techniques is developed to solve the coupled model of convection-dispersion with reactive transport. Section 2 states the mathematical models describing the saturated groundwater flow and contaminant transport in detail. For dealing with the varying water table, the finite volume method based on volume integral of partial differential equations (PDEs) is introduced. Section 3 describes the application of operator-splitting techniques to the reactive transport model. An implicit upwind finite difference scheme is established to solve the convection-dominant equation. Numerical solutions are calculated in an explicit way by taking advantages of the potential flow field. The dispersion flux term is discretized in a standard central-difference scheme. Either the PJCG method or the Thomas method based on local one-dimensional fashion is appropriate for solving the dispersion-dominant equation. Applications of the numerical solution method developed in the present study to the demonstration project of groundwater recharge with secondary effluent at Gaobeidian STP are described in Section 4. The details of this demonstration project are referred to [20].

2. Mathematical Models

The physical model of groundwater recharge is illustrated in Figure 1. The basic assumptions for groundwater flow in porous media are (a) Darcy’s law, impling that the flow velocity is proportional to the pressure gradient; (b) nondeformable skeleton of the porous media; (c) Bossinesq approximation for the density of groundwater; (d) isothermal and local equilibrium of chemical reactions, biodegradation; (e) immobile solid phases, Fickian dispersion and ideal mixing.

Figure 1 

Schematic diagram of groundwater recharge.

2.1. The Partial Differential Equations—PDEs

2.1.1. Governing Equations of Saturated Groundwater Flow

The governing equation for saturated groundwater flow model can be expressed as𝜇_𝑆𝜕𝐻𝜕𝑡+∇•[𝐤•∇𝐻]=𝑞_𝑊, (1)where 𝜇_𝑆=𝜙𝛽_𝐻+𝛼_𝐻 is the specific storage coefficient, 𝛽_𝐻 and 𝛼_𝐻 are the compressibility of groundwater and the compressibility of the skeleton, respectively; in a Cartesian coordinate system 𝐻=𝑃/(𝜌_𝑊𝑔)+𝑧 is the water head; k is the hydraulic conductivity tensor associated with 𝐤^𝐏 by 𝐤=(𝜌_𝑊𝑔)𝐤^𝑃/𝜇, where 𝐤^𝐏 is the intrinsic permeability tensor, μ is the viscosity of groundwater, g is the gravity acceleration, z is the vertical coordinate pointing upward from the datum below the water table. 𝑞_𝑊 is the source term in association with recharge or production.

Usually, by selecting the coordinate system to coincide with the main seepage direction, the hydraulic conductivity tensor can be simplified as𝑘𝐤=^xx000𝑘_yy000𝑘_zz. (2)

With the solution of the water head (1), the velocity field can be calculated using Darcy’s law𝐔=−𝐤•∇𝐻, (3)where 𝐔=(𝑢,𝑣,𝑤)^𝑇.

2.1.2. Governing Equations of Reactive Transport

The transport processes of solute contaminants simultaneously occur with the groundwater movement, which include convection-dispersion, soil adsorption/desorption, and/or biodegradation. For the component m of soluble matters, the reactive transport model can be expressed by the following equation:𝜕(𝜙𝐶_𝑚)𝜕𝑡+∇•[(𝐶_𝑚𝐔−𝐅_𝑚)]=𝐶^𝑊_𝑚𝑞_𝑊+Γ_ads+Γ_bio. (4)In above equation, the parameter 𝜙 is the porosity of rock/soil; 𝐶_𝑚 is the concentration of solute matter 𝑚;𝐶^𝑊_𝑚 is the concentration in the injection/production. 𝐅_𝑚 represents the mass flux of component m due to dispersion and can be expressed as𝐅_𝑚=𝜙𝐃_𝑚∇𝐶_𝑚, (5)in which 𝐃_𝑚 is the dispersion tensor (including molecular diffusion) and can be calculated as𝜙𝑑_𝛼𝛽𝑑=𝜙^𝑚𝜏𝛿_𝛼𝛽+(𝛼_𝐿−𝛼_𝑇)𝑢_𝛼𝑢_𝛽|𝑈|+𝛼_𝑇|𝑈|𝛿_𝛼𝛽(𝛼,𝛽=𝑥,𝑦,𝑧), (6) where 𝑑_𝑚 refers to molecular diffusion; 𝜏 is the tortuous factor; 𝛼_𝐿 and 𝛼_𝑇 are longitudinal/tangential dispersion coefficients; 𝛿_𝛼𝛽=1 (for 𝛼=𝛽) or = 0 (for 𝛼≠𝛽) (𝛼,𝛽=𝑥,𝑦,𝑧). Usually, the molecular diffusion is small enough to be safely excluded from the analysis ([21, 22]).

The term Γ_ads on the right-hand side of (4) represents the concentration change due to soil adsorption, which can be expressed in terms of aqueous concentrationΓ_ads=−𝜌_𝑏𝜕𝐶^𝑆_𝑚𝜕𝑡=−𝜌_𝑏𝐾^𝑑_𝑚𝜕𝐶_𝑚𝜌𝜕𝑡=−^𝑏𝐾^𝑑_𝑚𝜙𝜕𝜙𝐶_𝑚,𝜕𝑡 (7)where 𝜌_𝑏 is the soil dry bulk density; 𝐶^𝑆_𝑚 is the contaminant concentration of solid-phase adsorbed to the particle surface; and 𝐾^𝑑_𝑚 is the linear distribution coefficient 𝐶^𝑆_𝑚=𝐾^𝑑_𝑚𝐶_𝑚.

The biodegradation term can be written as Γ_bio=𝜕𝐶_𝑚𝜕𝑡=−𝑘^𝑖_𝑚(𝐶_𝑚−𝐶^𝑈_𝑚)(𝑖=1,2,3), (8) where 𝑘^𝑖_𝑚 is degradation coefficient, the superscripts 1, 2, and 3 represent aerobic mode, anoxic mode, or anaerobic mode, respectively. 𝐶^𝑈_𝑚 is called the undegradable concentration. For an organic contaminant, the biodegradation in specific time duration can be identified as one of the three degradation modes. Thus the superscript i will be neglected and practical values will be assigned to the degradation coefficient 𝑘_𝑚 according to the biodegradation process, aerobic, anoxic, or anaerobic.

By substituting (7) and (8) into (4), we have𝑅_𝑚𝜕(𝜙𝐶_𝑚)𝜕𝑡+∇•[(𝐶_𝑚𝐔−𝐅_𝑚)]+𝑘_𝑚𝐶_𝑚=Σ_𝑚, (9)where 𝑅_𝑚=1+(𝜌_𝑏𝐾^𝑑_𝑚)/𝜙 is called the retardation coefficient. For the nondegradable components, we have 𝑘_𝑚=0 while for nonadsorbable component, 𝑅_𝑚=1. Σ_𝑚=𝐶^𝑊_𝑚𝑞_𝑊+𝑘_𝑚𝐶^𝑈_𝑚 is new source term due to both recharge/production and biodegradation.

2.1.3. Initial and Boundary Conditions

For the particular problem that will be encountered in the recharge site, the initial gradient of water head can be neglected while the solute contaminants have a mean initial distribution 𝐻(𝑥,𝑦,𝑡)|_𝑡=0=𝐻₀(𝑥,𝑦)=𝐻₀𝐶,(𝑥,𝑦)∈Ω,_𝑚(𝑥,𝑦,𝑧,𝑡)|_𝑡=0=𝐶_𝑚0(𝑥,𝑦,𝑧)=𝐶_𝑚0,(𝑥,𝑦,𝑧)∈𝐺. (10a)(10b) where 𝑟<𝑟₁<𝑟₂<1.

The finite difference equations (29a), (29b), and (29c) can be expressed as 𝐵𝐶¹⋅𝐶^𝑛+𝑟1_𝑚+𝐵𝑊⋅𝐶^𝑛+𝑟1_{𝑚,𝑖−1}+𝐵𝐸⋅𝐶^𝑛+𝑟1_𝑚,𝑖+1=𝑅¹_𝑠,𝐵𝐶²⋅𝐶^𝑛+𝑟2_𝑚+𝐵𝑁⋅𝐶^𝑛+𝑟2_{𝑚,𝑗−1}+𝐵𝑆⋅𝐶^𝑛+𝑟2_𝑚,𝑗+1=𝑅²_𝑠,𝐵𝐶³⋅𝐶_𝑚^𝑛+1+𝐵𝑇⋅𝐶^𝑛+1_{𝑚,𝑘−1}+𝐵𝐵⋅𝐶^𝑛+1_𝑚,𝑘+1=𝑅³_𝑠, (30a)(30b)(30c)where𝐵𝑊=2(𝐴𝜙𝑑_𝑥𝑥)_𝑤(Δ𝑥₊₋₁+Δ𝑥_𝑖),𝐵𝐸=2(𝐴𝜙𝑑_𝑥𝑥)_𝑒(Δ𝑥_𝑖+1+Δ𝑥_𝑖),𝐵𝐶¹=𝜙𝑅_𝑚𝑉𝑅Δ𝑡−𝐵𝑊−𝐵𝐸,_1𝑠=𝜙𝑅_𝑚𝑉Δ𝑡(𝐶_𝑚)^𝑛+𝑟+[𝜙𝐴(𝑑_𝑥𝑦𝜕𝐶_𝑚𝜕𝑦+𝑑_𝑥𝑧𝜕𝐶_𝑚𝜕𝑧)]_𝑒^𝑛+𝑟−[𝜙𝐴(𝑑_𝑥𝑦𝜕𝐶_𝑚𝜕𝑦+𝑑_𝑥𝑧𝜕𝐶_𝑚𝜕𝑧)]_𝑤^𝑛+𝑟;𝐵𝑁=2(𝐴𝜙𝑑_𝑦𝑦)_𝑛(Δ𝑦_𝑗−1+Δ𝑦_𝑗),𝐵𝑆=2(𝐴𝜙𝑑_𝑦𝑦)_𝑠(Δ𝑦_𝑗+1+Δ𝑦_𝑗),𝐵𝐶²=𝜙𝑅_𝑚𝑉𝑅Δ𝑡−𝐵𝑁−𝐵𝑆,_2𝑠=𝜙𝑅_𝑚𝑉Δ𝑡(𝐶_𝑚)^𝑛+𝑟1+[𝜙𝐴(𝑑_𝑦𝑥𝜕𝐶_𝑚𝜕𝑥+𝑑_𝑦𝑧𝜕𝐶_𝑚𝜕𝑧)]^𝑛+𝑟1_𝑠−[𝜙𝐴(𝑑_𝑦𝑥𝜕𝐶_𝑚𝜕𝑥+𝑑_𝑦𝑧𝜕𝐶_𝑚𝜕𝑧)]^𝑛+𝑟1_𝑛,𝐵𝑇=2(𝐴𝜙𝑑_𝑧𝑧)_𝑡(Δ𝑧_𝑘−1+Δ𝑧_𝑘),𝐵𝐵=2(𝐴𝜙𝑑_𝑧𝑧)_𝑏(Δ𝑧_𝑘+1+Δ𝑧_𝑘),𝐵𝐶³=𝜙𝑅_𝑚𝑉𝑅Δ𝑡−𝐵𝑇−𝐵𝐵,_3𝑠=𝜙𝑅_𝑚𝑉Δ𝑡(𝐶_𝑚)^𝑛+𝑟2+[𝜙𝐴(𝑑_𝑧𝑧𝜕𝐶_𝑚𝜕𝑥+𝑑_𝑧𝑦𝜕𝐶_𝑚𝜕𝑦)]^𝑛+𝑟2_𝑏−[𝜙𝐴(𝑑_𝑧𝑧𝜕𝐶_𝑚𝜕𝑥+𝑑_𝑧𝑦𝜕𝐶_𝑚𝜕𝑦)]^𝑛+𝑟2_𝑡. (31) The solution of (30c) is the final solution of the reactive transport model: PDE (4) or FVM form (13b). To demonstrate the validity of the present method, we also computed the explicit solution of (13b). A target field of 250 m (NX=25)∗ 110 m (NY=11)∗ 20 m (NZ = 10) was selected to simulate the distribution of DO concentration during 180 days. The explicit solution with a very small time step of 0.1 day can be taken as the approximate of the real solution. Comparisons between the present implicit solution and the explicit solution show a good agreement.

4. Applications

The present solution procedure is a part of the numerical simulation work, which has been carried out for the demonstration project of groundwater recharge with secondary effluent at Gaobeidian STP [20]. The target site is on the Beijing alluvial plain, where there is a tight clay layer at 16 m below the surface, which can be set as the datum. The groundwater level is about 9.7 m. Three recharge pools were designed for alternative operations (Figure 6). One production well and three observation wells were designed to monitor groundwater movement and contaminant transport. The recharge implementation shows a very low recharge rate (20–30 t/day) for the recharge mode of recharge pool. With such a recharge rate, nothing can be monitored at the observation wells, nor can it be distinguished from the simulation results. The test values of hydraulic conductivity explain this phenomenon: 0.12 m/day (for 0–4 m depth); 0.48 m/day (for 4–8 m depth); 1.92 m /day (for 9–12 m depth); 3.84 m/day (for 13–16 m depth).

Figure 6 

Photo of the recharge pool.

Based on these practices, we suggest a recharge mode of injection well to recharge water directly into layers with high hydraulic conductivity and to give simulation results. The numerical model was established on a grid system 50×30×10 for the simulation domain of 1000 m × 600 m × 20 m. The initial conditions are groundwater level, 9.7 m deep; DOC concentration, 2.4 mg/L; oxygen concentration, 0.01 mg/L. The undegradable concentration is considered to be 1.2 mg/L. According to laboratory tests, the distribution coefficient for adsorption is 0.00012 [l Kg/L]. The degradation coefficients for aerobic environment (DO > 0.9 mg/L), anoxic environment (0.9 mg/L > DO > 0.01 mg/L), and anaerobic environment (0.01 mg/L > DO) are 0.139, 0.07, and 0.02, respectively. Hydraulic and reactive transport processes are simulated while some parameters are estimated, including hydraulic-steady-state; hydraulic-penetration; maximum-recharge; reused-water-quality. The results are partly presented in this paper.

The quality of groundwater is a determinative factor for the possible reuse of secondary effluent through subsurface recharge. Here, we select the DOC as the representative contaminant to assess the groundwater quality. Figure 7 shows the DOC distribution at 180 days with a recharge rate of 200 t/day and a production rate of 200 t/day.

Figure 7 

DOC distribution at 180 days with 200 t/day for both recharge and production.

It can be seen that DOC concentration becomes higher with time increasing, especially for a well closer to the drainage site. The contaminants reach the production well without sufficient biodegradation. For a well 300 m away, the increase of DOC concentration can be neglected, even after 720 days recharge. The water quality at the production well can be improved by enlarging its distance from the drainage site and/or by increasing DO concentration to enhance the biodegradation. The effect of DO concentration on the quality of water at production well 150 m away from the recharge site is shown in Table 1.

Table 1 

Variation of DOC with two different DO concentrations for a well of 150 m.

5. Conclusions and Discussions

Reuse of wastewater through subsurface recharge is an effective mode of water resource management. The fluid flow due to wastewater recharge always accompanies with reactive transport of solute contaminants. To avoid extraneous pollution to groundwater by such recharge, numerical simulation must be undertaken to assess the impact on the groundwater quality and the subsurface environments before the recharge implementation. The coupling of groundwater movement and reactive transport leads to complicated nonlinear models, which have been proved difficult to obtain numerical solution. Based on the operator-splitting techniques, the present study develops a new solution procedure to solve the reactive transport model of saturated flow during the shallow groundwater recharge with wastewater. For dealing with the varying water table, the finite volume integral of partial differential equations (PDEs) is introduced. By applying the operator-splitting techniques, the reactive transport model is decomposed into one convection-dominant equation plus a dispersion-dominant equation. An implicit upwind finite difference scheme is used to solve the convection-dominant equation. Numerical solutions are calculated in an explicit way by taking advantages of the potential flow field. A standard central-difference scheme is applied to the dispersion-dominant equation. Either the PJCG method or the Thomas method based on local one-dimensional is suitable for solving the dispersion-dominant equation. The explicit solution with a very small time step, which can be taken as the approximate of the real solution, is also computed to demonstrate the validity of the present method.

To simulate 2D convection-dominant reactive transport process, Cirpka et al. [9] developed a new method of streamline-oriented grids based on cell-centered finite volume scheme. In their method, the transport scheme is an adaption of the principal direction, thus the transverse flux diminishes. After the computation of the hydraulic model based on rectangular grids, the streamline-oriented grids containing quadrilateral cells should be generated, for which spatial gradient along the direction of the streamlines and orthogonal are required. The resulting model is solved by an operator-splitting approach, in which the convective transport is solved by a nonlinear explicit method while the dispersive transport is solved implicitly. Though this method shows its advantages (11a), (11b), it can only be used to two-dimensional problems because of the use of streamlines. For transient problems, the generation of stream-oriented grids for each time step will be time consuming. Furthermore, only the explicit solution can be obtained for the convective transport, which is always the dominant process in the reactive transport. Compared with the method of Cirpka et al. [9], the present solution method is not confined to the two-dimensional problems. It can be freely applied to three-dimensional problems of transient flow.

It should be mentioned that the solution procedure of high efficiency described in this study is functional for the reactive transport due to secondary effluent recharge under which the concentration of solute contaminants is relatively low. The concentration change due to adsorption and biodegradation can be approximated to Monod types. For more common cases of coupled nonlinear reactive transport problems, Simpson and Landman [23] developed split operator methods applied to reactive transport with Monod kinetics and analyzed the removal of errors from both the boundary and internal errors.