The control volume radial basis function method CV-RBF with Richardson extrapolation in geochemical problems

Computers & Geosciences 76 (2015) 151–163

Contents lists available at ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

The control volume radial basis function method CV-RBF with Richardson extrapolation in geochemical problems W.F. Florez a,n, M. Portapila c, A.F. Hill a, H. Power b, P. Orsini b, C.A. Bustamante a a

School of Engineering, Pontifical Bolivariana University, U.P.B., Circ. 1, 70-01 Medellin, Colombia School of Mechanical, Materials and Manufacturing Engineering, The University of Nottingham, Nottingham NG7 2RD, UK c CIFASIS, University of Rosario, Ocampo y Esmeralda, S2000EZP Rosario, Argentina b

art ic l e i nf o

a b s t r a c t

Article history: Received 2 September 2014 Received in revised form 10 December 2014 Accepted 2 January 2015 Available online 5 January 2015

The aim of this paper is to present how to implement a control volume approach improved by Hermite radial basis functions (CV-RBF) for geochemical problems. A multi-step strategy based on Richardson extrapolation is proposed as an alternative to the conventional dual step sequential non-iterative approach (SNIA) for coupling the transport equations with the chemical model. Additionally, this paper illustrates how to use PHREEQC to add geochemical reaction capabilities to CV-RBF transport methods. Several problems with different degrees of complexity were solved including cases of cation exchange, dissolution, dissociation, equilibrium and kinetics at different rates for mineral species. The results show that the solution and strategies presented here are effective and in good agreement with other methods presented in the literature for the same cases. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Hermite radial basis functions RBF Richardson extrapolation Chemical equilibrium Kinetics Groundwater flow

1. Introduction The radial basis functions (RBFs) have been widely used to solve partial differential equations either on their own as a pointwise meshless method (Jumarhon et al., 2000) or combined with other traditional methods to improve their accuracy (Orsini et al., 2009; Wright and Fornberg, 2006; Sladek et al., 2005). The RBF interpolation method is considered as an optimal numerical technique for interpolating multidimensional scattered data (Powell, 1994; Madych and Nelson, 1990). However, the application of RBF transcends the interpolation because they can also be used as the basis of several meshless collocation approaches for solving partial differential equations (PDEs) (see Kansa, 2000, and Jumarhon et al., 2000). The most tested approaches so far are global methods that approximate the differential operator over the whole domain (Bustamante et al. , 2013), local methods that approximate the PDE locally on a stencil made of some nodes (Bustamante et al., 2014), unsymmetrical methods that directly approximate the solution variable (Kansa, 2000), symmetric methods that achieve a Hermite interpolation with the RBF and the differential operator (Stevens et al., 2009) and finally some methods that approximate the particular solution by RBF (Chen et al., 2012; Bustamante et al., 2014). Although the global formulation becomes unpractical when n

Corresponding author. E-mail address: whady.fl[email protected] (W.F. Florez).

http://dx.doi.org/10.1016/j.cageo.2015.01.001 0098-3004/& 2015 Elsevier Ltd. All rights reserved.

the number of collocation points is relative large, the local implementation can be used for the improvement of classical numerical methods such as the control volume method CV (Versteeg and Malalasekera, 2007). Specifically, a conservative RBF interpolation (Orsini et al., 2008) can be used to improve the flux calculation on the faces of each volume (Ollivier-Gooch and Altena, 2002). There are certain specific advantages in using Hermite RBF for solving the transport differential equations, such as the cell surface fluxes are reconstructed using local collocation that besides interpolating the nodal values of the variable also satisfies the governing equation at some auxiliary points in each volume; the use of Hermite interpolation that approximates the governing equation is a form of analytical upwinding and finally, the boundary conditions are enforced explicitly on the interpolation formula at those volumes that have some nodes on the boundary and so there is no need for artificial node schemes. The above advantages are not included in traditional CV formulations (Zabka and Sembera, 2011). The idea of employing RBF interpolations to improve the accuracy of a classical numerical scheme has been proposed by Wright and Fornberg (2006). These authors worked with a Hermite RBF interpolation to remove the symmetry constraint required to achieve high-order approximation in the finite difference scheme. Nguyen-Van et al. (2007) introduce a mesh-free conforming nodal integration into the finite element method.

152

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

