Convergence of iterative split-operator approaches for approximating nonlinear reactive transport problems

Advances in Water Resources 26 (2003) 247–261 www.elsevier.com/locate/advwatres

Convergence of iterative split-operator approaches for approximating nonlinear reactive transport problems Joseph F. Kanney a, Cass T. Miller

a,*

, C.T. Kelley

b

a

b

Department of Environmental Sciences and Engineering, Center for the Advanced Study of the Environment, University of North Carolina, Chapel Hill, NC 27599-7400, USA Department of Mathematics, Center for Research in Scientific Computation, North Carolina State University, Raleigh, NC 27695-8205, USA Received 6 March 2002; received in revised form 22 October 2002; accepted 15 November 2002

Abstract Numerical solutions to nonlinear reactive solute transport problems (NRTPs) are often computed using split-operator (SO) approaches, which separate the transport and reaction processes. This uncoupling introduces an additional source of numerical error, known as the splitting error. The iterative split-operator (ISO) algorithm removes the splitting error through iteration. Although the ISO algorithm is often used, there has been very little analysis of its convergence behavior. This work uses theoretical analysis and numerical experiments to investigate the convergence rate of the ISO approach for solving NRTPs. We show that under certain assumptions regarding smoothness, the convergence rate of the ISO algorithm applied NRTPs is OðDt2 Þ. We demonstrate that the theoretical convergence rate can be achieved in practice if the numerical solution of the transport and reaction steps are carried out with sufficient accuracy. We also show that accurate estimation of the lagged operator in each step is crucial to obtaining the theoretical convergence rate. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Reactive transport; Nonlinear mass transfer; Numerical solution; Operator splitting; Iterative split-operator; Convergence rate

1. Introduction Nonlinear reactive transport problems (NRTPs), i.e. problems involving advective and diffusive/dispersive solute transport combined with nonlinear homogeneous and/or heterogeneous reactions, are important in many application areas. Examples include transport of contaminants, dissolved minerals, and microorganisms in groundwater systems [1,13,15,16,28,43,47,56,63,74]; contaminant and nutrient transport in oceans, lakes and rivers [19,20,34,44,60,78]; pollution transport in the atmosphere [9,29,36,66,75]; as well as a variety of industrial applications, e.g. [3,10,12,37,38,49,55]. The governing system of equations for NRTPs typically consist of partial differential equations (PDEs) for transport coupled with algebraic equations (AEs) and/or ordinary differential equations (ODEs) describing chemical equilibrium and chemical kinetics, respectively. These systems seldom admit analytical solutions, so *

Corresponding author. Fax: +1-919-966-2583. E-mail addresses: [email protected] (J.F. Kanney), casey_ [email protected] (C.T. Miller), [email protected] (C.T. Kelley).

numerical solution techniques are commonly used. Problems that are advective dominated or involve stiff chemical reactions are difficult to solve, even numerically. For large-scale problems and/or problems involving multiple species, the computational expense of obtaining accurate solutions can be substantial [80]. There are two fundamentally different algorithmic approaches to the numerical solution of NRTPs: (1) the fully coupled (FC) approach in which the discrete forms of the governing equations are solved as a single system; and (2) the split-operator (SO) approach, which uses a time-splitting method to uncouple and separately solve the discretized governing equations for transport and reaction. The SO algorithm is very attractive for large, multidimensional problems with multiple species because the memory requirements of the individual substeps of the uncoupled problem are significantly smaller than the FC approach [5,68,80]. The SO approach allows one to combine a transport-only computer code with a chemistry-only computer code with minimal effort [8,33]. The SO approach also lends itself to parallel computation, since the reaction step may be computed independently at each node in the computational grid

0309-1708/03/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 3 0 9 - 1 7 0 8 ( 0 2 ) 0 0 1 6 2 - 8

248

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

Nomenclature Roman Letters C volume averaged aqueous phase solute concentration [ML3 ] Dx dispersion coefficient in x direction [L2 T1 ] e difference of two general functions [1] E error in a numerical solution [1] f forcing function in simple ODE [–] F general operator acting on Rn -valued functions [–] ka aqueous phase reaction rate coefficient [Mð1ra Þ L3ðra 1Þ T1 ] kf rapidly sorbing solid phase fraction reaction rate coefficient [T1 ] ks slowly sorbing solid phase fraction reaction rate coefficient [T1 ] K general constant in algorithm analysis [–] Kf rapidly sorbing solid phase fraction Freundlich coefficient [L3nf Mnf ] Ks slowly sorbing solid phase fraction Freundlich coefficient [L3ns Mns ] Lt transport operator [–] Lr1 reaction operator [–] Lr2 reaction operator [–] m number of chemical species [–] M general constant in algorithm analysis [–] ne number of equations [–] nf rapidly sorbing solid phase fraction Freundlich exponent [–] ns slowly sorbing solid phase fraction Freundlich exponent [–] N number of nodes in spatial discretization [–] O asymptotic order [–] Q rate of convergence of ISO iterations [–] ra aqueous phase reaction order [–] rf rapidly sorbing solid phase fraction reaction order [–] rs slowly sorbing solid phase fraction reaction order [–] R general operator acting on Rn -valued functions [–] R0 bound for general operator acting on Rn valued functions [–] Rf retardation factor [–] Rn n-dimensional Euclidean space [Ln ] t temporal coordinate [T] ts simulation length [T] Dt time increment [T] u general Rn -valued function [–] v dependent variable in simple ODE [–] vx pore velocity component in x direction in PDE [LT1 ] w dependent variable in simple ODE [–]

x xl Dx

spatial coordinate [L] one-dimensional spatial domain length [L] space increment [L]

Greek Letters a first order mass transfer coefficient for mass transfer between aqueous phase and slowly sorbing solid phase fraction [T1 ] k
kk discrete vector norm of index k, where k ¼ 1; 2; 1 [–] x mass averaged solid phase solute concentration [MM1 ] c Lipschitz constant [–] q solid phase particle density [ML3 ] U denotes the solution ðC; xs Þ of the system of transport equations [–] qb solid phase bulk density [ML3 ] s temporal discretization error of numerical solution [–] h solid phase porosity (void space volume fraction) [–] Subscripts a qualifier denoting quantity in aqueous phase f qualifier denoting quantity in rapidly sorbing solid phase fraction i iteration index in ISO algorithm r1 operator qualifier denoting reaction r2 operator qualifier denoting reaction s qualifier denoting quantity in slowly sorbing solid phase fraction t operator qualifier denoting transport 0 qualifier denoting value of quantity at x ¼ 0, the inlet of one-dimensional domain Superscripts  qualifier denoting an exact solution 0 qualifier denoting quantity evaluated at t ¼ 0 1=2 qualifier denoting the quantity evaluated after the transport solve in ISO algorithm þ qualifier denoting the quantity evaluated after the reaction solve in ISO algorithm c qualifier denoting the quantity evaluated using the current approximate solution e qualifier denoting quantity measured at thermodynamic equilibrium qualifier denoting estimate obtained using intermediate values of dependent variables in split-operator algorithms Abbreviations 1D one-dimensional 2D two-dimensional 3D three-dimensional

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

AE algebraic equation ADRE advection–dispersion–reaction equation ASO alternating split-operator numerical solution algorithm DSA direct substitution numerical solution algorithm FC fully coupled numerical solution algorithm GIA global implicit numerical solution algorithm ISO iterative split-operator numerical solution algorithm

[80]. Due to these advantages, SO approaches have been widely used to solve NRTPs [39,68]. A significant drawback of the SO approach is that decoupling the governing equations introduces an additional source of numerical error, referred to as splitting error [73]. The splitting error has been shown to depend upon the order of the solution of the transport and reaction subproblems [5,73]. Furthermore, the splitting error may be removed by iteration [27], which we refer to as the iterative split-operator (ISO) algorithm. The applicability of ISO depends mainly on the stability requirements and the rate of convergence of the iterations [27,80]; others have reported difficulties in getting the iterative procedure to converge [24,30,68]. Despite this difficulty, and despite the wide application of SO approaches, there are relatively few studies of the convergence, computational accuracy, and efficiency of the method for general nonlinear problems. There are even fewer studies of the stability and convergence properties of the ISO algorithm. The overall goal of the work reported here is to investigate the convergence of the ISO algorithm for NRTPs. The specific objectives we pursue are (1) to perform a theoretical analysis of the convergence rate for the ISO algorithm solution of a general NRTP form; and (2) to investigate the implications of the analytical results through numerical experiments using a simplified ODE model and a model NRTP with features relevant to a wide class of applications in water resources.

