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

Authors: Shutang Zhu; Gang Zhao; Yuanle Ma
Journal: Modelling and Simulation in Engineering
 (2008)

---

## 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.


---

## Body

## 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. Section2 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 Figure1. 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)μS∂H∂t+∇•[k•∇H]=qW,where μS=ϕβH+αH is the specific storage coefficient, βH and αH are the compressibility of groundwater and the
compressibility of the skeleton, respectively; in a Cartesian coordinate system H=P/(ρWg)+z is the water head; k is the hydraulic
conductivity tensor associated with kP
by k=(ρWg)kP/μ,
where kP  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. qW 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(2)k=[kxx000kyy000kzz].With the solution of the water head
(1), the
velocity field can be calculated using Darcy’s law(3)U=−k•∇H,where U=(u,v,w)T.


#### 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 componentm of soluble
matters, the reactive transport model can be expressed by the  following
equation:(4)∂(ϕCm)∂t+∇•[(CmU−Fm)]=CmWqW+Γads+Γbio.In above equation, the parameter ϕ is the porosity of rock/soil; Cm is the concentration  of solute matter m;CmW  is the concentration in the
injection/production. Fm
represents the mass flux of
component m due to dispersion and can be expressed as(5)Fm=ϕDm∇Cm,in which Dm is the dispersion tensor (including molecular
diffusion) and can be calculated as(6)ϕdαβ=ϕdmτδαβ+(αL−αT)uαuβ|U|+αT|U|δαβ(α,β=x,y,z), where dm refers to molecular diffusion; τ is the tortuous factor; αL and αT are longitudinal/tangential dispersion
coefficients; δαβ=1 (for α=β) or = 0 (for α≠β) (α,β=x,y,z).
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(7)Γads=−ρb∂CmS∂t=−ρbKmd∂Cm∂t=−ρbKmdϕ∂ϕCm∂t,where ρb is the soil dry bulk density; CmS is the contaminant concentration of
solid-phase adsorbed to the particle surface; and Kmd is the linear distribution coefficient CmS=KmdCm.The biodegradation term can be written as(8)Γbio=∂Cm∂t=−kmi(Cm−CmU)(i=1,2,3), where kmi is degradation coefficient, the superscripts 1,
2, and 3 represent aerobic mode, anoxic mode, or anaerobic mode, respectively. CmU 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 km according to the biodegradation process,
aerobic, anoxic, or anaerobic.By substituting (7) and
(8) into (4), we have(9)Rm∂(ϕCm)∂t+∇•[(CmU−Fm)]+kmCm=Σm,where Rm=1+(ρbKmd)/ϕ is called the retardation coefficient. For the
nondegradable components, we have km=0 while for nonadsorbable component, Rm=1. Σm=CmWqW+kmCmU 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(10a)H(x,y,t)|t=0=H0(x,y)=H0,(x,y)∈Ω,(10b)Cm(x,y,z,t)|t=0=Cm0(x,y,z)=Cm0,(x,y,z)∈G.Dirichlet boundary conditions are assigned to both
the hydraulic model and the contaminants transport model(11a)H(x,y,t)|Γ1=H1(x,y,t),(x,y)∈Γ1,t>0,(11b)C(x,y,z,t)|Γ1=C1(x,y,z,t),(x,y)∈Γ1,t>0.


### 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 anNX∗NY∗NZ central-valued grid system.The integral form of the (1) is(12)∭ViμS∂H∂tdV+∭Vi∇•[k•∇H]dV=∭ViqWdV.The 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
(13a)∭ViRm∂(ϕCm)∂tdV+∭Vi∇•[(CmU−Fm)]dV+∭VikmCmdV=∭ViΣmdV.By applying the volume integral∭ViAdV=A∭VidV for accumulated term and the Gauss formula ∭Vi∇•VdV=∯SiV•ndS for the flux term, where n is the unit
vector in the normal direction of S, (13a) becomes(13b)Rm∂(ϕCm)∂t∭VidV+∯S(CmU)•ndS−∯SFm•ndS+kmCm∭VidV=Σm∭VidV.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(14a)∭VidV=ΔxiΔyiΔzi=Vi,where Δxi, Δyi,
and Δzi are dimensions of grid i in x, y, z direction, respectively; Vi is the geometry volume of grid i. Meanwhile, the circular integral can be
expanded as(14b)∯Si[V]•ndS=−[Vx]wAw+[Vx]eAe−[Vy]nAn+[Vy]sAs−[Vz]tAt+[Vz]bAb,
where the vector V=[Vx,Vy,Vz]T represents either [CmU] or [Fm];Aη(η=w,e,n,s,t,b)   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.
Figure2 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(15a)∭VidV=Vie=ΔxiΔyihi,where Vie=ΔxiΔyihi is the effective volume involving the
groundwater movement and hi is the height from grid bottom to the water
table. It is clear that hi/Δzi≤1 should be kept. Correspondingly, the enclosure surface
integral of the flux term can be expanded as(15b)∯Si[V]•ndS=−[Vx]wA^w+[Vx]eA^e−[Vy]nA^n+[Vy]sA^s+[Vz]bAb.
In above equation A^η(η=w,e,n,s)  are effective section areas for west, east,
north, and south faces of the target grid under the water table and should be
calculated as
(16)A^w=Δyhi−1/2,A^e=Δyhi+1/2,A^n=Δxhj−1/2,A^s=Δxhj+1/2, where the subscript ±1/2 means the interfacial position between two
adjacent grids.