Regarding the CV method, Moroney and Turner (2006) and Liu et al. (2002) used RBF interpolation functions instead of Lagrange polynomial shape functions to reconstruct the field variables and their derivatives. Their approach relies on a local RBF interpolation of the field variable used to obtain the surface fluxes where the CV centers are considered as collocation points. Moroney and Turner ascribed the ability of the CV-RBF scheme to approximate the solution on relatively coarse meshes to the high accuracy of the RBF interpolation for evaluating derivatives (Fornberg and Flyer, 2005). For strongly advective problems, Turner et al. (Liu et al., 2002) RBF approach requires some kind of up-winding scheme in order to avoid spurious oscillations. It is important to point out that in the aforementioned work a classical RBF interpolation is employed. About the simulation of reactive transport in groundwater, Saaltink et al. (2001) declare that there are two types of approaches namely the Direct substitution Approach (DSA), based on the Newton–Raphson method and the Picard or Sequential Iteration approach (SIA). They concluded that DSA is more robust than the SIA, i.e. that its convergence is less sensitive to time step size, but due to the size of the problem the computation time and memory requirements tend to be less favorable. Xu et al. (1999) compared the numerical performance of the SIA approach and the sequential non-iterative approach (SNIA), and they concluded that the accuracy of SNIA depends on space and time discretization as well as on the nature of the chemical reaction. The SNIA method does not solve the actual coupled problem of transport and chemical reactions, but a sequence of solution based on operator splitting. Nevertheless, Samper et al. (2009) propose that the accuracy of the traditional SNIA can be improved by iteration between transport and chemical equations only in those areas where there is a large mass transfer between solid and liquid phases. These authors also point out that the SNIA method could produce larger errors for cation exchange problems. Kanney et al. (2003) mention that each time step in the SNIA method may include one or more steps, e.g. the reaction step may be solved using multiple steps. Carrayrou et al. (2004) stated that the classical SNIA approach can be modified by using two time steps. Nardi et al. (2014) presented the development, verification and application of an efficient interface based on SNIA to couple two standalone simulation programs: the Finite Element framework COMSOL and the geochemical simulator PHREEQC. Also, they emphasize on numerical robustness, efficiency and parallel implementation for 3D cases. Wissmeier and Barry (2011) showed that an advantage of using PHREEQC to couple with a diffusion convection code is that the soil solution is speciated from its element composition according to thermodynamic mass action equations. Therefore, the element concentrations after using PHREEQC can be passed to each cell node on the transport equation domain. Although, in their work they solved some 2D problems for irrigation, they do not present any adaptive sequential reduction of the step to improve SNIA. Kolditz et al. (2012) used OpenGeoSys, a complete self-contained software/environment for numerical simulation of thermo-hydromechanical-chemical processes in porous media, based on the Finite Element Method (FEM) for solving multifield problems in porous media. Zabka and Sembera (2011) present a one-dimensional advection–diffusion model based on the Finite Volume Method. Their results are compared with experimental data. Mugler et al. (2004) proposed to combine 3D transport codes such as CASTEM and MT3D with some available chemical packages like PHREEQC and CHESS, to obtain a software platform for the simulation of nuclear waste storage. The coupling algorithm they used is sequential iterative SIA and not SNIA. In SIA the concentrations are calculated by solving successively and iteratively until convergence of the transport equation and the chemical equation.

Nasir et al. (2014) used PHREEQC to solve geochemical reactions while the transport equations are solved by COMSOL, but the coupling between transport and chemistry is made through MATLAB. Their code includes five coupled processes: thermal and mechanical processes, saturated flow, solute transport and chemical reactions. The obtained numerical results include the evolution of the rock porosity and permeability which are important for nuclear waste. Long term processes are also simulated and the results are compared with a couple of known benchmark solutions. Patel et al. solved the reaction step by PHREEQC and their contribution was to link the Lattice–Boltzmann method (LB) with the external geochemical solver. Hiorth et al. (2013) derived the LB boundary condition for nonlinear surface reaction kinetics. The chemical model they used is based on a state equation and in the SUPCRT geochemical database implemented in software EqAlt. They showed that the LB chemical simulations are nearly identical to PHREEQC and they illustrated their LB model by simulating the alteration produced when seawater is injected into chalk. In this paper the CV-RBF scheme proposed by Orsini et al. (2009), where the cell flux on the control volumes is reconstructed by a local Hermite interpolation, is extended for the first time, to study reactive groundwater flow problems with chemical equilibrium, kinetics and cation exchange. Currently, the feasibility of using local RBF methods for this kind of problems has not been assessed. Thus, it is the aim of this paper to solve some benchmark one dimensional problems to show how the RBF methods perform in cases with reactive groundwater flow. Also, this work introduces an improvement to the traditional SNIA method for coupling transport and reaction by means of a high order Richardson extrapolation. Additionally, this paper illustrates how to use the PHREEQC subroutines (Charlton and Parkhurst, 2011) to add geochemical reaction capabilities to groundwater transport models with an operator splitting scheme. The numerical solution for three standard problems are presented with several degrees of complexity: the groundwater flow with cation exchange in chemical equilibrium, the chemical kinetic flows with calcite dissolution and the flow with cation exchange together with equilibrium or kinetics. All of the numerical results were compared with other methods from the literature, showing the validity of both the proposed new Richardson high order method and the transport scheme. The obtained data shows that the CV-RBF method for this type of geochemical problems performs at least as well as other traditional methods like finite elements, finite volumes and finite differences. Also, the proposed Richardson multi-step approach is a convergent and numerically efficient procedure to diminish the error and so, it becomes a valid alternative to simple SNIA methods.

2. RBF interpolation of the field variables In the recent years, the theory of RBF interpolation has been successfully applied to approximate variables and functions in several meshless methods and also in traditional methods that require local approximation of the variables. The radial function Ψ (∥ x − ξ j ∥) depends upon the separation distances of a subset of collocation centers {ξ j ∈ R n; j = 1, 2, … , N } and a field point x ∈ R n where N is the number of collocation points. The name radial comes from the RBF spherical symmetry around the centers ξ j , and the distances ∥ x − ξ j ∥ are usually the Euclidean metric. In the RBF interpolation it is usual to select the collocation nodes and the field nodes as the same set of points; however, this is not necessary in principle. The most used RBFs applied to the solution of differential equations so far are listed in Table 1, where M is an positive integer and r = ∥ x − ξ j ∥ is the Euclidean distance

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

153

be interpolated, as follows:

Table 1 Table of common radial basis functions.

N

Thin plate spline

Generalized multiquadric

Gaussian

r 2M − 2 log r

(r 2 + c 2) M /2

exp( − r /c)

ϕ (x) =

NP

∑ α j Ψ (∥ x − ξ j ∥) + ∑ α j + N P Mj − 1(x) j=1

j=1

(1)

where some polynomial terms have been included, and they have to satisfy some algebraic constraints (Orsini et al., 2009); αj with j = 1, … , N , N + 1, … , N + NP , are real coefficients to be found by solving a collocation system of linear equations as explained in Bustamante et al.(2013); Bustamante et al. (2014), Orsini et al. (2008) and Kansa (2000), Ψ is an RBF and NP is the number of polynomial terms PM − 1. Although Micchelli (1986) proved that the collocation matrix resulting from the RBF interpolation is always non-singular, its condition number might be large (Schaback, 1995) for the interpolation over the whole domain, and so a stencil should be used to circumvent this problem (Orsini et al., 2009; Bustamante et al., 2014). Additionally, the RBF approximation on a stencil can be used to reconstruct the derivatives locally. Regarding the application of RBF to partial differential equations, it is possible to approximate the solution of a boundary value problem (BVP) by an Hermite interpolation of the following form (Zongmin, 1992; Fasshauer, 1997): nb