2. Background 2.1. Algorithms A wide variety of solution methods have been formulated and applied to solve NRTPs, and summarizing them all is beyond the scope of this paper. More detailed discussions of solution methods can be found in [4,32,46,52,68]. As mentioned previously, solution methods may be classified algorithmically as either FC or SO approaches. The FC approach is sometimes called

MOL NRTP ODE PDE SIA SO SNIA SSO

249

method of lines numerical solution approach nonlinear reactive transport problem ordinary differential equation partial differential equation sequential iterative approach split-operator numerical solution approach sequential noniterative approach sequential split-operator numerical solution algorithm

the direct substitution approach (DSA) [61,80], the global implicit approach (GIA) [52,68], or the one-step approach [27,68]. The fundamental idea of the FC approach is to solve the PDEs describing transport simultaneously with the AEs or ODEs describing the equilibrium or kinetic reaction system, respectively. Although mathematically rigorous, practical implementation of the FC approach is limited by available computer memory, since the size of the coefficient matrix in the discretized system grows as a product of the number of spatial discretization intervals and the number of species. In general, the FC approach involves solving a N m system of nonlinear AEs at each time step, where N is the number nodes in the spatial mesh and m is the number of chemical components. SO approaches are commonly used in a wide variety of disciplines where transport of reacting species is important. Notable examples include the study of transport in subsurface systems [2,17,24,25,28,42,48,53,56, 57,62,65,69,74,77,79,81], in surface water systems [35,51,70], and in the atmosphere [11,33,46,67,72]. The basic concept of the SO approach is to solve the transport and reaction equations separately. The decoupled problem consists of m transport equations, each with N unknowns, and a system of m reaction equations at each of the N nodes. Thus, each of the subproblems has a smaller memory requirement than the FC approach. The SO approach provides a framework that facilitates combining existing transport solvers with different reaction modules. The SO approach also lends itself to parallel solution, since the reaction subproblem at each node can be solved independently of the other nodes in the system. Application of SO approaches in the subsurface environment has been reviewed elsewhere [48,80]. Several variations of the SO approach are possible depending upon the order in which the subproblems are solved and whether or not iteration is involved. The magnitude of the error introduced by splitting depends in part on the ordering and the iteration procedure. Variants of the SO approach include the sequential SO (SSO) algorithm, the alternating SO (ASO) algorithm, and the ISO algorithm.

250

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

SSO, also known as the sequential noniterative approach (SNIA) [80], is typically formulated by (1) solving the transport portion of the problem over a time interval––neglecting reactions and mass transfer, and (2) solving for the reactions and mass transfer over a time interval of the same length. Each portion of the solution over the time interval used may include one or more substeps, e.g., the reaction step may be solved using multiple steps of an ODE solver [48]. ASO, which is closely related to Strang splitting [71] for hyperbolic PDEs, is typically formulated by (1) solving the transport portion of the problem over onehalf of the time interval, (2) solving for mass transfer and reaction terms over the entire time interval, and (3) solving the transport portion of the problem over the second half of the time interval. After an initial transport solve over a one-half interval of time, this algorithm reduces to time marching by solving for the transport and reaction steps in an alternating fashion over an entire time interval, with a half-step transport solve at time points where a solution is needed. As with the SSO algorithm, each portion of the solution may include one or more substeps. ISO, also known as the sequential iterative approach (SIA) [27,68,80], is typically formulated by (1) solving the transport portion of the problem over a full time interval, assuming the reaction contributions are known; (2) solving the reaction portion of the problem over a full time interval, assuming the transport contributions are known; and (3) iterating over the first two steps in the algorithm until a convergence criterion is satisfied. 2.2. Error analysis In the water resources literature, several researchers have investigated the convergence of the SO approach applied to advection–dispersion–reaction equations (ADREs). Although work on hyperbolic problems, has shown that time splitting can introduce numerical error for certain classes of operators [40,71], convergence of the SO formulation is expected as the numerical time step approaches zero. Wheeler, Dawson and co-workers [21,22,76] have constructed formal proofs for the convergence of the SO approach applied to nonlinear ADREs arising in contaminant transport. Herzer and Kinzelbach [27] used Taylor series analysis and numerical experiments applied to a linear ADRE to show that the error introduced by time splitting in the SSO algorithm can be viewed as an additional source of numerical dispersion. Valocchi and Malmstead [73] used a heuristic analysis for a linear ADRE to show that the magnitude of splitting error for SSO is OðDtÞ, where Dt is the time step in the numerical method, which agreed with results in Herzer and Kinzelbach. Numerical experiments indicated that the splitting error for the ASO algorithm could be an order

of magnitude less than that of SSO for this problem. Barry et al. showed that the splitting error of the ASO algorithm is OðDt2 Þ for ADREs with linear reactions [7], but that this result does not hold for nonlinear reactions [5]. Subsequent analysis of SSO and ASO applied to linear and nonlinear ADREs [31,50] confirmed earlier results [5,7,6,73] and highlighted the dependence of the splitting error on the magnitude of the reaction rates. 2.3. Comparison of ISO with other methods The practical value of applying ISO to NRTPs depends upon the rate at which the splitting error is reduced. If the rate is rapid enough, ISO will be computationally efficient compared to alternative methods. Thus, formal analysis of the convergence rate of ISO applied to NRTPs and/or numerical comparisons of ISO to alternative solution algorithms are useful. We are not aware of any formal analysis on the convergence rate of the ISO algorithm applied to NRTPs, but numerical comparisons of the ISO to other SO approaches and to FC approaches have appeared in the literature over the last decade [27,61,68,80,83]. Herzer and Kinzelbach [27] applied ISO to a linear ADRE and showed it to be a computationally efficient approach for reducing the error when compared to simply reducing the time step in the SSO algorithm. Yeh and Tripathi [80] analyzed computer memory requirements and estimated computational time requirements, for both the ISO and the FC approach, for twodimensional (2D) and three-dimensional (3D) NRTPs. They concluded that only ISO is viable for realistic applications. The FC approach was judged to require excessive memory and computational time. Later work on biodegradation in groundwater systems [82,83] advocated the use of ASO for NRTPs with kinetic chemistry and ISO for NRTPs with equilibrium chemistry. The most extensive comparison that we know of is that of Steefel and MacQuarrie [68], which considered example problems ranging from simple decay and equilibrium sorption reactions to more highly nonlinear Monod kinetics and multicomponent sorption problems. In some cases ISO exhibited the greatest efficiency, while in other cases SSO was more efficient. The FC approach was the least efficient from a computational viewpoint. But recent work comparing the ISO approach to the FC approach for several NRTP test cases reports that ISO is more efficient than the FC approach only for large, chemically simple problems [61]. These results seem to contradict, to some extent, the results of other researchers. The results reported above indicate clear evidence that, at least in some cases, the convergence of the ISO algorithm is sufficiently rapid to make it computationally efficient compared to other SO approaches and the FC approach. But the results of these investigations

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

251

indicate that the optimal algorithm may be problem dependent, and formal analysis of the convergence behavior of the ISO algorithm has yet to appear in the literature.

Having computed U1=2 , then Uþ is computed as the solution of

3. Analysis

oxþ þ s ¼ Lr2 ðC þ ; xþ ð7Þ s Þ ¼ Lr2 ðU Þ ot with initial condition Uþ ðt0 Þ ¼ ðC 0 ; x0s Þ. Here e L t is an approximation. In this work, we use e ð8Þ L t ¼ Lt ðC 1=2 Þ:

An overview of the mathematical formulation of reactive transport problems in subsurface applications is given by Rubin [59] and will not be repeated here. But in order to describe the ISO algorithm and to motivate the following analysis, we consider a multicomponent reactive transport problem written in a general form using operator notation Rf ðCÞ

oC ¼ Lt ðCÞ þ Lr1 ðC; xs Þ; ot

ð1Þ

Rf ðC þ Þ

oC þ e þ e ¼ L t þ Lr1 ðC þ ; xþ s Þ ¼ L t þ Lr1 ðU Þ; ot

ð6Þ

As long as the spatial discretization is fixed over the course of the ISO iterations, we can ignore the spatial variables because they have no role in the analysis and the results will hold in one-dimensional (1D), 2D and 3D domains. 3.2. Preliminary estimate