## 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)μS∂H∂t+∇•[k•∇H]=qW,where μS=ϕβH+αH is the specific storage coefficient, βH and αH are the compressibility of groundwater and the
compressibility of the skeleton, respectively; in a Cartesian coordinate system H=P/(ρWg)+z is the water head; k is the hydraulic
conductivity tensor associated with kP
by k=(ρWg)kP/μ,
where kP  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. qW 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(2)k=[kxx000kyy000kzz].With the solution of the water head
(1), the
velocity field can be calculated using Darcy’s law(3)U=−k•∇H,where U=(u,v,w)T.


### 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 componentm of soluble
matters, the reactive transport model can be expressed by the  following
equation:(4)∂(ϕCm)∂t+∇•[(CmU−Fm)]=CmWqW+Γads+Γbio.In above equation, the parameter ϕ is the porosity of rock/soil; Cm is the concentration  of solute matter m;CmW  is the concentration in the
injection/production. Fm
represents the mass flux of
component m due to dispersion and can be expressed as(5)Fm=ϕDm∇Cm,in which Dm is the dispersion tensor (including molecular
diffusion) and can be calculated as(6)ϕdαβ=ϕdmτδαβ+(αL−αT)uαuβ|U|+αT|U|δαβ(α,β=x,y,z), where dm refers to molecular diffusion; τ is the tortuous factor; αL and αT are longitudinal/tangential dispersion
coefficients; δαβ=1 (for α=β) or = 0 (for α≠β) (α,β=x,y,z).
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(7)Γads=−ρb∂CmS∂t=−ρbKmd∂Cm∂t=−ρbKmdϕ∂ϕCm∂t,where ρb is the soil dry bulk density; CmS is the contaminant concentration of
solid-phase adsorbed to the particle surface; and Kmd is the linear distribution coefficient CmS=KmdCm.The biodegradation term can be written as(8)Γbio=∂Cm∂t=−kmi(Cm−CmU)(i=1,2,3), where kmi is degradation coefficient, the superscripts 1,
2, and 3 represent aerobic mode, anoxic mode, or anaerobic mode, respectively. CmU 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 km according to the biodegradation process,
aerobic, anoxic, or anaerobic.By substituting (7) and
(8) into (4), we have(9)Rm∂(ϕCm)∂t+∇•[(CmU−Fm)]+kmCm=Σm,where Rm=1+(ρbKmd)/ϕ is called the retardation coefficient. For the
nondegradable components, we have km=0 while for nonadsorbable component, Rm=1. Σm=CmWqW+kmCmU 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(10a)H(x,y,t)|t=0=H0(x,y)=H0,(x,y)∈Ω,(10b)Cm(x,y,z,t)|t=0=Cm0(x,y,z)=Cm0,(x,y,z)∈G.Dirichlet boundary conditions are assigned to both
the hydraulic model and the contaminants transport model(11a)H(x,y,t)|Γ1=H1(x,y,t),(x,y)∈Γ1,t>0,(11b)C(x,y,z,t)|Γ1=C1(x,y,z,t),(x,y)∈Γ1,t>0.


## 2.1.1. Governing Equations of Saturated  Groundwater Flow

The governing equation for saturated groundwater flow
model can be expressed as(1)μS∂H∂t+∇•[k•∇H]=qW,where μS=ϕβH+αH is the specific storage coefficient, βH and αH are the compressibility of groundwater and the
compressibility of the skeleton, respectively; in a Cartesian coordinate system H=P/(ρWg)+z is the water head; k is the hydraulic
conductivity tensor associated with kP
by k=(ρWg)kP/μ,
where kP  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. qW 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(2)k=[kxx000kyy000kzz].With the solution of the water head
(1), the
velocity field can be calculated using Darcy’s law(3)U=−k•∇H,where U=(u,v,w)T.