ϕ (x) =

N

∑ α j Bξ Ψ (∥ x − ξ j ∥) + ∑ j=1

α j L ξ Ψ (∥ x − ξ j ∥)

j = nb + 1

NP

+

∑ α j + N P Mj − 1(x)

(2)

j=1

Fig. 1. RBF interpolation on the faces of each control volume.

Table 2 Initial and boundary conditions for the cation exchange problem. Boundary conditions (units mmol/kg) pH pe Ca Cl

7.0 12.5 0.6 1.2

Initial conditions pH pe Na K N

7.0 12.5 1.0 0.2 1.2

where nb and N − nb are the number of boundary and internal nodes respectively, L ξ is the differential operator on the domain Ω and Bξ is operator on the boundary Γ associated with the BVP. After substitution of (2) into the boundary condition and into the differential equation, a symmetric and non-singular collocation linear system of equations is obtained (see Stevens et al., 2009; Zongmin, 1992; Pearson, 2013 for details). This system of equations is the starting point to define a local interpolation based on the governing equations on the faces of each control volume in the domain (Orsini et al., 2009; Stevens et al., 2009; Florez et al., 2012). The nodes around each control volume are used to construct a radial basis expansion that allows us to approximate the value of the variables and their derivatives on the CV faces e.g., the values at the nodes of each stencil can be used to interpolate the values on the boundaries, this process is outlined in Fig. 1.

3. CV-RBF formulation for the solution of the transport equations between node x and collocation points ξ. The Gaussian and the inverse multiquadric are positive-definite functions, while the thin-plate spline TPS and the multiquadric are conditionally positive-definite functions which require the addition of a polynomial along with an homogeneous constraint condition (Powell, 1994; Madych and Nelson, 1990) in order to obtain a non-singular interpolation matrix. Although any RBF could be used for the approximation of variables in a partial differential equation, they differ in their computational performance. For instance, the TPS functions only converge linearly (Powell, 1994) while the multiquadric (MQ) functions exhibit an exponential convergence (Madych and Nelson, 1990). Also, some RBFs have adjustable parameters such as c in the MQ functions. The setup of these parameters is not trivial and it still is an open problem (Wright and Fornberg, 2006). In a typical interpolation problem, there are N pairs of data points that are samples of the unknown function Φ to

The governing equation for the transport of a component in an nd-dimensional space can be written as

∂Ui ϕ ∂ϕ ∂ ⎛ ∂ϕ ⎞ ⎜⎜Dij ⎟⎟ − = , ∂t ∂xi ⎝ ∂x j ⎠ ∂xi

i, j = 1, nd (3)

where ϕ is a variable to be transported, D is the diffusivity and Uj is the advective velocity along the j-direction. The initial and boundary conditions are

ϕ (x, 0) = ϕ0

on Ω

B (ϕ) = g (x)

on Γ

(4)

By using an implicit time stepping scheme to discretize the time derivative and integrating on each control volume, the following integral equation is obtained:

154

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

Fig. 2. Comparison between CV-RBF and other benchmark solutions for the cation exchange problem.

Fig. 3. Effect of Richardson extrapolation on the cation exchange results.

VP ϕ¯t +Δt − Δt

∫S

⎛ ∂ϕt +Δt ⎞ ⎜⎜Dij − Ui ϕt +Δt ⎟⎟ ni dS = ϕ¯t VP ∂x j ⎝ ⎠

(5)

In the above equation, the bar indicates the cell average value over the control volume, and this value is considered to be located at the CV center. The volume of the cell is VP over which the integration is performed by the midpoint integration formula. The integral over the control volume surface can be divided into Ns integrals over each cell face, and each of these surface integrals can be evaluated by using the midpoint quadrature. After this approximation, the following discrete expression is obtained:

Ns ⎛ ⎞ ∂ϕt +Δt − Ui ϕt +Δt ⎟⎟ nil Sl = ϕ¯t VP VP ϕ¯t +Δt − Δt ∑ ⎜⎜Dij ∂ x j ⎠l l= 1 ⎝

(6)

In Eq. (6) the value of ϕ and its derivatives ∂ϕ/∂x j at the center of each cell face (⁎) l are obtained by the Hermite interpolation explained in Orsini et al. (2009), from the values at the collocation nodes surrounding each volume face. In contrast to other weak formulations such as finite elements (FEM) (Kolditz et al., 2012), the CV-RBF presents a strong formulation where the field variables are interpolated by Hermite radial basis functions. Therefore in CV-RBF the interpolation

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

155

Fig. 4. Effect of the Richardson extrapolation order on the cation exchange results.

Fig. 5. Effect of the number of cells for the cation exchange case.

functions include the differential operators i.e. the interpolation error of the approximation of the differential equations at the collocation nodes is zero. Additionally, the boundary conditions are not approximated but enforced exactly on the interpolation in each cell.

aqueous solution through a porous media such as underground water flow. The transport equation includes diffusion, advection and the reaction terms. In this paper the mass balance is stated for each element of the different species present in the fluid without including the generation term. The reactive transport equation of a total element concentration Ej in an aqueous solution has the following expression:

4. Transport and chemical reaction coupling for groundwater flow

∂E j

The following procedure is proposed in order to predict the concentration of the different species that are been added to a

where the L (E j ) is the lineal transport operator with diffusion and advection with velocity Ui:

∂t

= L (E j ) + R j

(7)

156

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

Fig. 6. RMS deviation as a function of the mesh density.

Table 3 RMS difference between CV-RBF and other benchmark solutions. Comparison