oxs ¼ Lr2 ðC; xs Þ; ot

ð2Þ

where Rf ðCÞ, Lt ðCÞ; Lr1 ðC; xs Þ; Lr2 ðC; xs Þ are operators. C and xs are species concentrations and mass fractions, respectively. In a typical application, Lt ðCÞ is a linear operator containing spatial derivatives associated with advective and dispersive processes, while Rf ðCÞ; Lr1 ðC; xs Þ; Lr2 ðC; xs Þ are nonlinear operators associated with reaction and mass transfer processes.

We begin by defining the norms used throughout the analysis, then introducing a lemma which will be used twice. The norms and the lemma will be expressed in terms of a general function, u, and general operators R; F1 , and F2 . Approximations to R and F2 are denoted e and Fe2 , respectively. by R If u : ½0; Dt ! Rn is continuous, then we define kuk1 ¼ max kuðtÞk;

ð9Þ

0 6 t 6 Dt

3.1. ISO solution algorithm

where k  k is the Euclidean norm on Rn .

The ISO algorithm can be expressed as a two-step procedure in which the first step solves an approximation to Eq. (1), which is typically a transport equation, and the second step solves an approximation of a local system formed by Eqs. (1) and (2), typically reaction equations at each location in space. The initial conditions for each step of the ISO algorithm is the solution from the previous time step, or the initial conditions for the first time step, although other choices are possible and may lead to a more efficient algorithm. Specifically, we consider the time interval t 2 ½t0 ; t0 þ Dt with initial conditions Cðt0 Þ ¼ C 0 and xs ðt0 Þ ¼ x0s . We adopt the notation U ¼ ðC; xs Þ for the solution to the system of equations. Any stage of the iteration will move from a current approximate solution Uc to a new one Uþ through an intermediate step U1=2 . In the first step, U1=2 ¼ ðC 1=2 ; xs1=2 Þ is computed as the solution of

Lemma 3.1. Assume that u and u are twice Lipschitz continuously differentiable Rn -valued functions in [0; Dt] and

1=2 e f oC ¼ Lt ðC 1=2 Þ þ e R L r1 ; ot

xs1=2

¼

C 1=2 ðt0 Þ ¼ C 0 ;

xcs ;

ð3Þ

0 < t < Dt; u ð0Þ ¼ u0 ð10Þ

and e ut ¼ F1 ðuÞ þ Fe2 ; R e 1

0 < t < Dt; uð0Þ ¼ u0 : 1

ð11Þ

k 6 R1 0

for some R0 > 0. Then if Assume that k R k; kR R; F1 , and F2 are Lipschitz continuous and Dt is sufficiently small, then there is a K > 0 such that   ek : ku  u k1 6 KDt kF2 ðu Þ  Fe2 k1 þ kRðu Þ  R 1 ð12Þ Proof. Set e ¼ u  u , so that et ¼ ut  ut . Subtract Eq. (10) from Eq. (11) to obtain e 1 ðF1 ðuÞ þ Fe2 Þ  Rðu Þ1 ½F1 ðu Þ þ F2 ðu Þ : et ¼ R

ð4Þ

e f and e L r1 are simple approximations. In this where R analysis, we use e f ¼ Rf ðC c Þ and e R L r1 ¼ Lr1 ðUc Þ:

Rðu Þut ¼ F1 ðu Þ þ F2 ðu Þ;

ð5Þ

e 1



ð13Þ



Adding and subtracting R ½F1 ðu Þ þ F2 ðu Þ and rearranging yields e 1 ½F1 ðuÞ  F1 ðu Þ þ ð R e 1  Rðu Þ1 ÞF1 ðu Þ et ¼ R e 1 ½ Fe2  F2 ðu Þ þ ð R e 1  Rðu Þ1 ÞF2 ðu Þ: þR

ð14Þ

252

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

Combining the second and fourth terms in Eq. (14) yields  e e 1 ½F1 ðuÞ  F1 ðu Þ þ Rðu Þ  R et ¼ R e Rðu Þ R e 1 ½ Fe2  F2 ðu Þ :

½F1 ðu Þ þ F2 ðu Þ þ R ð15Þ Let c1 be the Lipschitz constant of F1 . By assumption  1  R e ½F1 ðuÞ  F1 ðu Þ  6 R1 c1 kek1 ; ð16Þ 0 1 and therefore, letting      e M ¼ R2 0 kF1 ðu Þk1 þ kF2 ðu Þk1 kRðu Þ  R k1 þ

 R1 0 kF2 ðu Þ

 Fe2 k1 ;

ð17Þ

we have ket ðtÞk 6 R1 0 c1 kek1 þ M:

ð18Þ

Hence, for all 0 < t < Dt, Z t  Z t    et ðsÞ ds ½R1 keðtÞk ¼  0 c1 kek1 þ M ds: 6 0

ð19Þ

0

ð24Þ

In order to apply Lemma 3.1 to Eqs. (6) and (7) we must incorporate the truncation error into the Fe2 term. Theorem 3.3. Let Uþ be the solution to Eqs. (6) and (7), and let sþ be the truncation error in the integration. Then there is a K þ such that kEUþ k1 6 K þ DtðDtkEUc k þ s1=2 Þ þ sþ :

ð25Þ

Proof. We will apply Lemma 3.1 to the complete system. Here  Rf 0 ; ð26Þ R¼ 0 I 

Lr1 ðUÞ ; L ðUÞ  r1 1=2 Lt C Fe2 ¼ : 0

F1 ¼

So,

 F2 ¼

Lt C þ ; 0

and ð27Þ

The error in F2 does not depend on Ex , in fact

kek1 6 DtðR1 0 c1 kek1 þ MÞ: Assuming that 1  K¼

1=2

kEC k1 6 K 1=2 DtkEUc k1 þ s1=2 :

DtR1 0 c1

 maxðR2 0 ðkF1 ðu Þk1

completes the proof.

ð20Þ < 1=2 and setting

ð28Þ

þ kF2 ðu Þk1 Þ; R1 0 Þ

ð21Þ



3.3. Applications of the estimate We now apply the estimate given in Lemma 3.1 to the two steps of the ISO algorithm described in Eqs. (3) and (4), (6) and (7), respectively. In each application, the general functions and operators of the lemma are replaced by specific functions and operators of the ISO step in question. We use Lemma 3.1 to estimate kUþ  U k1 . We begin by showing that the error in C 1=2 is, up to the truncation error of the solver for Eq. (3), OðDtÞ smaller than the error in Uc . We let EC ¼ C  C  ;

1=2 kF2 ðU Þ  Fe2 k1 6 MkEC k1 6 MðK 1=2 DtkECc k1 þ s1=2 Þ:

Ex ¼ xs  xs ;

Eq. (25) follows directly from the lemma and Eq. (28).  As was the case with s1=2 , sþ reflects the truncation error of the solution. To summarize, Theorem 3.3 says that, within the truncation error of the solution methods, the error in Uþ should be OðDt2 Þ smaller than the error in Uc . 4. Numerical experiments To confirm the results of the theoretical analysis, we conducted numerical investigations on a series of ODE systems and on a model PDE problem that reflects a realistic model of reactive transport in the subsurface environment.

EU ¼ maxðEC ; Ex Þ ð22Þ

Lemma 3.2. Assume that Eq. (5) holds. Let C 1=2 be the solution to Eq. (3); then there is K 1=2 > 0 such that 1=2

kEC k1 6 K 1=2 DtkEUc k1 :

ð23Þ

ef ¼ R e; Proof. We apply Lemma 3.1 to Eq. (3) with R Lt C ¼ F1 ; Lr1 ðC  ; xs Þ ¼ F2 ðC  Þ, and e L r1 ¼ Fe2 .  1=2

Lemma 3.2 does not imply that the solution C that is computed in practice satisfies Eq. (23). If s1=2 is the truncation error in the solution of Eq. (3) then

4.1. A simple ODE We first looked at a simple ODE system 0

v ðtÞ ¼ f1 ðtÞ  5vðtÞ þ vðtÞwðtÞ;

ð29Þ

w0 ðtÞ ¼ f2 ðtÞ  vðtÞw2 ðtÞ;

ð30Þ

vð0Þ ¼ 0;

ð31Þ

wð0Þ ¼ 1;