## 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 componentm of soluble
matters, the reactive transport model can be expressed by the  following
equation:(4)∂(ϕCm)∂t+∇•[(CmU−Fm)]=CmWqW+Γads+Γbio.In above equation, the parameter ϕ is the porosity of rock/soil; Cm is the concentration  of solute matter m;CmW  is the concentration in the
injection/production. Fm
represents the mass flux of
component m due to dispersion and can be expressed as(5)Fm=ϕDm∇Cm,in which Dm is the dispersion tensor (including molecular
diffusion) and can be calculated as(6)ϕdαβ=ϕdmτδαβ+(αL−αT)uαuβ|U|+αT|U|δαβ(α,β=x,y,z), where dm refers to molecular diffusion; τ is the tortuous factor; αL and αT are longitudinal/tangential dispersion
coefficients; δαβ=1 (for α=β) or = 0 (for α≠β) (α,β=x,y,z).
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(7)Γads=−ρb∂CmS∂t=−ρbKmd∂Cm∂t=−ρbKmdϕ∂ϕCm∂t,where ρb is the soil dry bulk density; CmS is the contaminant concentration of
solid-phase adsorbed to the particle surface; and Kmd is the linear distribution coefficient CmS=KmdCm.The biodegradation term can be written as(8)Γbio=∂Cm∂t=−kmi(Cm−CmU)(i=1,2,3), where kmi is degradation coefficient, the superscripts 1,
2, and 3 represent aerobic mode, anoxic mode, or anaerobic mode, respectively. CmU 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 km according to the biodegradation process,
aerobic, anoxic, or anaerobic.By substituting (7) and
(8) into (4), we have(9)Rm∂(ϕCm)∂t+∇•[(CmU−Fm)]+kmCm=Σm,where Rm=1+(ρbKmd)/ϕ is called the retardation coefficient. For the
nondegradable components, we have km=0 while for nonadsorbable component, Rm=1. Σm=CmWqW+kmCmU 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(10a)H(x,y,t)|t=0=H0(x,y)=H0,(x,y)∈Ω,(10b)Cm(x,y,z,t)|t=0=Cm0(x,y,z)=Cm0,(x,y,z)∈G.Dirichlet boundary conditions are assigned to both
the hydraulic model and the contaminants transport model(11a)H(x,y,t)|Γ1=H1(x,y,t),(x,y)∈Γ1,t>0,(11b)C(x,y,z,t)|Γ1=C1(x,y,z,t),(x,y)∈Γ1,t>0.