Na

PHREEQC

7.96 × 10−5

6.74 × 10−5

1.89 × 10−4

9.61 × 10−5

PHT3D

7.54 × 10−6

1.48 × 10−5

1.29 × 10−5

8.79 × 10−6

FD

6.69 × 10−4

1.11 × 10−4

5.55 × 10−4

7.72 × 10−4

L (E j ) =

Cl

∂ (Ui E j ) ∂ ⎛ ∂E j ⎞ ⎜⎜Dij ⎟⎟ − ∂xi ⎝ ∂x j ⎠ ∂xi

K

Ca

(8)

In general, the reaction term Rj depends on the concentration of individual species of the system. In geochemical problems the reaction term usually includes mineral precipitation, dissolution and sorption. As in Nardi et al. (2014), the mass balance equations are applied to the total concentration of each element and so no chemical sink-source terms for the liquid phase are necessary at the solute transport step. The solution of Eq. (7) requires to solve a coupled non-linear system of equations with numerical approaches such as Saaltink et al. (2001). A simpler alternative is to split the reactive and transport operators, and solve the two of them separately:

∂E j ∂t ∂E j ∂t

= L (E j )

= Rj

(9)

(10)

The solution of Eq. (9) requires the global solution of the transport equations by the CV-RBF method, i.e. considering the entire computational domain, while Eq. (10), even though involving non-linear calculations, it is only applied for each computational cell. As a result of the decoupling, concentrations are assumed independent of reactions in the solute transport computations (Wissmeier and Barry, 2011). In this paper two kinds of chemical cases are considered, problems in equilibrium and problems with kinetics. In order to solve the partial differential equations and to find the concentration of

the different elements with time, the transport phenomena are uncoupled from the kinetic and the equilibrium phenomena. The governing equations for the diffusive and convective phenomena are first solved by the CV-RBF with Hermite interpolation. Then, the local chemical equilibrium calculation in each cell stated in Eq. (10) is solved by means of PHREEQC (Charlton and Parkhurst, 2011). The geochemical subroutines PHREEQC can simulate a wide range of equilibrium reactions between water and minerals, ion exchangers, surface complexes, solid solutions, and gases. These subroutines also have a general kinetic formulation that allows modeling of non-equilibrium mineral dissolution and other kinetic reactions. As additional advantage, PHREEQC provides a framework to calculate other chemical variables that depend on the solution speciation such as pH, saturation indices, and ion activities (Wissmeier and Barry, 2011). PHREEQC can complement the CV-RBF transport code by adding geochemical reaction capabilities to the groundwater transport model. By doing so, it is possible to store and manipulate solution compositions and reaction information for any number of cells within the domain and to feed back this chemical data results to the transport model. The algorithm couples the transport diffusive and convective terms, of the transport equation, with the chemical kinetics by means of PHREEQC. One of the most widely used iterative methods for coupling the chemistry and transport equations is the sequential non-iterative approach SNIA (Nardi et al., 2014; Carrayrou et al., 2004). In this method, at each time step the solution of the transport equations (9) is followed by the solution of chemistry equations (10) before advancing to the next time. Usually, in equilibrium problems with SNIA, the reactions are assumed to be instantaneous compared to the transport time scale. For chemical kinetics problems the time step of the kinetic equations must be accounted for, thus the time advance is accomplished by a two step sequence: First, the concentrations Ej are obtained after solving the transport equations (9). That is, the solution values are fed to PHREEQC which implicitly solves Eq. (10) and updates the solution to the next time step. Specifically, the calculation of concentrations Ej is a two step algorithm where the transport equations are first solved from time t to t + Δt/2, and then the kinetics equations are solved to advance the solution from time t + Δt/2 to t. Because the chemical kinetics

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

157

does not really take place in a discrete two step process, it is convenient to improve the approximation by including more transport-reaction steps within the same time interval by using Richardson extrapolation. 4.1. Richardson extrapolation to improve the operator splitting SNIA method Richardson extrapolation is a method to accelerate the convergence of a sequence. If a method converges as O (Δt m), the error can be reduced by a factor of Δt to produce a technique that converges as O (Δt m + 1), and this can be used iteratively to produce other approaches that converge at any chosen rate. As a result, Richardson extrapolation can be used to obtain a higher order approximation from the sequence of transport-reaction steps in SNIA. Let Cj be the exact value of the concentrations that have a numerical approximation E j (Δt) that depends on the time step Δt , i.e

C j = E j (Δt) + c1Δt 2 + c2 Δt 3 + c3 Δt 4 + ... = E j (Δt) + O (Δt 2)

(11)

Eq. (11) means that E j (Δt) is accurate to O (Δt 2). Another numerical approximation can be obtained with a smaller step Δt/2:

C j = E j (Δt /2) + c1Δt 2/4 + c2 Δt 3/8 + c4 Δt 4/16 + ...

(12)

Note that in Eq. (12) the total interval Δt is divided into 2 steps Δt/2 in order to obtain the new numerical approximation E j (Δt /2) at the end of the total interval. Multiplying Eq. (12) by 4 and subtracting from the result, Eq. (11) gives

4 1 1 1 E j (Δt /2) − E j (Δt) − c2 Δt 3 − c3 Δt 4 + ... 3 3 6 4 4E j (Δt /2) − E j (Δt) = + O (Δt 3) 3

Cj =

(13)

Eq. (13) is a new approximation formula for the concentrations accurate to O (Δt 3), that is one order higher that the initial approximation (Eq. (11)). It can be seen that to obtain the new higher order approximation (Eq. (13)), the problem in a time interval must be solved twice, first using one step Δt and then with two smaller steps Δt/2. Also, it must be pointed out that each step in the solution sequence is made of a two stage transport-reaction calculation as was explained in Section 4. An even higher order approximation can be inferred following the same procedure but replacing Δt by Δt/2 in Eq. (13), to obtain the following expression:

Fig. 7. Dissolution of calcite.

Table 4 Initial and boundary conditions for calcite dissolution problems (cases CAL and WAD). Problem case

Primary Species

Initial concentrations (log mol l
1)

Influx concentrations (log mol l
1)

Secondary species

Rate constant (κ) (mol m
2 s
1)

s

Calcite dissolution

Hþ


7.978


5.496

OH−, CO2, CO23−, CaHCO+ 3

4.64 × 10−7

HCO− 3


3.018


5.421

Ca2 þ


3.019


4.398

6.8 × 10−5 (CAL-1) 6.8 × 10−4 (CAL-2) 6.8 × 10−3 (CAL-3)

H+


7.998


7.249

4.64 × 10−7

4.9 × 10−5 (WAD-1)

Na þ


0.320


7.000

Ca2 þ


2.022


2.770

Mg2 þ


1.323


7.000

Cl



0.247
5.175


6.700
5.456

6.8 × 10−2 (CAL-4) Calcite dissolution with cation exchange

CO2− 3

OH−, HCO− 3 , CO2 , XNa , X 2Ca, X2 Mg

4.9 × 10−4 (WAD-2) 4.9 × 10−3 (WAD-3) 4.9 × 10−2 (WAD-4)

158

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

Fig. 8. Calcite dissolution examples.

Fig. 9. Effect of Richardson extrapolation on the calcite dissolution problems.

Cj =

32 E (Δt /4) 21 j

4

− 7 E j (Δt /2) +

1 E (Δt) 21 j

+ O (Δt 4)

(14)

This equation defines an approximation that is one order higher than Eq. (13), and it is obtained sequentially from three successive solutions in the same time interval, which are calculated with different steps Δt , Δt/2 and Δt/4 . Also, each time step has to be further divided into two sequential stages of transportreaction. The result of the previous division procedure is that the calculation on each time interval is an alternating sequence of transport and reaction stages, which renders a solution that is closer to an actual coupled solution where the transport and kinetics take place simultaneously. While Carrayrou et al. (2004) modified the SNIA by using two time steps as it is done in the Strang-splitting method, the Richardson extrapolation can be considered a further improvement to SNIA, because it uses more than two time steps at each iteration. Therefore, in Richardson method there are more chemical correction steps for the concentrations after the transport stage of each iteration. In other coupling algorithms such as the sequential iterative approach (SIA) (Mugler et al., 2004), the concentrations are calculated by solving the transport and the chemical equations

successively and iteratively until convergence. While methods like SNIA have the inherent error associated with explicit algorithms relative to implicit ones, they can be improved by a number of alternating steps or by Richardson extrapolation.

5. Numerical results Several problems were solved using the CV-RBF method coupled with PHREEQC for the chemistry calculations. The examples encompass several cases which include simple chemical diffusion, chemical equilibrium, cation exchange and kinetics. The results also show the effect of Richardson extrapolation for time stepping on the coupling of the transport and chemical calculations. 5.1. Transport and chemical equilibrium with cation exchange The first problem to compare the proposed numerical solution is the advective transport with cation exchange (Appelo and Postma, 1993), where the chemical composition of the effluent from a column containing a cation exchanger is simulated. Initially, the column contains a sodium–potassium–nitrate solution in equilibrium with a cation exchanger. The column is then flushed

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

159

Fig. 10. RMS difference as a function of Richardson extrapolation order for an equilibrium problem.

Fig. 11. RMS difference as a function of Richardson extrapolation order for a chemical kinetics case (CAL-2).

with three pore volumes of calcium chloride solution, and calcium, potassium and sodium react up to equilibrium with the exchanger at all times. The number of pore volumes is equivalent to a characteristic time and it can be considered as the time that the flow takes to go from the top of the column to the bottom i.e. PV ¼L/u where PV is the number of pore volumes, L = 0.08 m is the length of the column and u = 2.78 × 10−6 m/s the flow velocity. The initial composition in the column and the influx boundary conditions are presented in Table 2. The number of exchange sites in each cell is 1.1 mmol, and the initial composition of the exchanger is calculated such that it is in equilibrium with the initial composition in the column. The numerical solution obtained by the proposed CV-RBF method is compared with other three results: the solution reported by Appelo and Postma (1993), finite differences (FD) and the reactive multi-component model (PHT3D) by Henning

Prommer (Appelo and Rolle, 2010). The column has 40 cells to be consistent with one of the PHREEQC runs described by Appelo. Fig. 2 shows the concentrations at the exit of the column plotted against the pore volume number. Chloride is a conservative solute and it exits at about one pore volume. The sodium exchanges with the calcium and it is eluted as long as sodium is present in the cation exchanger. Potassium is released after sodium because the log K value in the exchange is larger for the potassium reaction. When all the potassium has been released, the calcium concentration rises to a steady state value that is the same as in the infilling solution. The solutions by CV-RBF, FD and PHT3D are similar, although there is a small difference between these solutions and PHREEQC. This small lag might be caused by the different transport scheme used in PHREEQC. In the FD, CV-RBF and PHT3D methods the transport is accomplished by discretizing the differential equation by using local interpolation schemes, whereas in the PHREEQC

160

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

Fig. 12. Dissolution with simultaneous cation exchange. Case without calcite.

Fig. 13. Dissolution of calcite with simultaneous cation exchange. Case of chemical equilibrium with calcite.

solution the traditional control volume method is used together with a strategy of cell mixing (Parkhurst and Appelo, 2013). Also, Carrayrou et al. (2004) pointed out that additional numerical diffusion is introduced by SNIA schemes if the initial concentrations in the chemical step are those coming from the beginning of the time step. This conclusion could support the time lag present in Fig. 2. The difference with PHREEQC can be diminished by using a higher order Richardson extrapolation as can be seen in Fig. 3. Fig. 3 depicts the effect of Richardson extrapolation of different orders for the time advancing i.e. truncation errors of the order O (Δt 2), O (Δt 3) and O (Δt 4). As it can be seen the concentration values are higher and closer to the benchmark solution values as the order of the method increases.