ð32Þ

in which the temporal derivative term is linear in v. The forcing functions fj are constructed to make the solution of Eqs. (29) and (30)

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

v ðtÞ ¼ sinðtÞ;

ð33Þ

w ðtÞ ¼ expðtÞ:

ð34Þ

So, in the language of Section 3.1, we have Rf ðvÞ ¼ 1;

ð35Þ

Lt ðvÞ ¼ 5v;

ð36Þ

Lr1 ðv; wÞ ¼ f1 ðtÞ þ vw;

ð37Þ

Lr2 ðv; wÞ ¼ f2 ðtÞ  vðtÞw2 ðtÞ:

ð38Þ

For U ¼ ðv; wÞ, we can express the error in the ith iterate as Ei ¼ EðUi Þ ¼ maxðjvi ðDtÞ  v ðDtÞj; jwi ðDtÞ  w ðDtÞjÞ: ð39Þ From the results in Section 3.3, one would expect that, for sufficiently small Dt and sufficiently small truncation error, the convergence rate Qi ðDtÞ ¼ Ei =Ei1 ¼ OðDt2 Þ

ð40Þ

and hence that Qi ðDtÞ=Qi ð2DtÞ  0:25:

ð41Þ

We can verify this prediction through numerical experiments since we can compute the error Ei and convergence rate Qi as a function of the iteration index for several values of Dt, and look at the ratios Qi ðDtÞ=Qi ð2DtÞ. The initial value problem given by Eqs. (29) and (30) was solved in MATLAB using the ODE15s code [64] with relative and absolute error tolerances set at 1013 .

253

We directed the code to output the solution at t ¼ Dt and used Hermite interpolation to obtain solution estimates for t 2 ½0; Dt used in computing estimates of the lagged operators in each step. The initial iterate was v0 ðtÞ ¼ 0;

w0 ðtÞ ¼ 1;

for t 2 ð0; DtÞ:

ð42Þ

In Table 1 we show the error in the initial iterate and three subsequent iterates for six values of Dt (Dtj ¼ 2j ; j ¼ 2; 3; . . . ; 7). In Table 2 we show the convergence rates Qi ; i ¼ 1; 2; 3. In Table 3 we see the ratio of convergence rates predicted by the theory and Eq. (41). The interpolation error for Hermite interpolation is OðDt4 Þ, which only begins to affect the convergence rates on the 3rd iteration. 4.2. A simple ODE with errors Here we repeat the computations from Section 4.1 under conditions more likely to be present in actual computations, i.e., the integration of the split equations and the operator estimates contain error. We introduce several types of error by: (1) replacing Hermite interpolation with linear interpolation; (2) increasing the tolerance of the ODE solver in both split equations; and (3) reducing the order and the accuracy of only the transport solve. The effect of any of these errors will be to increase the truncation error effects in Eqs. (24) and (25) and make Eq. (41) fail when the truncation error terms dominate the right sides of Eqs. (24) and (25). 4.2.1. Interpolation error We repeat the computations in Section 4.1, except now we use linear interpolation to estimate e L r1 and e Lt

Table 1 ODE with Hermite interpolation: Errors Ei i n Dt

2.5000e01

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

0 1 2 3

2.4740e01 1.8552e02 1.2395e03 8.4003e05

1.2467e01 2.3828e03 4.3203e05 7.9894e07

6.2459e02 3.0155e04 1.4199e06 6.9112e09

3.1245e02 3.7921e05 4.5465e08 5.7793e11

1.5624e02 4.7542e06 1.4381e09 3.1157e13

7.8124e03 5.9516e07 4.5367e11 1.6342e13

Table 2 ODE with Hermite interpolation: Ratios Qi ¼ Ei =Ei1 i n Dt

2.5000e01

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

1 2 3

7.4985e02 6.6811e02 6.7774e02

1.9112e02 1.8131e02 1.8493e02

4.8279e03 4.7086e03 4.8675e03

1.2137e03 1.1990e03 1.2711e03

3.0428e04 3.0250e04 2.1665e04

7.6181e05 7.6226e05 3.6022e03

Table 3 ODE with Hermite interpolation: Estimated rates Qi ðDtÞ=Qi ð2DtÞ i n Dt

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

1 2 3

2.5488e01 2.7138e01 2.7286e01

2.5261e01 2.5970e01 2.6321e01

2.5138e01 2.5463e01 2.6115e01

2.5072e01 2.5230e01 1.7044e01

2.5036e01 2.5199e01 1.6627e+01

254

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

Table 4 ODE with linear interpolation: Errors Ei i n Dt

2.5000e01

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

0 1 2 3

2.4740e01 1.2823e02 2.7850e04 8.7598e04

1.2467e01 1.9587e03 2.3675e05 5.3147e05

6.2459e02 2.7272e04 2.0909e06 3.2574e06

3.1245e02 3.6042e05 1.6016e07 2.0131e07

1.5624e02 4.6343e06 1.1136e08 1.2503e08

7.8124e03 5.8759e07 7.3461e10 7.7869e10

Table 5 ODE with linear interpolation: Estimated rates Qi ðDtÞ=Qi ð2DtÞ i n Dt

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

1 2 3

3.0312e01 5.5652e01 7.1370e01

2.7793e01 6.3427e01 6.9401e01

2.6418e01 5.7962e01 8.0680e01

2.5713e01 5.4074e01 8.9329e01

2.5357e01 5.2029e01 9.4407e01

over the transport and reaction solves, respectively (Tables 4 and 5). The error in linear interpolation is OðDt2 Þ, the same order as the reduction in errors. As expected, we see interpolation effects in the second iteration, shown in the second row of Table 5. 4.2.2. Truncation error We also investigate the effect on the ISO convergence rate of increased truncation error in the ODE solver (Tables 6 and 7). We consider the case when both steps are solved to the same tolerance, and also the case where one step is solved less accurately than the other. First we repeat the computations of Section 4.1 using relative and absolute error tolerances set at 104 in the ODE solver ODE15s. The first row of Table 7 indicates that as ðDtÞ is reduced, the truncation error will begin to dominate, and second-order convergence is not evident.

The ratio Qi ðDtÞ=Qi ð2DtÞ actually increases as ðDtÞ is decreased. FC solutions to NRTPs are typically implemented using low-order finite difference or finite element methods which are second-order accurate in space and time. Although typical ISO implementations solve the reaction step using methods which are higher order in time, the transport step is often solved using methods that are lower order in time and often use fixed time steps. To model this scenario, we repeat the calculations of Section 4.1 but solve the first step of the ISO algorithm (Eq. (3)) using a second-order accurate integration method and adjust the solver tolerance such that the solution over the timestep Dt is completed in only one or two substeps (Tables 8 and 9). The second step of the ISO algorithm (Eqs. (6) and (7)) is computed using the same adaptive time stepping and tolerance used previously in

Table 6 ODE with solver tolerance ¼ 104 : Errors Ei i n Dt

2.5000e01

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

0 1 2 3

2.4740e01 1.8148e02 1.1167e03 1.3209e04

1.2467e01 2.3925e03 1.8007e04 1.7829e04

6.2459e02 3.3655e04 1.5880e04 1.5875e04

3.1245e02 1.0494e04 1.0582e04 1.0582e04

1.5624e02 4.8641e05 4.8756e05 4.8756e05

7.8124e03 2.0856e05 2.0866e05 2.0866e05

Table 7 ODE with solver tolerance ¼ 104 : Estimated rates Qi ðDtÞ=Qi ð2DtÞ i n Dt

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

1 2 3

2.6161e01 1.2232eþ00 8.3700eþ00

2.8079e01 6.2691eþ00 1.0097eþ00

6.2335e01 2.1369eþ00 1.0003eþ00

9.2687e01 9.9410e01 1.0000eþ00

8.5751e01 9.9814e01 1.0000eþ00

Table 8 ODE with one low-order substep: Ei i n Dt

2.5000e01

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

0 1 2 3

2.4740e01 1.8446e02 3.7074e04 1.4908e03

1.2467e01 2.8784e03 5.3095e05 1.1341e04

6.2459e02 4.2037e04 5.4449e06 8.1405e06

3.1245e02 5.7744e05 4.4909e07 5.5323e07

1.5624e02 7.6062e06 3.2553e08 3.6211e08

7.8124e03 9.7747e07 2.1970e09 2.3185e09

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

255