## 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 anNX∗NY∗NZ central-valued grid system.The integral form of the (1) is(12)∭ViμS∂H∂tdV+∭Vi∇•[k•∇H]dV=∭ViqWdV.The 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
(13a)∭ViRm∂(ϕCm)∂tdV+∭Vi∇•[(CmU−Fm)]dV+∭VikmCmdV=∭ViΣmdV.By applying the volume integral∭ViAdV=A∭VidV for accumulated term and the Gauss formula ∭Vi∇•VdV=∯SiV•ndS for the flux term, where n is the unit
vector in the normal direction of S, (13a) becomes(13b)Rm∂(ϕCm)∂t∭VidV+∯S(CmU)•ndS−∯SFm•ndS+kmCm∭VidV=Σm∭VidV.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(14a)∭VidV=ΔxiΔyiΔzi=Vi,where Δxi, Δyi,
and Δzi are dimensions of grid i in x, y, z direction, respectively; Vi is the geometry volume of grid i. Meanwhile, the circular integral can be
expanded as(14b)∯Si[V]•ndS=−[Vx]wAw+[Vx]eAe−[Vy]nAn+[Vy]sAs−[Vz]tAt+[Vz]bAb,
where the vector V=[Vx,Vy,Vz]T represents either [CmU] or [Fm];Aη(η=w,e,n,s,t,b)   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.
Figure2 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(15a)∭VidV=Vie=ΔxiΔyihi,where Vie=ΔxiΔyihi is the effective volume involving the
groundwater movement and hi is the height from grid bottom to the water
table. It is clear that hi/Δzi≤1 should be kept. Correspondingly, the enclosure surface
integral of the flux term can be expanded as(15b)∯Si[V]•ndS=−[Vx]wA^w+[Vx]eA^e−[Vy]nA^n+[Vy]sA^s+[Vz]bAb.
In above equation A^η(η=w,e,n,s)  are effective section areas for west, east,
north, and south faces of the target grid under the water table and should be
calculated as
(16)A^w=Δyhi−1/2,A^e=Δyhi+1/2,A^n=Δxhj−1/2,A^s=Δxhj+1/2, where the subscript ±1/2 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(17)∂(ϕCm)∂t≅(ϕCm)n+1−(ϕCm)nΔt,(Δt=tn+1−tn).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
(18a)Rm(ϕCm)n+r−(ϕCm)nΔt∭VidV+∯S(CmU)•ndS+kmCm∭VidV=Σm∭VidV,(18b)Rm(ϕCm)n+1−(ϕCm)n+rΔt∭VidV−∯SFm•ndS=0,
where 0<r<1.
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 gridijk):
(19a)∯S(CmU)•ndS=−[Cmu]wAw+[Cmu]eAe−[Cmu]nAn+[Cmu]sAs−[Cmu]tAt+[Cmu]bAb=12{(Au)e(Cm)i+1+[(Au)e−(Au)w]Cm−(Au)w(Cm)i−1}+12{(Av)s(Cm)j+1+[(Av)s−(Av)n]Cm−(Av)n(Cm)j−1}+12{(Aw)b(Cm)k+1+[(Aw)b−(Aw)t]×Cm−(Aw)t(Cm)k−1}=12((Au)e−|(Au)e|)(Cm)i+1−12((Au)w+|(Au)w|)(Cm)i−1+12((Av)s−|(Av)s|)(Cm)j+1−12((Av)n+|(Av)n|)(Cm)j−1+12((Aw)b−|(Aw)b|)(Cm)k+1−12((Aw)t+|(Aw)t|)(Cm)k−1+((Au)e−(Au)w+|(Au)e|+|(Au)w|+(Av)s−(Av)n+|(Av)s|+|(Av)n|)Cm2+((Aw)b−(Aw)t+|(Aw)b|+|(Aw)t|)Cm2+|(Au)e|(Cm)i+1−Cm2+|(Av)s|(Cm)j+1−Cm2+|(Aw)b|(Cm)k+1−Cm2−|(Au)w|Cm−(Cm)i−12−|(Av)n|Cm−(Cm)j−12−|(Aw)t|Cm−(Cm)k−12.                    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),(19b)∯S(CmU)•ndS=−[Cmu]wA^w+[Cmu]eA^e−[Cmv]nA^n+[Cmv]sA^s+[Cmw]bAb.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 vn=0,
the volume integral term is replaced by (15a) while the effective areas, A^w,A^e,A^n,A^s,
are used instead of the section areas Aw,Ae,An,As. 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 ijk) (20)AC⋅Cmn+r+AW⋅(Cm)i−1n+r+AE⋅(Cm)i+1n+r+AN⋅(Cm)j−1n+r+AS⋅(Cm)j+1n+r+AT⋅(Cm)k−1n+r+AB⋅(Cm)k+1n+r=Rs,  where(21)AW=−12((Au)w+|(Au)w|),AE=12((Au)e−|(Au)e|),AN=−12((Av)n+|(Av)n|),AS=12((Av)s−|(Av)s|),AT=−12((Aw)t+|(Aw)t|),AB=12((Aw)b−|(Aw)b|),AC=(ϕRm)ijkΔt−AWi+1−AEi−1−ANj+1−ASj−1−ATk+1−ABk−1,Rs=(ΣmV)ijk+(ϕRmCm)ijknΔt−{|(Au)e|[(Cm)i+1−Cm2]n+|(Av)s|[(Cm)j+1−Cm2]n+|(Aw)b|[(Cm)k+1−Cm2]n−|(Au)w|[Cm−(Cm)i−12]n−|(Av)n|[Cm−(Cm)j−12]n−|(Aw)t|[Cm−(Cm)k−12]n}.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 Cm can be independently derived:(22)Cmn+r=1AC[(ϕVCm)n(Δt)+Rm].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(i+1,j,k), its
solution can be calculated as(23)Cmn+ri+1,j,k=1[AC]i+1,j,k{[(ϕVkCkW)n(Δt)+Rm]i+1,j,k−[AW]i+1,j,k•Cmn+ri,j,k}.Thus we can define an index array[ID(n)=l,n=1,NX∗NY∗NZ] to record the abovementioned upwind-downwind
sequence of the grids l so that all
the values of Cmn+r(l) can be explicitly calculated one by one along
the index array ID(n). That is to say, as the solution
procedure is implemented from 1 to n+1, the value Cmn+r of grid [ID(i),i∈1,n] is already known or the grid [ID(i),i∈1,n] is NOT the upwind of l=ID(n).Now, let us discuss how to determine the index array[ID(n)=l,n=1,NX∗NY∗NZ]. Without loss of the generality, we define six
auxiliary arrays, W,E,N,S,T,B, to indicate the states of west, east,
south, north, top, and bottom boundaries for a grid ijk. For each of the six boundaries, say the west,
we assign W(ijk)=1  for influx and W(ijk)=0 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, sayijk, if the following criterion is
met: (24)W(ijk)=0,E(ijk)=0,N(ijk)=0,S(ijk=0,T(ijk)=0,B(ijk)=0, we call such grid as the starting grid and indicate it as ID(l)=ijk, 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,(25)W(ijk+1)=0,E(ijk−1)=0,N(ijk+NX)=0,S(ijk−NX)=0,T(ijk+NX∗NY)=0,B(ijk−NX∗NY)=0.Then we examine boundary states for adjacent grids
of the starting grid. If any, sayijk+1, meets the condition (24), we indicate
that ID(2)=ijk+1, as illustrated in
Figure 5. Following this procedure, we can indicate all the grids and determine
the array ID. We then calculate the
solution for all the grids indexed by the array ID 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 term(26)∯S[Fm]•n→dS=∯S[Fm]D•n→dS+∯S[Fm]N•n→dS,where(27)[Fm]D=ϕ[dxx∂Cm∂x,dyy∂Cm∂y,dzz∂Cm∂z]T,[Fm]N=ϕ[dxy∂Cm∂y+dxz∂Cm∂z,dyx∂Cm∂x+dyz∂Cm∂z,dzx∂Cm∂x+dzy∂Cm∂y]T. Applying (14a),
(14b), and
(26) to
(18b) leads to(28)Rm(ϕCm)n+1−(ϕCm)n+rΔt+[ϕdxx∂Cm∂x]wAw−[ϕdxx∂Cm∂x]eAe+[ϕdyy∂Cm∂y]nAn−[ϕdyy∂Cm∂y]sAs+[ϕdzz∂Cm∂z]tAt−[ϕdzz∂Cm∂z]bAb=∯S[Fm]N•ndS.    With the solution of Cmn+r,
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
(29a)ϕRmV(Cm)n+r1−(Cm)n+rΔt−[(Aϕdxx)e2(Cm,i+1−Cm)n+r1(Δxi+1+Δxi)−(Aϕdxx)w2(Cm−Cm,i−1)n+r1(Δxi−1+Δxi)]=[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]en+r−[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]wn+r,(29b)ϕRmV(Cm)n+r2−(Cm)n+r1Δt−[(Aϕdyy)s2(Cm,j+1−Cm)n+r2(Δyj+1+Δyj)−(Aϕdyy)n2(Cm−Cm,j−1)n+r2(Δyj−1+Δyj)]=[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]sn+r1−[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]nn+r1,(29c)ϕRmV(Cm)n+1−(Cm)n+r2Δt−[(Aϕdzz)b2(Cm,k+1−Cm)n+1(Δzk+1+Δzk)−(Aϕdzz)t2(Cm−Cm,k−1)n+1(Δzk−1+Δzk)]=[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]bn+r2−[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]tn+r2,
where r<r1<r2<1.The finite difference equations (29a), (29b), and
(29c) can be expressed as
(30a)BC1⋅Cmn+r1+BW⋅Cm,i−1n+r1+BE⋅Cm,i+1n+r1=Rs1,(30b)BC2⋅Cmn+r2+BN⋅Cm,j−1n+r2+BS⋅Cm,j+1n+r2=Rs2,(30c)BC3⋅Cmn+1+BT⋅Cm,k−1n+1+BB⋅Cm,k+1n+1=Rs3,where(31)BW=2(Aϕdxx)w(Δx+−1+Δxi),BE=2(Aϕdxx)e(Δxi+1+Δxi),BC1=ϕRmVΔt−BW−BE,Rs1=ϕRmVΔt(Cm)n+r+[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]en+r−[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]wn+r;BN=2(Aϕdyy)n(Δyj−1+Δyj),BS=2(Aϕdyy)s(Δyj+1+Δyj),BC2=ϕRmVΔt−BN−BS,Rs2=ϕRmVΔt(Cm)n+r1+[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]sn+r1−[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]nn+r1,BT=2(Aϕdzz)t(Δzk−1+Δzk),BB=2(Aϕdzz)b(Δzk+1+Δzk),BC3=ϕRmVΔt−BT−BB,Rs3=ϕRmVΔt(Cm)n+r2+[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]bn+r2−[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]tn+r2.     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.


## 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 gridijk):
(19a)∯S(CmU)•ndS=−[Cmu]wAw+[Cmu]eAe−[Cmu]nAn+[Cmu]sAs−[Cmu]tAt+[Cmu]bAb=12{(Au)e(Cm)i+1+[(Au)e−(Au)w]Cm−(Au)w(Cm)i−1}+12{(Av)s(Cm)j+1+[(Av)s−(Av)n]Cm−(Av)n(Cm)j−1}+12{(Aw)b(Cm)k+1+[(Aw)b−(Aw)t]×Cm−(Aw)t(Cm)k−1}=12((Au)e−|(Au)e|)(Cm)i+1−12((Au)w+|(Au)w|)(Cm)i−1+12((Av)s−|(Av)s|)(Cm)j+1−12((Av)n+|(Av)n|)(Cm)j−1+12((Aw)b−|(Aw)b|)(Cm)k+1−12((Aw)t+|(Aw)t|)(Cm)k−1+((Au)e−(Au)w+|(Au)e|+|(Au)w|+(Av)s−(Av)n+|(Av)s|+|(Av)n|)Cm2+((Aw)b−(Aw)t+|(Aw)b|+|(Aw)t|)Cm2+|(Au)e|(Cm)i+1−Cm2+|(Av)s|(Cm)j+1−Cm2+|(Aw)b|(Cm)k+1−Cm2−|(Au)w|Cm−(Cm)i−12−|(Av)n|Cm−(Cm)j−12−|(Aw)t|Cm−(Cm)k−12.                    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),(19b)∯S(CmU)•ndS=−[Cmu]wA^w+[Cmu]eA^e−[Cmv]nA^n+[Cmv]sA^s+[Cmw]bAb.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 vn=0,
the volume integral term is replaced by (15a) while the effective areas, A^w,A^e,A^n,A^s,
are used instead of the section areas Aw,Ae,An,As. 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 ijk) (20)AC⋅Cmn+r+AW⋅(Cm)i−1n+r+AE⋅(Cm)i+1n+r+AN⋅(Cm)j−1n+r+AS⋅(Cm)j+1n+r+AT⋅(Cm)k−1n+r+AB⋅(Cm)k+1n+r=Rs,  where(21)AW=−12((Au)w+|(Au)w|),AE=12((Au)e−|(Au)e|),AN=−12((Av)n+|(Av)n|),AS=12((Av)s−|(Av)s|),AT=−12((Aw)t+|(Aw)t|),AB=12((Aw)b−|(Aw)b|),AC=(ϕRm)ijkΔt−AWi+1−AEi−1−ANj+1−ASj−1−ATk+1−ABk−1,Rs=(ΣmV)ijk+(ϕRmCm)ijknΔt−{|(Au)e|[(Cm)i+1−Cm2]n+|(Av)s|[(Cm)j+1−Cm2]n+|(Aw)b|[(Cm)k+1−Cm2]n−|(Au)w|[Cm−(Cm)i−12]n−|(Av)n|[Cm−(Cm)j−12]n−|(Aw)t|[Cm−(Cm)k−12]n}.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 Cm can be independently derived:(22)Cmn+r=1AC[(ϕVCm)n(Δt)+Rm].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(i+1,j,k), its
solution can be calculated as(23)Cmn+ri+1,j,k=1[AC]i+1,j,k{[(ϕVkCkW)n(Δt)+Rm]i+1,j,k−[AW]i+1,j,k•Cmn+ri,j,k}.Thus we can define an index array[ID(n)=l,n=1,NX∗NY∗NZ] to record the abovementioned upwind-downwind
sequence of the grids l so that all
the values of Cmn+r(l) can be explicitly calculated one by one along
the index array ID(n). That is to say, as the solution
procedure is implemented from 1 to n+1, the value Cmn+r of grid [ID(i),i∈1,n] is already known or the grid [ID(i),i∈1,n] is NOT the upwind of l=ID(n).Now, let us discuss how to determine the index array[ID(n)=l,n=1,NX∗NY∗NZ]. Without loss of the generality, we define six
auxiliary arrays, W,E,N,S,T,B, to indicate the states of west, east,
south, north, top, and bottom boundaries for a grid ijk. For each of the six boundaries, say the west,
we assign W(ijk)=1  for influx and W(ijk)=0 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, sayijk, if the following criterion is
met: (24)W(ijk)=0,E(ijk)=0,N(ijk)=0,S(ijk=0,T(ijk)=0,B(ijk)=0, we call such grid as the starting grid and indicate it as ID(l)=ijk, 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,(25)W(ijk+1)=0,E(ijk−1)=0,N(ijk+NX)=0,S(ijk−NX)=0,T(ijk+NX∗NY)=0,B(ijk−NX∗NY)=0.Then we examine boundary states for adjacent grids
of the starting grid. If any, sayijk+1, meets the condition (24), we indicate
that ID(2)=ijk+1, as illustrated in
Figure 5. Following this procedure, we can indicate all the grids and determine
the array ID. We then calculate the
solution for all the grids indexed by the array ID in an explicit way.Figure 5
Schematic diagram of the upwind-downwind
grids.


## 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 gridijk):
(19a)∯S(CmU)•ndS=−[Cmu]wAw+[Cmu]eAe−[Cmu]nAn+[Cmu]sAs−[Cmu]tAt+[Cmu]bAb=12{(Au)e(Cm)i+1+[(Au)e−(Au)w]Cm−(Au)w(Cm)i−1}+12{(Av)s(Cm)j+1+[(Av)s−(Av)n]Cm−(Av)n(Cm)j−1}+12{(Aw)b(Cm)k+1+[(Aw)b−(Aw)t]×Cm−(Aw)t(Cm)k−1}=12((Au)e−|(Au)e|)(Cm)i+1−12((Au)w+|(Au)w|)(Cm)i−1+12((Av)s−|(Av)s|)(Cm)j+1−12((Av)n+|(Av)n|)(Cm)j−1+12((Aw)b−|(Aw)b|)(Cm)k+1−12((Aw)t+|(Aw)t|)(Cm)k−1+((Au)e−(Au)w+|(Au)e|+|(Au)w|+(Av)s−(Av)n+|(Av)s|+|(Av)n|)Cm2+((Aw)b−(Aw)t+|(Aw)b|+|(Aw)t|)Cm2+|(Au)e|(Cm)i+1−Cm2+|(Av)s|(Cm)j+1−Cm2+|(Aw)b|(Cm)k+1−Cm2−|(Au)w|Cm−(Cm)i−12−|(Av)n|Cm−(Cm)j−12−|(Aw)t|Cm−(Cm)k−12.                    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),(19b)∯S(CmU)•ndS=−[Cmu]wA^w+[Cmu]eA^e−[Cmv]nA^n+[Cmv]sA^s+[Cmw]bAb.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 vn=0,
the volume integral term is replaced by (15a) while the effective areas, A^w,A^e,A^n,A^s,
are used instead of the section areas Aw,Ae,An,As. 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 ijk) (20)AC⋅Cmn+r+AW⋅(Cm)i−1n+r+AE⋅(Cm)i+1n+r+AN⋅(Cm)j−1n+r+AS⋅(Cm)j+1n+r+AT⋅(Cm)k−1n+r+AB⋅(Cm)k+1n+r=Rs,  where(21)AW=−12((Au)w+|(Au)w|),AE=12((Au)e−|(Au)e|),AN=−12((Av)n+|(Av)n|),AS=12((Av)s−|(Av)s|),AT=−12((Aw)t+|(Aw)t|),AB=12((Aw)b−|(Aw)b|),AC=(ϕRm)ijkΔt−AWi+1−AEi−1−ANj+1−ASj−1−ATk+1−ABk−1,Rs=(ΣmV)ijk+(ϕRmCm)ijknΔt−{|(Au)e|[(Cm)i+1−Cm2]n+|(Av)s|[(Cm)j+1−Cm2]n+|(Aw)b|[(Cm)k+1−Cm2]n−|(Au)w|[Cm−(Cm)i−12]n−|(Av)n|[Cm−(Cm)j−12]n−|(Aw)t|[Cm−(Cm)k−12]n}.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 Cm can be independently derived:(22)Cmn+r=1AC[(ϕVCm)n(Δt)+Rm].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(i+1,j,k), its
solution can be calculated as(23)Cmn+ri+1,j,k=1[AC]i+1,j,k{[(ϕVkCkW)n(Δt)+Rm]i+1,j,k−[AW]i+1,j,k•Cmn+ri,j,k}.Thus we can define an index array[ID(n)=l,n=1,NX∗NY∗NZ] to record the abovementioned upwind-downwind
sequence of the grids l so that all
the values of Cmn+r(l) can be explicitly calculated one by one along
the index array ID(n). That is to say, as the solution
procedure is implemented from 1 to n+1, the value Cmn+r of grid [ID(i),i∈1,n] is already known or the grid [ID(i),i∈1,n] is NOT the upwind of l=ID(n).Now, let us discuss how to determine the index array[ID(n)=l,n=1,NX∗NY∗NZ]. Without loss of the generality, we define six
auxiliary arrays, W,E,N,S,T,B, to indicate the states of west, east,
south, north, top, and bottom boundaries for a grid ijk. For each of the six boundaries, say the west,
we assign W(ijk)=1  for influx and W(ijk)=0 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, sayijk, if the following criterion is
met: (24)W(ijk)=0,E(ijk)=0,N(ijk)=0,S(ijk=0,T(ijk)=0,B(ijk)=0, we call such grid as the starting grid and indicate it as ID(l)=ijk, 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,(25)W(ijk+1)=0,E(ijk−1)=0,N(ijk+NX)=0,S(ijk−NX)=0,T(ijk+NX∗NY)=0,B(ijk−NX∗NY)=0.Then we examine boundary states for adjacent grids
of the starting grid. If any, sayijk+1, meets the condition (24), we indicate
that ID(2)=ijk+1, as illustrated in
Figure 5. Following this procedure, we can indicate all the grids and determine
the array ID. We then calculate the
solution for all the grids indexed by the array ID 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 term(26)∯S[Fm]•n→dS=∯S[Fm]D•n→dS+∯S[Fm]N•n→dS,where(27)[Fm]D=ϕ[dxx∂Cm∂x,dyy∂Cm∂y,dzz∂Cm∂z]T,[Fm]N=ϕ[dxy∂Cm∂y+dxz∂Cm∂z,dyx∂Cm∂x+dyz∂Cm∂z,dzx∂Cm∂x+dzy∂Cm∂y]T. Applying (14a),
(14b), and
(26) to
(18b) leads to(28)Rm(ϕCm)n+1−(ϕCm)n+rΔt+[ϕdxx∂Cm∂x]wAw−[ϕdxx∂Cm∂x]eAe+[ϕdyy∂Cm∂y]nAn−[ϕdyy∂Cm∂y]sAs+[ϕdzz∂Cm∂z]tAt−[ϕdzz∂Cm∂z]bAb=∯S[Fm]N•ndS.    With the solution of Cmn+r,
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
(29a)ϕRmV(Cm)n+r1−(Cm)n+rΔt−[(Aϕdxx)e2(Cm,i+1−Cm)n+r1(Δxi+1+Δxi)−(Aϕdxx)w2(Cm−Cm,i−1)n+r1(Δxi−1+Δxi)]=[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]en+r−[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]wn+r,(29b)ϕRmV(Cm)n+r2−(Cm)n+r1Δt−[(Aϕdyy)s2(Cm,j+1−Cm)n+r2(Δyj+1+Δyj)−(Aϕdyy)n2(Cm−Cm,j−1)n+r2(Δyj−1+Δyj)]=[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]sn+r1−[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]nn+r1,(29c)ϕRmV(Cm)n+1−(Cm)n+r2Δt−[(Aϕdzz)b2(Cm,k+1−Cm)n+1(Δzk+1+Δzk)−(Aϕdzz)t2(Cm−Cm,k−1)n+1(Δzk−1+Δzk)]=[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]bn+r2−[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]tn+r2,
where r<r1<r2<1.The finite difference equations (29a), (29b), and
(29c) can be expressed as
(30a)BC1⋅Cmn+r1+BW⋅Cm,i−1n+r1+BE⋅Cm,i+1n+r1=Rs1,(30b)BC2⋅Cmn+r2+BN⋅Cm,j−1n+r2+BS⋅Cm,j+1n+r2=Rs2,(30c)BC3⋅Cmn+1+BT⋅Cm,k−1n+1+BB⋅Cm,k+1n+1=Rs3,where(31)BW=2(Aϕdxx)w(Δx+−1+Δxi),BE=2(Aϕdxx)e(Δxi+1+Δxi),BC1=ϕRmVΔt−BW−BE,Rs1=ϕRmVΔt(Cm)n+r+[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]en+r−[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]wn+r;BN=2(Aϕdyy)n(Δyj−1+Δyj),BS=2(Aϕdyy)s(Δyj+1+Δyj),BC2=ϕRmVΔt−BN−BS,Rs2=ϕRmVΔt(Cm)n+r1+[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]sn+r1−[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]nn+r1,BT=2(Aϕdzz)t(Δzk−1+Δzk),BB=2(Aϕdzz)b(Δzk+1+Δzk),BC3=ϕRmVΔt−BT−BB,Rs3=ϕRmVΔt(Cm)n+r2+[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]bn+r2−[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]tn+r2.     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.


## 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 term(26)∯S[Fm]•n→dS=∯S[Fm]D•n→dS+∯S[Fm]N•n→dS,where(27)[Fm]D=ϕ[dxx∂Cm∂x,dyy∂Cm∂y,dzz∂Cm∂z]T,[Fm]N=ϕ[dxy∂Cm∂y+dxz∂Cm∂z,dyx∂Cm∂x+dyz∂Cm∂z,dzx∂Cm∂x+dzy∂Cm∂y]T. Applying (14a),
(14b), and
(26) to
(18b) leads to(28)Rm(ϕCm)n+1−(ϕCm)n+rΔt+[ϕdxx∂Cm∂x]wAw−[ϕdxx∂Cm∂x]eAe+[ϕdyy∂Cm∂y]nAn−[ϕdyy∂Cm∂y]sAs+[ϕdzz∂Cm∂z]tAt−[ϕdzz∂Cm∂z]bAb=∯S[Fm]N•ndS.    With the solution of Cmn+r,
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
(29a)ϕRmV(Cm)n+r1−(Cm)n+rΔt−[(Aϕdxx)e2(Cm,i+1−Cm)n+r1(Δxi+1+Δxi)−(Aϕdxx)w2(Cm−Cm,i−1)n+r1(Δxi−1+Δxi)]=[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]en+r−[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]wn+r,(29b)ϕRmV(Cm)n+r2−(Cm)n+r1Δt−[(Aϕdyy)s2(Cm,j+1−Cm)n+r2(Δyj+1+Δyj)−(Aϕdyy)n2(Cm−Cm,j−1)n+r2(Δyj−1+Δyj)]=[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]sn+r1−[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]nn+r1,(29c)ϕRmV(Cm)n+1−(Cm)n+r2Δt−[(Aϕdzz)b2(Cm,k+1−Cm)n+1(Δzk+1+Δzk)−(Aϕdzz)t2(Cm−Cm,k−1)n+1(Δzk−1+Δzk)]=[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]bn+r2−[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]tn+r2,
where r<r1<r2<1.The finite difference equations (29a), (29b), and
(29c) can be expressed as
(30a)BC1⋅Cmn+r1+BW⋅Cm,i−1n+r1+BE⋅Cm,i+1n+r1=Rs1,(30b)BC2⋅Cmn+r2+BN⋅Cm,j−1n+r2+BS⋅Cm,j+1n+r2=Rs2,(30c)BC3⋅Cmn+1+BT⋅Cm,k−1n+1+BB⋅Cm,k+1n+1=Rs3,where(31)BW=2(Aϕdxx)w(Δx+−1+Δxi),BE=2(Aϕdxx)e(Δxi+1+Δxi),BC1=ϕRmVΔt−BW−BE,Rs1=ϕRmVΔt(Cm)n+r+[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]en+r−[ϕA(dxy∂Cm∂y+dxz∂Cm∂z)]wn+r;BN=2(Aϕdyy)n(Δyj−1+Δyj),BS=2(Aϕdyy)s(Δyj+1+Δyj),BC2=ϕRmVΔt−BN−BS,Rs2=ϕRmVΔt(Cm)n+r1+[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]sn+r1−[ϕA(dyx∂Cm∂x+dyz∂Cm∂z)]nn+r1,BT=2(Aϕdzz)t(Δzk−1+Δzk),BB=2(Aϕdzz)b(Δzk+1+Δzk),BC3=ϕRmVΔt−BT−BB,Rs3=ϕRmVΔt(Cm)n+r2+[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]bn+r2−[ϕA(dzz∂Cm∂x+dzy∂Cm∂y)]tn+r2.     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 system50×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. Figure7 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
Table1.Table 1
Variation of DOC with two different DO
concentrations for a well of 150 m.
Time (days)DO = 0.3 mg/LDO = 1.0 mg/L2401.2071.2053601.2401.2284801.3061.2676001.3881.3157201.4691.362


## 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.



---

*Source: 375868-2008-10-13.xml*