The convergence of Richardson method as its order increases can be more clearly appreciated in Fig. 4 that shows the RMS deviation between the numerical results and the solution obtained by PHREEQC alone. From this figure it is inferred that the convergence is uniform for each variable. Also, note that the deviation reduction when Richardson order is increased from 2 to 3 is higher than when the order is increased from 3 to 4, which suggests that increasing orders of Richardson can be regarded as series of successive decreasing corrections on the numerical solution. The reduction of the RMS deviation achieved when the order is raised suggests that the alternating steps introduced by Richardson extrapolation effectively improve the operator splitting coupling between the transport equations and the chemistry reaction

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

Table 5 RMS difference between CV-RBF and other benchmark solutions for the WAD problem.

Table 6 Computation time comparison for several cases. Problem

Comparison

Ca

Na

Case without calcite PHREEQC 2.67 × 10−03 FD 5.18 × 10−03 Saaltink et al. 7.20 × 10−02 Equilibrium case PHREEQC 2.50 × 10−03 FD 4.28 × 10−03 Saaltink et al. 6.70 × 10−02

Cl

∑ (Ci,CV − RBF i

6.59 × 10−05

2.88 × 10−04

3.15 × 10−03

1.20 × 10−04

5.62 × 10−04

6.21 × 10−03

7.50 × 10−02

7.15 × 10−02

7.89 × 10−02

5.60 × 10−05

2.65 × 10−04

3.65 × 10−03

1.15 × 10−04

4.30 × 10−04

5.30 × 10−03

7.50 × 10−02

6.20 × 10−02

6.70 × 10−02

− Ci,0 )2 /N

Computation time (s)

Mg

equations. As it was addressed by Nardi et al. (2014), the SNIA does not introduce global convergence problems, although a tight control on the time step is required in order to minimize the operator-splitting errors. This issue is precisely what is tackled with Richardson extrapolation. The convergence analysis of the solution can be seen in Fig. 5 where different meshes of 40, 100 and 200 cells are used. As expected the accuracy improves with higher number of volumes and the pore volume gap i.e. the time gap between the reference solution and the CV-RBF results becomes smaller as the mesh is denser. It is pointed out that when the volumes are smaller, the local RBF interpolation of the face values is supposed to improve and this positively affects the results from the finite volume method. Fig. 6 shows the RMS deviation between the numerical results and the solution obtained by PHREEQC alone. Observe that for mesh densities less than 100 volumes, the deviation increases and so for coarse discretizations it might be recommended to use more collocation nodes in the Hermite interpolation of the CV-RBF. Table 3 presents a comparison of the CV-RBF results with other methods such as PHREEQC, Finite differences and the PHT3D code (Appelo and Rolle, 2010). The RMS difference in the table is defined as

RMS =

161

(15)

Cation Exchange 120 CAL problem equilibrium 1320 CAL Problem kinetics 4380 WAD problem 10,800

where Ci,CV − RBF is the solution obtained by CV-RBF, Ci,0 is the other solution to compare with and N is the number of data. As can be concluded from the table, the CV-RBF result is in agreement with the other 3 solutions presented. 5.2. Dissolution of calcite The second numerical example is the dissolution of calcite by infiltrating water in a one-dimensional domain (see Fig. 7). The initial and boundary conditions are given in Table 4. This case consists of several sub-problems: one assuming dissolution in equilibrium CAL-E and other four cases with chemical kinetics, with varying dissolution rates i.e. problems CAL-1 to CAL-4, the first having the slowest rate and the last the fastest rate. A case without calcite CAL-0 was also solved. The dissolution rate (r) is calculated form the following expression:

r = κσ (1 − Ω)

(16)

where κ is the reaction rate constant, s the specific reactive surface and Ω the saturation index. Fig. 8 displays the solution pH, after one pore volume has been flushed, that is equivalent to the time for a solute to pass through the whole column. The pH rises due to the increase in the rate of dissolution of calcite, because CO2− 3 from dissolved calcite combines with H+ to form HCO−3 . Therefore, pH is higher for the equilibrium than for the case without calcite. The kinetic circumstances lie between these two extremes, hence a case with high rate kinetics approaches equilibrium. Note that Fig. 8 compares three different solutions: CV-RBF, the results presented by Saaltink et al. (2001) and the results obtained by PHREEQC as standalone

Fig. 14. Effect of time step for the problem of dissolution of calcite with simultaneous cation exchange.

162

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