Table 9 ODE with one low-order substep: Estimated rates Qi ðDtÞ=Qi ð2DtÞ i n Dt

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

1 2 3

3.0966e01 9.1776e01 5.3118e01

2.9151e01 7.0221e01 6.9994e01

2.7460e01 6.0044e01 8.2396e01

2.6341e01 5.5029e01 9.0297e01

2.5701e01 5.2517e01 9.4873e01

Section 4.1. The rates in Table 9 show that we deviate from the second-order convergence, but the deviation appears smaller as ðDtÞ is decreased. This behavior was also observed in Table 5, where we reduced the interpolation accuracy. 4.3. A more complicated ODE Many subsurface transport modeling efforts employ equilibrium mass transfer relationships that lead to an Rf that is not Lipschitz continuous near C ¼ 0 [54]. The theoretical analysis assumes that Rf is Lipschitz, but one might ask if this is a necessary condition. To investigate the sensitivity of the analysis to the Lipschitz assumppffiffiffi tion, we consider ODE systems with Rf ðvÞ ¼ 1= v and Rf ðvÞ ¼ 1=ð1 þ v2 Þ. In the first case, the Rf has a form

that arises in subsurface modeling applications and for which the assumption of Lipschitz continuity is violated. The second form is sufficiently smooth. One might expect that the first form degrades the ISO convergence rate the second does not. The results for Rf ðvÞ ¼ pffiffiwhile ffi 1= v in Tables 10 and 11 and Rf ðvÞ ¼ 1=ð1 þ v2 Þ in Tables 12 and 13 confirm this expectation. 4.4. A nonlinear reactive transport problem The previous numerical experiments using ODEs provide some insight into the convergence behavior of the ISO algorithm but do not guarantee that the same behavior will be observed in more realistic situations. Thus, we also consider the numerical solution of a model NRTP of the form described by Eqs. (1) and (2).

Table 10 pffiffiffi ODE with Rf ðvÞ ¼ 1= v: Errors Ei i n Dt

2.5000e01

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

0 1 2 3

2.4740e01 1.4964e02 1.9867e04 5.1079e06

1.2467e01 3.6734e03 2.7798e05 5.6116e06

6.2459e02 8.1460e04 3.8206e06 6.1460e06

3.1245e02 1.6783e04 5.0851e06 5.0852e06

1.5624e02 3.2646e05 6.3568e07 6.3569e07

7.8124e03 5.7186e06 7.9464e08 7.9465e08

Table 11 pffiffiffi ODE with Rf ðvÞ ¼ 1= v: Estimated rates Qi ðDtÞ=Qi ð2DtÞ i n Dt

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

1 2 3

4.8714e01 5.6996e01 7.8519eþ00

4.4265e01 6.1980e01 7.9684eþ00

4.1186e01 6.4600eþ00 6.2166e01

3.8898e01 6.4267e01 9.9999e01

3.5033e01 7.1363e01 9.9999e01

Table 12 ODE with Rf ðvÞ ¼ 1=ð1 þ v2 Þ: Errors Ei i n Dt

2.5000e01

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

0 1 2 3

2.4740e01 2.4247e02 2.7184e03 2.9909e04

1.2467e01 2.6466e03 5.9384e05 1.3457e06

6.2459e02 3.1588e04 1.6382e06 8.7217e09

3.1245e02 3.8757e05 4.8655e08 6.4374e11

1.5624e02 4.8048e06 1.4864e09 3.3619e13

7.8124e03 5.9827e07 4.6109e11 1.6334e13

Table 13 ODE with Rf ðvÞ ¼ 1=ð1 þ v2 Þ: Estimated rates Qi ðDtÞ=Qi ð2DtÞ i n Dt

1.2500e01

6.2500e02

3.1250e02

1.5625e02

7.8125e03

1 2 3

2.1660e01 2.0014e01 2.0596e01

2.3824e01 2.3113e01 2.3493e01

2.4528e01 2.4206e01 2.4851e01

2.4791e01 2.4642e01 1.7095e01

2.4902e01 2.4914e01 1.5662eþ01

256

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

4.4.1. Model formulation For convenience, we consider a model problem introduced in [32], which describes solute transport, reaction, and interphase mass transfer in a porous medium system. Mass transfer occurs between a fluid phase and two types of solid phase––a fraction that achieves equilibrium quickly and a fraction that achieves equilibrium slowly. Transformation reactions are considered to be of a general form and may be linear or nonlinear. Subsets of this general model occur routinely in application, and because of this it makes a good test problem for the solution issues of concern in this work. Because the major aspects of this work do not depend upon the spatial dimensionality of the system, we will restrict our efforts to one spatial dimension. We also limit ourselves to the case of transport with a known uniform flow field. This subset of the governing equations describing 1D transport in a uniform flow field is often employed for problems such as solute transport through soil columns, e.g. [18,23,41], fixed bed reactors, e.g. [3,37,49], and in many field-scale investigations, e.g. [1,58,63]. We formulate this model as oC qb oxf o2 C oC q þ  ka C ra  b ½aðxes  xs Þ ¼ Dx 2  vx ot ox ox h ot h ð43Þ þ kf xrf f ; in ½0; Lx ½0; ts ; oxs ¼ aðxes  xs Þ  ks xrss ; in ½0; Lx ½0; ts ; ot

centration, i.e., mass fraction; xf and xs are solid-phase solute mass fractions for the rapidly sorbing and slowly sorbing solid phase fractions, respectively; and xe ; xef and xes are equilibrium solute mass fractions for the total, the rapidly sorbing and slowly sorbing solid phase fraction, respectively, in equilibrium with the aqueousphase solute concentration C e . Equilibrium mass transfer between the fluid and solid phases is described by a combination of Freundlich equilibrium expressions xe ¼ xef þ xes ;

ð49Þ

xe ¼ Kf ðC e Þnf þ Ks ðC e Þns ;

ð50Þ

where Kf and Ks are Freundlich sorption capacity coefficients for the rapidly sorbing and slowly sorbing solid-phase fraction, respectively; and nf and ns are Freundlich exponents for the rapidly sorbing and slowly sorbing fraction, respectively. The formulation summarized by Eqs. (43)–(50) is a coupled system of nonlinear PDEs in two primary dependent variables, C and xs , if xf is expressed as a function of C by invoking the assumption of local equilibrium. Eq. (43) is frequently written in nonconservative form (e.g. [14,26]) by assuming local equilibrium for the fast sites and using the chain rule to yield Rf

ð44Þ

oC o2 C oC q ¼ Dx 2  vx  ka C ra  b ½aðxes  xs Þ þ kf xrf f ; h ot ox ox ð51Þ in ½0; Lx ½0; ts ;

subject to boundary conditions Cð0; tÞ ¼ C0 ;

ð45Þ

oC ðLx ; tÞ ¼ 0 ox

ð46Þ

and initial conditions  px Cðx; 0Þ ¼ sin ; Lx xs ðx; 0Þ ¼ 0; 1

Rf ¼ 1 þ

qb oxf ; h oC

ð52Þ

ð47Þ

where Rf is termed the retardation factor. We use this form of the equation in this work. In terms of the operator notation introduced Section 3, we have q Rf ðCÞ ¼ 1 þ b Kf nf C nf 1 ; ð53Þ h

ð48Þ

Lt ðCÞ ¼ Dx

where ½0; Lx  R is the spatial domain; ½0; ts is the temporal domain; x and t are the spatial and temporal coordinates, respectively; C is an aqueous phase solute concentration; t is time; Dx is a hydrodynamic dispersion coefficient; vx is a mean pore velocity; ka , kf , and ks are reaction rate constants for the aqueous-phase, the rapidly sorbing solid phase, and the slowly sorbing solid phase, respectively; ra , rf , and rs are reaction order exponents for the aqueous phase, the rapidly sorbing solid phase, and the slowly sorbing solid phase, respectively; qb is a solid-phase bulk density; a is a first-order mass transfer rate coefficient for mass transfer between the aqueous phase and the slowly sorbing solid phase; h is porosity; x is a mass averaged solid-phase solute con-

o2 C oC ;  vx 2 ox ox

Lr1 ðC; xs Þ ¼ ka C ra 

qb ½aðKs C ns  xs Þ þ kf C rf ; h

Lr2 ðC; xs Þ ¼ aðKs C ns  xs Þ  ks xrss :

ð54Þ ð55Þ ð56Þ

