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 aswhere 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

With the solution of the water head (1), the velocity field can be calculated using Darcy’s lawwhere .

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: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 asin which is the dispersion tensor (including molecular diffusion) and can be calculated as where refers to molecular diffusion; is the tortuous factor; and are longitudinal/tangential dispersion coefficients; (for ) or = 0 (for ) (). Usually, the molecular diffusion is small enough to be safely excluded from the analysis ([21, 22]).

The term on the right-hand side of (4) represents the concentration change due to soil adsorption, which can be expressed in terms of aqueous concentrationwhere 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 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 havewhere is called the retardation coefficient. For the nondegradable components, we have while for nonadsorbable component, . 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

Dirichlet boundary conditions are assigned to both the hydraulic model and the contaminants transport model

2.2. Finite Volume Method for PDEs

Because the equations describing reactive processes and biodegradations are based on concentrations, the most common way to solve these equations for spatially variable domains is to divide the domain into control volumes and calculate reactive processes independently in each of them. For the unconfined groundwater flow, the water table varies with the recharge and/or production, intersects some grids, and gradually changes its vertical level. Thus the FVM, based on the integral form of PDE’s with respect to grid volume, is employed for both the hydraulic model and the reactive transport model. To establish the numerical model, the region to be simulated is divided into an central-valued grid system.

The integral form of the (1) isThe water head (12) can be solved using the preconditioning Jacobi conjugate gradient method or several other existing methods. The detailed discussion on the simulation model can be referred to [20]. Hereafter, we will concentrate our interests on the solution procedure of the reactive transport model (9).

Integrating the PDE (9) with respect to grid volume yields

By applying the volume integral for accumulated term and the Gauss formula for the flux term, where n is the unit vector in the normal direction of S, (13a) becomes

There exist three types of grids: (a) fully submerged grids below the water table; (b) partly submerged grids crossing water table; (c) empty grids above water table. For the saturated flow model, type (c) grids can be neglected because they do not involve groundwater movement. We will only deal with the former two types of grids and focus our discussion on the type (b) grids.

For type (a) grids, the volume integral can be expressed as where , , and are dimensions of grid i in x, y, z direction, respectively; is the geometry volume of grid i. Meanwhile, the circular integral can be expanded as where the vector represents either or are the section areas corresponding to west, east, north, south, top, and bottom faces of the target grid, respectively.

It will be somehow complicated to deal with type (b) grids because of the varying water table. Figure 2 illustrates the groundwater flow pattern for such grids.

Figure 2: Groundwater flow for type (b) grids.

Because only part of such grid contains groundwater and involves the fluid flux, the volume integral of accumulation term can be expanded as where is the effective volume involving the groundwater movement and is the height from grid bottom to the water table. It is clear that should be kept. Correspondingly, the enclosure surface integral of the flux term can be expanded as In above equation are effective section areas for west, east, north, and south faces of the target grid under the water table and should be calculated as where the subscript means the interfacial position between two adjacent grids.

3. Operator-Splitting Techniques

In the present study, a peculiar operator-splitting method is derived to achieve high-efficiency solution.

The temporal differential term can be discretized in a forward finite difference scheme

By splitting this term into two parts, (13b) is decomposed into two equations with different characteristics: one convection-dominated equation plus a dispersion-dominant equation where . The convection-dominant equation (18a) should be solved at the first stage. With the numerical solution of (18a) as the initial condition, the dispersion-dominant equation (18b) can be solved.

3.1. Implicit Upwind Scheme of Convection-Dominant Equation

To solve the convection-dominant equation (18a) efficiently, we derived a special solution procedure by using an implicit upwind scheme. As mentioned above, there exist two different types of grids: type (a) and type (b). Correspondingly, both the expression of accumulated integral term and the treatment of enclosure surface integral term are slightly different for the two types of grids. We will discuss the discretization procedure for each type of grids separately.

3.1.1. Fully Submerged Grids

Using implicit upwind finite difference scheme, the convection term for a fully submerged grid can be expanded as follows (for grid ): The last six terms are treated explicitly in the present study.

3.1.2. Partially Submerged Grids