transport-reaction program. In Fig. 8, CAL-0 refers to the case without calcite, CAL-E to the equilibrium case, and CAL-1 to CAL-4 correspond to different kinetics problems with varying reaction rate as specified in Table 4. In Fig. 9 it can be seen the effect of Richardson extrapolation on this problem. Again, by using Richardson method it is possible to reduce the difference between the CV-RBF solution and the benchmark solution, because with Richardson method there are several chemistry calculations within each time step instead of just one chemical calculation as stated in the traditional SNIA approach. Fig. 10 shows the overall RMS difference between CV-RBF and PHREEQC for different orders of Richardson extrapolation in the equilibrium problem. As can be seen there is a considerable reduction of the error between the case of simple SNIA compared to Richardson extrapolation of order 2, and this trend continues as the order is increased. For instance, note that at the inlet the initial error was equal to 0.014 with simple SNIA, while the error with an order 2 Richardson extrapolation was reduced to 0.007 and even further down to 0.003 with an order 4. Similarly, an equivalent error trend can be observed for the kinetic case depicted in Fig. 11. Finally, it is important to mention that for pH calculations it is necessary to satisfy the charge balance (Hiorth et al., 2013). Thus, in the CV-RBF approach, PHREEQC performs the charge balance at each cell and the total H+ and O− concentrations are advected which allows for the accurate calculation of pH. 5.3. Dissolution of calcite with simultaneous cation exchange The last example solves the flushing of saline water by fresh water in lake Waddenzee (the Netherlands), in a one-dimensional domain (Appelo and Postma, 1993, example 5.8). This problem is named WAD for its reference and it includes both cation exchange and dissolution of calcite. Figs. 12 and 13 display the time variation of the breakthrough concentrations for the cases without calcite and equilibrium with calcite respectively. In Fig. 12 different numerical solutions obtained by several methods are being compared. The methods plotted and compared on the figure are the CV-RBF, finite differences FD, the coupled solution presented by Saaltink et al. (2001) and PHREEQC. Although the solutions agree well, Saaltink et al. solution was obtained by a coupled non-linear approach with Newton's method and not by the SNIA. Initially saline water comes out. Next Na is exchanged by Ca and Mg, while which Na concentration stays relatively high in comparison to Mg and Ca concentrations. Then, adsorbed Mg is exchanged by Ca, leading to an increase in Mg concentration in the effluent. Finally, Ca has replaced the rest of the exchanged cations and water comes out with the same composition as the recharging water. The case of calcite in equilibrium shows approximately the same behavior. Only that during the equilibrium reaction the Na species remains high for a shorter time and the solution contains more Na in response to the addition of Ca by the dissolution of calcite. Table 5 shows the RMS difference between CV-RBF with Richardson method of order 4 and other three solutions. The small difference between CV-RBF and the results by Saaltink et al. (2001) might be explained by the fact that the latter authors used their own chemistry calculation code, while in the present work all the chemistry equations were solved by PHREEQC. An indication of this fact may be that the smallest difference in the Table 5 corresponds to the comparison with PHREEQC. It is also worth pointing out that Saaltink et al. used finite elements for solving the transport equations, while the Hermite CV-RBF collocation approach is used in the present work. The effect of time step on the solution of this problem is shown

in Fig. 14. As it should be expected, the results improve as the time step is reduced. However, a smaller time step leads to higher computational cost, because in the operator splitting technique the coupled problem is not exactly solved, but it is approximated by a sequence of transport steps followed by reaction steps. Table 6 shows the computation time for some problems solved in a desktop Intel Triple Core 2 GHz. As already pointed out by Zabka and Sembera (2011), the coupling of reactive chemistry with transport significantly extends the computation time. The above is particularly true for the kinetic cases where there are several Runge Kutta steps for solving the chemistry equation at each time step.

6. Conclusions A control volume Hermite RBF (CV-RBF) scheme has been implemented and tested for groundwater reactive transport problems that include cation exchange, chemical equilibrium and chemical kinetics. With this method the fluxes on the faces of each volume are calculated by using an Hermite RBF interpolation that includes the differential operator in the interpolation functions themselves. Therefore, the up-winding is implicit and it is not necessary to recourse to any artificial numerical upwind schemes that characterize the traditional control volume methods. The solution of the system of transport equations for each chemical species together with the chemical reaction is accomplished by an operator splitting method following a sequential non-iterative procedure SNIA. The standard SNIA approach was improved by a Richardson extrapolation method which produces an interspersed sequence of chemical reaction and transport calculations at each time step. In this work all the chemistry calculations at each control volume were done by the PHREEQC subroutines, so this work also exemplifies for the first time how to couple the well documented PHREEQ modules with a transport solution obtained by a local RBF method. The results presented in the examples suggest that the transport solution by CV-RBF is effective and in agreement with other numerical methods presented in the literature for the same problems. Also, the proposed Richardson extrapolation is an effective strategy to increase the order of the approximation when coupling the transport and the chemistry equations. Further development of the proposed approach is in progress, in particular to include three dimensional cases and the coupled solution of the Darcy equations together with the chemical equations and the diffusion–advection transport model.

References Appelo, C.A.J., Postma, D., 1993. Geochemistry, Groundwater and Pollution. Balkema, A.A., Brookfield, VT. Appelo, C., Rolle, M., 2010. Pht3d: a reactive multicomponent transport model for saturated porous media. Groundwater 48 (5), 627–632. Bustamante, C.A., Power, H., Florez, W.F., 2013. A global meshless collocation particular solution method for solving the two-dimensional Navier–Stokes system of equations. Comput. Math. Appl. 65, 1939–1955. Bustamante, C.A., Power, H., Florez, W.F., 2014. An efficient accurate local method of approximate particular solutions for solving convection–diffusion problems. Eng. Anal. Bound. Elem. 47, 32–37. Carrayrou, J., Mose, R., Behra, P., 2004. Operator-splitting procedures for reactive transport and comparison of mass balance errors. J. Contam. Hydrol. 68 (3–4), 239–268. Charlton, S.R., Parkhurst, D.L., 2011. Modules based on the geochemical model phreeqc for use in scripting and programming languages. Comput. Geosci. 37 (10), 1653–1663. Chen, C., Fan, C., Wen, P., 2012. The method of approximate particular solutions for solving certain partial differential equations. Numer. Methods Partial Differ. Equ. 28 (2), 506–522. Fasshauer, G.E., Surface fitting and multiresolution methods. In: Solving Partial Differential Equations by Collocation with Radial Basis Functions. University

W.F. Florez et al. / Computers & Geosciences 76 (2015) 151–163