4.4.2. Numerical solution We implemented the ISO algorithm using an approach in which the transport and reaction substeps are both solved using an ODE solver with formal error control. The transport substep is solved using a method of lines (MOL) approach in which the spatial derivatives are approximated using a centered finite difference to produce an ODE system for the concentrations at the

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

nodes. This ODE system is integrated over ½0; Dt to a prescribed error tolerance using the stiff ODE solver ODE15s. The reaction substep consists of solving a local reaction problem over ½0; Dt at each node, again using ODE15s to produce a solution with a prescribed level of error. We implemented several methods for computing the reaction operator estimate e L r1 used in the transport substep, and the transport operator estimate e L t used in the reaction substep. In the simplest approach these estimates were held constant over the respective substeps. Before the transport solve, the reaction operator estimate e L r1 is computed using the current approximate solution. It is then kept constant during the transport solve. Before the reaction solve, the transport operator estimate e L t is computed using the currrent approximate solution. It is then kept constant during the reaction solve. We also implemented two approaches which represented e L t and e L r1 as functions of time. One approach uses linear interpolation between the current approximate solution and the solution at the begining of the time step. Another approach uses Hermite interpolation. Based on the results in Section 4.3 for the nonLipschitz Rf , we considered a technique that uses a smooth approximation of the Freundlich isotherm. The n approach approximates xef ¼ Kf ðC e Þ f with a cubic spline, using about 1000 points to cover the interval C 2 ½0; C0 . Then Rf was computed using the derivative of the spline representation of xef .

257

4.4.3. Numerical experiments In this instance no analytical solution is available. Therefore, our exact solution is a numerical one computed with a MOL approach using the stiff ODE solver ODE15s. We selected this approach because it allows us to control the temporal error. We use this exact solution to compute the error in the ISO solution at each iteration. We compare the convergence over a single time step. Table 14 shows the parameter set used in the model problem. Numerical experiments were conducted for a single time step over a wide range of spatial and temporal discretizations. We considered the effect of interpolation accuracy for e L t and e L r1 on convergence by using constant, linear, and Hermite interpolation. We observed results similar to those for the simple ODE problem, so they will not be shown here. In summary, we found that constant and linear interpolation both degraded the convergence rate to OðDtÞ, while Hermite interpolation of the operator estimates preserved the OðDt2 Þ convergence rate predicted by the analysis. Tables 15–18 show the errors and convergence rate for the unsplined and splined Rf approaches, respectively, for a 50 cell spatial discretization. Hermite interpolation was used in both approaches, and Dt ¼ 1=2k for k ¼ 6; 7; . . . ; 11. The exact solution and the substeps of the ISO solution were solved in MATLAB using the ODE15s code with relative and absolute error tolerances set at 1013 . We observe that using the

Table 14 Model parameters Dx

vx

qb

h

a

Kf

Ks

ka

kf

0.01

1

1.59

0.4

0.1

0.126

0.252

1

1

nf

ns

ra

rf

rs

Lx

ts

C0

C

0.7

0.7

1

1

1

1

1

1

0

ks 1 0

x0s 0

Table 15 Model PDE with unsplined Rf approach: Errors Ei i n Dt

1.5625e02

7.8125e03

3.9062e03

1.9531e03

9.7656e04

4.8828e04

0 1 2 3

1.3289e01 6.1164e03 1.6587e03 1.6622e03

6.8916e02 1.1270e03 2.4611e04 2.4636e04

3.4932e02 2.1160e04 3.6634e05 3.7587e05

1.7561e02 4.1161e05 6.7019e06 6.7780e06

8.8004e03 8.1896e06 1.2379e06 1.2441e06

4.4048e03 1.6388e06 2.2518e07 2.2566e07

Table 16 Model PDE with unsplined Rf approach: Estimated rates Qi ðDtÞ=Qi ð2DtÞ i n Dt

7.8125e03

3.9062e03

1.9531e03

9.7656e04

4.8828e04

1 2 3

3.5530e01 8.0524e01 9.9892e01

3.7041e01 7.9281e01 1.0250eþ00

3.8696e01 9.4046e01 9.8570e01

3.9702e01 9.2836e01 9.9368e01

3.9980e01 9.0900e01 9.9723e01

258

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

Table 17 Model PDE with splined Rf approach: Errors Ei i n Dt

1.5625e02

7.8125e03

3.9062e03

1.9531e03

9.7656e04

4.8828e04

0 1 2 3

1.3289e01 1.4244e03 1.6587e03 1.6622e03

6.8916e02 2.2490e04 2.4611e04 2.4636e04

3.4932e02 3.1539e05 3.5803e05 3.5812e05

1.7561e02 4.1739e06 4.8263e06 4.8266e06

8.8004e03 5.3679e07 6.2645e07 6.2646e07

4.4048e03 6.8058e08 7.9793e08 7.9794e08

Table 18 Model PDE with splined Rf approach: Estimated rates Qi ðDtÞ=Qi ð2DtÞ i n Dt

7.8125e03

3.9062e03

1.9531e03

9.7656e04

4.8828e04

1 2 3

3.0445e01 9.3971e01 9.9892e01

2.7666e01 1.0374eþ00 9.9924e01

2.6326e01 1.0186eþ00 9.9980e01

2.5662e01 1.0093eþ00 9.9995e01

2.5331e01 1.0046eþ00 9.9999e01

splined Rf approach preserves the OðDt2 Þ convergence rate.

5. Discussion The analysis demonstrates that under fairly reasonable smoothness assumptions, one can show that the theoretical convergence rate for ISO applied to our model NRTPs is OðDt2 Þ, where Dt is the time step of the numerical method. This result is independent of the class of spatial discretization employed and is valid in 1D, 2D, or 3D domains. The theory does assume that the relevant temporal integration of the split equations can be done as accurately as needed. Obviously, this assumption must be relaxed in practice and may result in slower convergence in real applications. Numerical experiments on a simple ODE system indicate that several factors may reduce the convergence rate of ISO in practice. We observed that reducing the accuracy used in the integration of the split equations has an adverse effect on the rate of convergence. We also observed that reducing the accuracy of the interpolation used in the operator estimates can lower the convergence rate. Numerical experiments using the simple ODE also demonstrate that in order to preserve the second-order convergence of the ISO algorithm, one must have accurate estimates for the reaction and transport operators in the transport solve and reaction solve, respectively. Our results for the simple ODE showed that linear interpolation effects were seen after one iteration. Results not tabulated here showed that using constants for these operator estimates completely destroyed the secondorder convergence for the simple ODE. Since using a constant approximation of the these operators is a typical approach [45,62,74], it has significant implications for solving realistic problems. Numerical experiments using an ODE with a model Rf term lend insight into the importance of the Lipschitz

continuity assumption used in the analysis. The analysis treats Lipschitz continuity as a sufficient condition, but the numerical experiments indicate that it may be a necessary condition. This has significant practical implications for subsurface transport investigations, since the Freundlich model for sorption equilibrium often used in subsurface transport simulations results in an Rf which is non-Lipschitz for very low concentrations. The Freundlich model does not reduce to a linear relationship between C and x at low fluid-phase concentrations. Instead, ox=oC tends to infinity as the fluid phase solute concentration approaches zero. The results of the numerical experiments on the model NRTP show that one can implement the ISO algorithm in a manner that preserves the second-order convergence predicted by the theory. Simple constant or linear approximations for the transport and reaction operator estimates destroy the second-order convergence rate, while Hermite interpolation along with a smooth approximation for the retardation factor preserve second-order convergence. We stress that we looked at convergence rate, and we have not systematically investigated the computational efficiency of this approach. If the work for each iteration of the ISO method is excessive, it may not be competitive with other SO approaches or with the best FC approaches. This is an area which deserves further research.

6. Conclusions Based upon the analysis and numerical experiments described in the previous sections, we draw several conclusions regarding the use of the ISO algorithm for NRTPs: (1) Theoretical analysis shows that, given certain assumptions regarding smoothness, the rate of convergence for ISO solution of NRTPs is OðDt2 Þ.

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