For a partially submerged grid, the expansion of the bounding surface integral term may be somewhat different from that for a fully submerged grid because no flux exists on the top boundary. The section areas of west, east, north, and south boundaries should be replaced by the effective areas, as illustrated in (15b),

To make the FVM expansion of integral term identical for both the fully submerged grids and partially submerged grids, we also express (19b) as the same scheme as (19a) with exceptions that , the volume integral term is replaced by (15a) while the effective areas, are used instead of the section areas In the following section, we will use the universal form to discuss the solution procedure.

3.1.3. Solution Procedure

Substituting (14a) and (19a) into the convection-dominant equation (18a), we have (for grid ) where

The present study employs an efficient solution procedure to compute the implicit solution in an explicit way by taking advantages of potential characteristics of the flow field and the implicit upwind scheme. In a potential field of discrete grids, there must exist at least one grid from which no influx (except for the source) occurs at each boundary. For such a grid, we can see from the finite difference equation (20) that the upwind scheme results in zero values of coefficients AW, AE, AN, AS, AT, and AB. The solution of can be independently derived:

We identify and indicate such a special grid as the starting grid. For the adjacent grids of the starting grid, we can also find at least one grid whose adjacent grids except for the starting grid are NOT its upwind grids. By assuming that this grid is denoted as , its solution can be calculated as

Thus we can define an index array to record the abovementioned upwind-downwind sequence of the grids l so that all the values of can be explicitly calculated one by one along the index array . That is to say, as the solution procedure is implemented from 1 to , the value of grid is already known or the grid is NOT the upwind of .

Now, let us discuss how to determine the index array . Without loss of the generality, we define six auxiliary arrays, , to indicate the states of west, east, south, north, top, and bottom boundaries for a grid . For each of the six boundaries, say the west, we assign for influx and for other states, as illustrated in Figure 3.

Figure 3: Determination of the boundary states.

After assigning values of the six auxiliary arrays for each grid, we can identify the starting grid in the following procedure.

For a grid, say , if the following criterion is met: we call such grid as the starting grid and indicate it as , as shown in Figure 4.

Figure 4: Illustration of the starting grid.

The solution for this grid can be calculated with (22). For its adjacent grids, the starting grid can be taken as the upwind grid and its solution is known. Therefore, we reassign boundary state to be 0 for the interfacial boundaries between the starting grid and its adjacent grids,

Then we examine boundary states for adjacent grids of the starting grid. If any, say , meets the condition (24), we indicate that , as illustrated in Figure 5. Following this procedure, we can indicate all the grids and determine the array . We then calculate the solution for all the grids indexed by the array in an explicit way.

Figure 5: Schematic diagram of the upwind-downwind grids.

3.2. Solution of the Dispersion-Dominated Equation

With the solution of convection-dominant equation as the initial condition, the dispersion-dominated equation (18b) can be solved by either the PJCG method or Thomas method based on local-one-dimensional scheme.

3.2.1. Central-Difference Scheme and PJCG Method

In the present study, the dispersion-dominant term was also approximated by standard centered difference scheme and the equation be solved using the PJCG method. The details of standard central-difference scheme were given in Russo et al. [17, 18] while the description on PJCG method can be found in many references. Thus we will concentrate our discussion on the treatment of local-one-dimensional scheme.

3.2.2. Thomas Method Based on Local-One-Dimensional Scheme

For most models of groundwater problem, nondiagonal terms of the hydraulic conductivity tensor are neglected by selecting the coordinate system to coincide with the main seepage direction while nondiagonal terms of the diffusion tensor are neglected. In the present study, we split the diffusion flux into two parts corresponding to diagonal terms and nondiagonal terms of the dispersion tensor, respectively. The nondiagonal terms are treated as explicit termwhere Applying (14a), (14b), and (26) to (18b) leads to With the solution of , now, we begin to solve (28). By splitting (28) into three local one-dimensional finite difference scheme in x, y, and z, respectively, we can use Thomas method, which is the most efficient method for solving a tridiagonal matrix system as long as the diagonal elements are dominant, to obtain the solution where .

The finite difference equations (29a), (29b), and (29c) can be expressed as where 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 110 m 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 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.