Press, Nancy, 1997, pp. 131–138. Florez, W., Bustamante, C., Giraldo, M., Hill, A., 2012. Local mass conservative hermite interpolation for the solution of flow problems by a multi-domain boundary element approach. Appl. Math. Comput. 218 (11), 6446–6457. Fornberg, B., Flyer, N., 2005. Accuracy of radial basis function interpolation and derivative approximations on 1-d infinite grids. Adv. Comput. Math. 23 (1–2), 5–20. Hiorth, A., Jettestuen, E., Cathles, L.M., 2013. Precipitation, dissolution, and ion exchange processes coupled with a lattice Boltzmann advection diffusion solver. Geochim. Cosmochim. Acta 104, 99–110. Jumarhon, B., Amini, S., Chen, K., 2000. The Hermite collocation method using radial basis functions. Eng. Anal. Bound. Elem. 24 (7–8), 607–611. Kanney, J.F., Miller, C.T., Kelley, C., 2003. Convergence of iterative split-operator approaches for approximating nonlinear reactive transport problems. Adv. Water Resour. 26 (3), 247–261. Kansa, E.J., Hon, Y.C., 2000. Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations. Comput. Math. Appl. 39, 123–137. Kolditz, O., Bauer, S., Bilke, L., 2012. Opengeosys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (thm/c) processes in porous media. Environ. Earth Sci. 67 (2), 589–599. Liu, F., Turner, I., Anh, V., 2002. An unstructured mesh finite volume method for modelling saltwater intrusion into coastal aquifers. Korean J. Comput. Appl. Math. 9 (2), 391–407. Madych, W.R., Nelson, S.A., 1990. Multivariate interpolation and conditionally positive definite functions. Math. Comput. 54, 211–230. Micchelli, C.A., 1986. Interpolation of scattered data: distance matrices and conditionally positive definite functions. Constr. Approx. 2(1), 11–22. Moroney, T., Turner, I., 2006. A finite volume method based on radial basis functions for two-dimensional nonlinear diffusion equations. Appl. Math. Model. 30 (10), 1118–1133. Mugler, C., Montarnal, P., Dimier, A., 2004. Reactive transport modelling on the alliances software platform. In: Miller, C., Farthing, M., Gray, W. (Eds.), 15th International Conference on Computational Methods in Water Resources, vol. 55, pp. 1103–1115. Nardi, A., Idiart, A., Trinchero, P., de Vries, L., Molinero, J., 2014. Interface comsolphreeqc (icp), an efficient numerical framework for the solution of coupled multiphysics and geochemistry. Comput. Geosci. 69, 10–21. Nasir, O., Fall, M., Evgin, E., 2014. A simulator for modeling of porosity and permeability changes in near field sedimentary host rocks for nuclear waste under climate change influences. Tunn. Undergr. Space Technol. 42, 122–135. Nguyen-Van, H., Mai Duy, N., Tran-Cong, T., 2007. A simple and accurate four-node quadrilateral element using stabilized nodal integration for laminated plates. Comput. Mater. Contin. 6, 159–176. Ollivier-Gooch, C., Altena, M.V., 2002. A high-order-accurate unstructured mesh finite-volume scheme for the advection–diffusion equation. J. Comput. Phys.

163

181 (2), 729–752. Orsini, P., Power, H., Morvan, H., 2008. Improving volume element method by meshless radial basis function. Comput. Model. Eng. Sci. 23 (3), 187–207. Orsini, P., Morvan, Power H., Lees, M., 2009. An implicit upwinding volume element method based on meshless radial basis function techniques for modelling transport phenomena. Int. J. Numer. Methods Eng. 81 (1), 1–27. Parkhurst, D., Appelo, C.A., 2013. PHREEQC Version 3–A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations, 〈http://wwwbrr.cr.usgs.gov/projects/GWC_coupled/ phreeqc/phreeqc3-html/phreeqc3.htm〉. Pearson, J., 2013. A radial basis function method for solving pde-constrained optimization problems. Numer. Algorithms 64 (3), 481–506. Powell, M.J.D., 1994. The uniform convergence of thin plate spline interpolation in two dimensions. Numer. Math. 68, 107–128. Saaltink, M.W., Carrera, J., Ayora, C., 2001. On the behavior of approaches to simulate reactive transport. J. Contam. Hydrol. 48 (3–4), 213–235. Samper, J., Xu, T., Yang, C., 2009. A sequential partly iterative approach for multicomponent reactive transport with core2d. Comput. Geosci. 13 (3), 301–316. Schaback, R., Multivariate interpolation and approximation by translates of a basis function, in: Approximation Theory VIII, World Scientific Publishing Co, Singapore, 1995, pp. 1–35. Sladek, V., Sladek, J., Tanaka, M., 2005. Local integral equations and two meshless polynomial interpolations with application to potential problems in nonhomogeneous media. Comput. Model. Eng. Sci. 7 (1), 69–84. Stevens, D., Power, H., Lees, M., Morvan, H., 2009. The use of PDE centres in the local RBF Hermitian method for 3D convective–diffusion problems. J. Comput. Phys. 228 (12), 4606–4624. Versteeg, H.K., Malalasekera, W., 2007. An Introduction to Computational Fluid Dynamics: The Finite Volume Method, 2nd Ed. Pearson Education Limited, Harlow, Essex, U.K. Wissmeier, L., Barry, D.A., 2011. Simulation tool for variably saturated flow with comprehensive geochemical reactions in two- and three-dimensional domains. Environ. Model. Softw. 26 (2), 210–218. Wright, G., Fornberg, B., 2006. Scattered node compact finite difference-type formulas generated from radial basis functions. J. Comput. Phys. 212 (1), 99–123. Xu, T., Samper, J., Ayora, C., Manzano, M., Custodio, E., 1999. Modeling of non-isothermal multi-component reactive transport in field scale porous media flow systems. J. Hydrol. 214 (1–4), 144–164. Zabka, V., Sembera, J., 2011. Coupled reactive transport modeling - the program transport. In: Papadrakakis, M., Onate, E., Schrefler, B. (Eds.), 4th International Conference on Computational Methods for Coupled Problems in Science and Engineering (COUPLED PROBLEMS), pp. 1400–1407. Zongmin, W., 1992. Hermite–Birkhoff interpolation of scattered data by radial basis functions. Approx. Theory Appl. 8 (2), 1–10.