(2) Numerical experiments using a simple ODE system indicate that the theoretical convergence rate can be degraded by discretization errors introduced into the numerical solution and by coarse estimates of the reaction and transport operators. Experiments on ODE models also imply that Lipschitz continuity is an important assumption in practice. (3) Numerical experiments on a model PDE system using typical finite difference implementations of the ISO approach (i.e., with constant and linear estimates of the reaction and transport operators) and using the Freundlich sorption equilibrium model agree with the results of ODE experiments. (4) Our model PDE results show that by using Hermite interpolation to compute the transport and reaction operator estimates, and using a smooth approximation to the Freundlich sorption model, one can obtain the theoretical rate of convergence. (5) Additional research on the computational efficiency of the the ISO approach introduced here is needed before it is proposed for routine subsurface simulation projects.

Acknowledgements The UNC aspects of this work have been supported by the National Institute of Environmental Health Sciences (NIEHS) grant 5 P42 ES05948-02 and the National Science Foundation (NSF) grant DMS-0112653. The work at NCSU has been supported by NSF grants DMS-0070641 and DMS-0112542, and the US Army Research Office (ARO) grant DAAD19-99-1-0186. Computational support was provided by the North Carolina Supercomputer Center.

References [1] Abrams RH, Loague K. A compartmentalized solute transport model for redox zones in contaminated aquifers, 2. Field-scale simulations. Water Resour Res 2000;36(8):2015–29. [2] Appelo CAJ, Willemsen A. Geochemical calculations and observations on salt water intrusions, I. A combined geochemical/ mixing cell model. J Hydrol 1987;94:313–30. [3] Avci AK, Trimm DL, Onsan ZI. Heterogeneous reactor modeling for simulation of catalytic oxidation and steam reforming of methane. Chem Eng Sci 2001;56(2):641–9. [4] Barry DA. Supercomputers and their use in modeling subsurface solute transport. Rev Geophys 1990;38(3):277–95. [5] Barry DA, Miller CT, Culligan PJ, Bajracharya K. Split operator methods for reactive chemical transport in groundwater. In: Proceeding of the International Conference on Modeling and Simulation 1995: MODSIM95. New South Walets, 1995. Modeling and Simulation Society of Australia, Inc., p. 125–32. [6] Barry DA, Miller CT, Culligan PJ, Bajracharya K. Analysis of split operator methods for nonlinear and multispecies groundwater chemical transport models. Math Comput Simulat 1997;43:331–41.

259

[7] Barry DA, Miller CT, Culligan-Hensley PJ. Temporal discretization errors in non-iterative split-operator approaches to solving chemical reaction/groundwater transport models. J Contam Hydrol 1996;22(1/2):1–17. [8] Barry DA, Bajracharya K, Crapper M, Prommer H, Cunningham CJ. Comparison of split-operator methods for solving coupled chemical non-equilibrium reaction/groundwater transport models. Math Comput Simulat 2000;53:113–27. [9] Bey I, Jacob DJ, Yantosca RM, Logan JA, Field BD, Fiore AM, et al. Global modeling of tropospheric chemistry with assimilated meteorology: Model description and evaluation. J Geophys ResAtmosph 2001;106(D19):23073–95. [10] Bhuyan D, Lake LW, Pope GA. Mathematical modeling of highpH chemical flooding. SPE Reservoir Eng 1990;5:213–20. [11] Blom JG, Verwer JG. A comparison of integration methods for atmospheric transport-chemistry problems. J Computat Appl Math 2000;126(1–2):381–96. [12] Blunt M, Fayers FJ, Orr FM. Carbon-dioxide in enhanced oilrecovery. Energy Convers Mgmt 1992;34(9–1):1197–204. [13] Borden RC, Bedient PB, Lee MD, Ward CH, Wilson JT. Transport of dissolved hydrocarbons influenced by oxygenlimited biodegradation, 2. Field application. Water Resour Res 1986;22(13):1983–90. [14] Brusseau ML, Rao PSC. Sorption nonideality during organic contaminant transport in porous media. Crit Rev Environ Contr 1989;19(1):33–99. [15] Celia MA, Kindred JS, Herrera I. Contaminant transport and biodegradation 1. A numerical model for reactive transport in porous media. Water Resour Res 1989;25(6):1141–8. [16] Cheng HP, Yeh GT. Development and demonstrative application of a 3-D numerical model of subsurface flow, heat transfer, and reactive chemical transport: 3DHYDROGEOCHEM. J Contam Hydrol 1998;34(1–2):47–83. [17] Chiang CY, Dawson CN, Wheeler MF. Modeling of in-situ biorestoration of organic compounds in groundwater. Transport Por Med 1991;6:667–702. [18] Chrysikopoulos CV, Kitanidis PK, Roberts PV. Analysis of onedimensional solute transport through porous media with spatially variable retardation factor. Water Resour Res 1990;26(3):437–46. [19] Chu SP, Elliott H. Latitude versus depth simulations of ecodynamics and dissolved gas chemistry relationships in the central pacific. J Atmosph Chem 2001;40(3):305–33. [20] Cooper DM, Jenkins A, Skeffington R, Gannon B. Catchmentscale simulation of stream water chemistry by spatial mixing: Theory and application. J Hydrol 2000;233(1–4):4121–37. [21] Dawson CN, Van Duijin CJ, Wheeler MF. Characteristic-galerkin methods for contaminant transport with nonequilibrium adsorption kinetics. SIAM J Num Anal 1994;31(4):982–9. [22] Dawson CN, Wheeler MF. Time-splitting methods for advection– diffusion–reaction equations arising in contaminant transport. In: OÕMalley R, editor. Proceedings of the ICIAM Õ91. SIAM; 1992. p. 71–82. [23] Dykaar BB, Kitanidis PK. Macrotransport of a biologically reacting solute through porous media. Water Resour Res 1996;32(2):307–20. [24] Engesgaard P, Kipp KL. A geochemical transport model for redoxcontrolled movement of mineral fronts in groundwater flow systems A case of nitrate removal by oxidation of pyrite. Water Resour Res 1992;28(10):2829–43. [25] Engesgaard P, Traberg R. Contaminant transport at a waste residue deposit. 2. Geochemical transport modeling. Water Resour Res 1996;32(4):939–51. [26] Harmon TC, Ball WP, Roberts PV. Nonequilibrium transport of organic contaminants in groundwater. In: Sawhney BL, Brown K, editors. Reactions and Movement of Organic Chemicals in Soils. Madison, Wisconsin: Soil Science Society of America, Inc.; 1989. p. 405–38.

260

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

[27] Herzer J, Kinzelbach W. Coupling of transport and chemical processes in numerical transport models. Geoderma 1989;44:115– 27. [28] Hunter KS, Wang YF, Van Cappellen PM. Kinetic modeling of microbially-driven redox chemistry of subsurface environments: Coupling transport, microbial metabolism and geochemistry. J Hydrol 1998;209(1–4):53–80. [29] Jacobson MZ. Fundamentals of atmospheric modeling. Cambridge University Press; 1998. [30] Jennings AA, Kirkner DJ, Theis TL. Multicomponent equilibrium chemistry in groundwater quality models. Water Resour Res 1982;18(4):1089–96. [31] Kaluarachchi JJ, Morshed J. Critical assessment of the operator splitting technique in solving the advection dispersion reaction equation.1. First-order reaction. Adv Water Resour 1995;18(2):89–100. [32] Kanney JF, Miller CT, Barry DA. Comparison of fully coupled approaches for approximating nonlinear transport and reaction problems. Adv Water Resour, in press. [33] Kim J, Cho SY. Computation accuracy and efficiency of the timesplitting method in solving atmospheric transport-chemistry equations. Atmosph Environ 1997;31(15):2215–24. [34] Koelmans AA, Van der Heijde A, Knijff LM, Aalderink RH. Integrated modelling of eutrophication and organic contaminant fate and effects in aquatic ecosystems. A review. Water Res 2001; 35(15):3517–36. [35] Kopmann R, Markofsky M. Three-dimensional water quality modelling with TELEMAC-3D. Hydrolog Process 2000;14(13): 2279–92. [36] Koziol AS, Pudykiewicz JA. Global-scale environmental transport of persistent organic pollutants. Chemosphere 2001;45(8): 1181–200. [37] Krishnaswamy S, Gunn RD, Agarwal PK. Low-temperature oxidation of coal, 2. An experimental and modelling investigation using a fixed-bed isothermal flow reactor. Fuel 1996;75(3):344–52. [38] Lake LW. Enhanced oil recovery. Englewood Cliffs, NJ: Prentice Hall; 1989. [39] Lanser D, Verwer JG. Analysis of operator splitting for advection-diffusion-reaction problems from air pollution modelling. J Computat Appl Math 1999;111(1–2):201–16. [40] LeVeque RJ, Oliger J. Numerical methods based on additive splittings for hyperbolic partial differential equations. Math Computat 1983;40(162):469–97. [41] Linn DM. Sorption and degradation of pesticides and organic chemicals in soil. Soil Science Society of America, Inc., and American Society of Agronomy, Inc., Madison, WI, 1993. [42] Liu CW, Narasimhan TN. Redox-controlled multiple-species reactive chemical transport 1. Model development. Water Resour Res 1989;25(5):869–82. [43] MacQuarrie KT, Sudicky EA, Frind EO. Simulation of biodegradable organic contaminants in groundwater, 1. Numerical formulation in principal directions. Water Resour Res 1990; 26(2):207–22. [44] Mahadevan A, Archer D. Modeling the impact of fronts and mesoscale circulation on the nutrient supply and biogeochemistry of the upper ocean. J Geophys Res-Oceans 2000; 105(C1):1209–25. [45] Mcnab WW, Narasimhan TN. A multiple species transport model with sequential decay chain interactions in heterogeneous subsurface environments. Water Resour Res 1993;29(8):2737–46. [46] McRae GJ, Goodin WR, Seinfeld JH. Numerical solution of the atmospheric diffusion equation for chemically reacting flows. J Computat Phys 1982;45:1–42. [47] Miller CT, Christakos G, Imhoff PT, McBride JF, Pedit JA, Trangenstein JA. Multiphase flow and transport modeling in heterogeneous porous media: Challenges and approaches. Adv Water Resour 1998;21(2):77–120.

[48] Miller CT, Rabideau AJ. Development of split-operator, Petrov– Galerkin methods to simulate transport and diffusion problems. Water Resour Res 1993;29(7):2227–40. [49] Mizsey P, Cuellar A, Newson E, Hottinger P, Truong T, von Roth F. Fixed bed reactor modelling and experimental data for catalytic dehydrogenation in seasonal energy storage applications. Comput Chem Eng 1999;23:S379–82. [50] Morshed J, Kaluarachchi JJ. Critical assessment of the operator splitting technique in solving the advection dispersion reaction equation 2. Monod kinetics and coupled transport. Adv Water Res 1995;18(2):101–10. [51] Mossman DJ, AlMulki N. One-dimensional unsteady flow and unsteady pesticide transport in a reservoir. Ecologl Modell 1996;89(1–3):259–67. [52] Oran ES, Boris JP. Numerical simulation of reactive flow. New York: Elsevier; 1987. [53] Parkhurst DL. UserÕs Guide to PHREEQC: A computer program for speciation, reaction-path, advective-transport and diffusion problems. US Geological Survey, 1995. Water Resource Investigation Report 95-4227. [54] Pedit JA. An investigation of heterogeneous sorption by a natural solid, PhD thesis. Chapel Hill, NC: University of North Carolina; 1994. [55] Pope GA, Baviere M. Basic concepts in enhanced oil recovery processes. New York: Elsevier; 1991. [56] Prommer H, Barry DA, Davis GB. A one-dimensional reactive multi-component transport model for biodegradation of petroleum hydrocarbons in groundwater. Environ Modell Software 1999;14:213–23. [57] Prommer H, Davis GB, Barry DA. PHT3D––A three-dimensional biogeochemical transport model for modelling natural and enhanced remediation. In: Johnston CD, editor. Contaminated site remediation: challenges posed by urban and industrial contaminants. Fremantle, WA, 21–25 March 1999. p. 351– 58. [58] Rocha F, Walker A. Simulation of the persistence of atrazine in soil at different sites in portugal. Weed Res 1995;35:179–86. [59] Rubin J. Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions. Water Resour Res 1983;19(5):1231– 52. [60] Runkel RL, Bencala KE, Broshears RE, Chapra SC. Reactive solute transport in streams. 1. Development of an equilibriumbased model. Water Resour Res 1996;32(2):409–18. [61] Saaltink MW, Carrera J, Ayora C. On the behavior of approaches to simulate reactive transport. J Contam Hydrol 2001;48:213–35. [62] Schafer D, Schafer W, Kinzelbach W. Simulation of reactive processes related to biodegradation in aquifers. 1. Structure of the three-dimensional reactive transport model. J Contam Hydrol 1998;31(1-2):167–86. [63] Scheibe TD, Chien YJ, Radtke JS. Use of quantitative models to design microbial transport experiments in a sandy aquifer. Ground Water 2001;39(2):210–22. [64] Shampine LF, Reichelt MW. The MATLAB ODE suite. SIAM J Sci Comput 1997;18(1):1–22. [65] Simunek J, Suarez DL. Two-dimensional transport model for variably saturated porous media with major ion chemistry. Water Resour Res 1994;30(4):1115–33. [66] Singh HB. Composition, chemistry, and climate of the atmosphere. New York: John Wiley and Sons; 1995. [67] Sportisse B, Bencteux G, Plion P. Method of lines versus operator splitting for reaction–diffusion systems with fast chemistry. Environ Modell Software 2000;15(6–7):673–9. [68] Steefel CI, MacQuarrie KTB. Approaches to modeling of reactive transport in porous media. In: Lichtner PC, Steefel CI, Oelkers EH, editors. Reactive transport in porous media: general princi-

J.F. Kanney et al. / Advances in Water Resources 26 (2003) 247–261

[69]

[70]

[71] [72] [73]

[74]

[75]

[76]

ples and applications to geochemical processes. Washington, DC: Mineralogical Society of America; 1996. p. 83–129. Steefel CI, Yabusaki SB. OS3D/GIMRT, software for modelling multicomponent-multidimensional reactive transport, user manual and programmerÕs guide. Richland, Washington: Pacific Northwest Laboratory; 1995. Stefanovic DL, Stefan HG. Accurate two-dimensional simulation of advective–diffusive–reactive transport. J Hydraulic Eng-ASCE 2001;127(9):728–37. Strang G. On the construction and comparison of difference schemes. SIAM J Num Anal 1968;5(3):506–17. Sun N-Z. Mathematical modeling of groundwater pollution. New York: Springer-Verlag; 1996. Valocchi AJ, Malmstead M. Accuracy of operator splitting for advection–dispersion–reaction problems. Water Resour Res 1992;28(5):1471–6. Walter AL, Frind EO, Blowes DW, Ptacek CJ, Molson JW. Modeling of multicomponent reactive transport in groundwater 1. Model development and evaluation. Water Resour Res 1994;30(11):3137–48. Wang KY, Pyle JA, Shallcross DE. Formulation and evaluation of Ims, an interactive three-dimensional tropospheric chemical transport model 1. Model emission schemes and transport processes. J Atmosph Chem 2001;38(2):195–227. Wheeler MF, Dawson CN. An operator-splitting method for advection–diffusion–reaction problems. Technical Report 87-9.

[77]

[78] [79]

[80]

[81]

[82]

[83]

261

Rice University Department of Computational and Applied Mathematics, Houston, TX, 1987. Wheeler MF, Dawson CN, Bedient PB, Chiang CY, Borden RC, Rifai HS. Numerical simulation of microbial biodegradation of hydrocarbons in ground water. In: Proceedings of the NWWA Conference on Solving Ground Water Problems with Models. 1987. p. 92–109. Wood TM, Baptista AM. A model for diagnostic analysis of estuarine geochemistry. Water Resour Res 1993;29(1):51–71. Xu T, Samper J, Ayora C, Manzano M, Custidio E. Modeling of non-isothermal multi-component reactive transport in field scale porous media flow systems. J Hydrol 1999;214(1-4):144– 64. Yeh GT, Tripathi VS. A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components. Water Resour Res 1989;25(1):93–108. Yeh G-T, Tripathi VS. A model for simulating transport of reactive multispecies components: Model development and demonstration. Water Resour Res 1991;27(12):3075–94. Zysset A, Stauffer F. Modeling of microbial processes in groundwater infiltration systems. In: Russell TF, Ewing RE, Brebbia CA, Gray WG, Pinder GF, editors. Mathematical modeling in water resources, vol. 2. Southampton, UK: Computational Mechanics Publications; 1992. p. 325–32. Zysset A, Stauffer F, Dracos T. Modeling of chemically reactive groundwater transport. Water Resour Res 1994;30(7):2217–28.