The conservative splitting domain decomposition method for multicomponent contamination flows in porous media

Journal of Computational Physics 400 (2020) 108974

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

The conservative splitting domain decomposition method for multicomponent contamination flows in porous media Dong Liang a,∗ , Zhongguo Zhou b,∗ a b

Department of Mathematics and Statistics, York University, Toronto, Ontario, M3J 1P3, Canada School of Information Science and Engineering, Shandong Agricultural University, Taian, Shandong, 271018, China

a r t i c l e

i n f o

Article history: Received 30 January 2018 Received in revised form 12 August 2019 Accepted 19 September 2019 Available online 24 September 2019 Keywords: Mass conservative S-DDM Solution-flux technique Modified upwind Non-overlapping domain decomposition Nonlinear multicomponent contamination flow

a b s t r a c t In the paper, a new conservative splitting decomposition method (S-DDM) is developed for computing nonlinear multicomponent contamination flows in porous media over multiblock sub-domains. On each block-divided sub-domain, we take three steps to solve the coupled nonlinear system of water-head equation and multicomponent concentration equations in each time interval. The interface Darcy’s velocity and the interface global concentration fluxes are first predicted by the semi-implicit flux schemes, while the solutions of water-head and multicomponent concentrations, Darcy’s velocity and global concentration fluxes in the interiors of sub-domains are computed by one-directional splitting implicit solution-flux coupled schemes on staggered meshes, and finally the interface Darcy velocity and global concentration fluxes are corrected by the interior solutions. The significance of our scheme is that while it keeps the advantages of the non-overlapping domain decomposition and the splitting technique, it preserves mass on the whole domain of domain decompositions. Numerical experiments are presented to illustrate the excellent performance of our proposed conservative S-DDM approach for computing nonlinear multicomponent contamination flows in groundwater. The developed algorithm of the conservative S-DDM works efficiently over multiple block-divided subdomains, which can be applied in simulation of large scale multicomponent contamination flows in parallel computing. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Recently, groundwater resources have been seriously threatened by the leaching contaminants from industrial wastes, agricultural fertilizers and domestic sewages, etc. There were effects on characterizing the multicomponent transport and nonlinear reaction processes in porous media. In general, the multicomponent reactions can be described as kinetically controlled dissolution-precipitation reactions (see [2,14,16,33]) or the geochemical equilibrium reactions as hydrolysis aqueous complexation, oxidation-reduction, ion exchange, surface complexation, and gas dissolution-exsolution reaction (see [18,40, 41]). Mathematical models describing multicomponent contamination flows in porous media are the time-dependent and coupled nonlinear partial differential equations ([1,10,11,17,32,42]). Due to the complexity, the nonlinearity, and the large scale and long term prediction, it is important to develop efficient domain decomposition methods for solving multicomponent contamination flows in parallel computing.

*

Corresponding authors. E-mail addresses: [email protected] (D. Liang), [email protected] (Z. Zhou).

https://doi.org/10.1016/j.jcp.2019.108974 0021-9991/© 2019 Elsevier Inc. All rights reserved.

2

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

Many domain decomposition methods have been proposed for linear problems based on the overlapping domain decompositions (see [26,27,31,37–39], etc.), where the implicit schemes are used on overlapping subdomains and the iterations among subdomains are required such as the Schwarz alternating method. Since non-overlapping domain decomposition methods have the advantages of the noniterativity, the low computation and communication costs at each time step, the non-overlapping domain decomposition methods have attractively been developed in solving time dependent linear equations (see [3,6–9,23,35,44,47,48], etc.). Paper [6] proposed a mixed/hybrid scheme, where the interface solutions were first computed by the explicit schemes on interfaces while the interior solutions were then solved by the implicit schemes in sub-domains, but the scheme is conditionally stable. For improving stability condition, paper [9] proposed the explicitimplicit domain decomposition (EIDD) methods for solving linear parabolic problems, where a multi-step explicit scheme or a high-order explicit scheme was used on the interfaces. By introducing the implicit correction step on the interfaces, a class of corrected explicit-implicit domain decomposition (CEIDD) methods were further studied in papers [35,48], where the zigzag-shaped subdomains were analyzed theoretically the unconditional stability for solving two dimensional linear heat equations in [35]. But, the schemes in [6,9,35,48] work well only for the stripe-divided subdomains along one variable or for the zigzag-shaped subdomains. Papers [34,47] proposed the time three levels schemes for solving linear equations, where the inner boundary solutions were firstly computed by the values of previous two time levels at the interface points, then the values in sub-domains were computed by the fully implicit schemes, and the inner boundary solutions were updated by the fully implicit schemes. Over multi-block subdomains, by combining the non-overlapping domain decomposition and the splitting technique, papers [8,23] developed an efficient explicit-implicit splitting domain decomposition method (S-DDM) for solving the parabolic equations and the compressible contamination fluid flows in porous media, where an efficient local multilevel scheme was used for computing the values on the interfaces of sub-domains, and the interior values were computed by one splitting implicit scheme. However, these previous explicit-implicit domain decomposition methods [7–9,23,34,35,47,48] do not satisfy the physical law of mass conservation over the whole domain. To obtain conservative domain decomposition schemes, there were a few works developed for solving time-dependent linear problems (see [5,46,45]). Paper [5] presented an explicit-implicit conservative domain decomposition procedure for linear parabolic equations but over only stripe-divided domain decompositions along one variable direction in two dimensions, where the fluxes at the subdomain interfaces were calculated by an average operator from the solutions at the previous time level. The conditional stability and the error estimate were carried out over the stripe-divided subdomains in two dimensions. Paper [46] studied the cell centered finite difference domain decomposition procedure for linear heat equations with constant coefficients in one dimension, but it can not be applied to solve high dimensional equations over multi-block domain decompositions. The stripe-divided domain decomposition methods only work well for application problems with a thin domain, where the whole domain is required to be divided only along one variable direction and the size of the whole domain is very small in the other direction. By combining the operator splitting technique with the solution-flux staggered mesh, paper [45] studied the mass-preserving splitting domain decomposition method for solving linear parabolic equations over multiple block-divided domain, where the interface fluxes were first computed by the semiimplicit (explicit) flux scheme while the interior values were computed by the one-dimensional splitting implicit scheme, and further the interface fluxes were corrected on interfaces. However, to our knowledge, there is no conservative domain decomposition methods developed for nonlinear multicomponent contamination flows in porous media. The mass conservation is an important physical law in multicomponent contamination flows in porous media. Preserving mass of numerical domain decomposition schemes over the whole domain is important and required for solving the time dependent and coupled nonlinear partial differential equations in parallel computing, especially for long term simulations and large scale problems. Therefore, it is a challenging and practically significant task to study the mass conservative domain decomposition methods for solving nonlinear multicomponent contamination flows in porous media. In this paper, we propose a new mass conservative splitting domain decomposition method over multi-block subdomains for nonlinear multicomponent contamination flows in porous media, where the non-overlapping block-divided domain decomposition, operator splitting, modified upwind, solution-flux technique on staggered mesh are used to solve the coupled nonlinear system of water-head equation and multicomponent concentration equations over multiple subdomains. On each block-divided sub-domain, our scheme proposes that firstly, the interface x-directional Darcy’s velocity and the interface x-directional global concentration fluxes are completed by the semi-implicit flux schemes, and the intermediate solutions and fluxes of water-head and multicomponent concentrations in the interiors of sub-domains are then computed by x-dimensional splitting implicit solution-flux coupled schemes, where the second-order modified upwind technique is proposed to approximate the multicomponent contamination concentration equations in the interiors of subdomains, and the interface Darcy’s velocity and the global concentration fluxes are corrected by the interior intermediate solutions; Secondly, the interface y-directional Darcy’s velocity and the interface y-directional global concentration fluxes are completed by the flux schemes, and the solutions and fluxes of water-head and multicomponent concentrations in the interiors of subdomains are then computed by y-dimensional splitting implicit solution-flux coupled schemes, and the interface Darcy’s velocity and global concentration fluxes are finally corrected by the interior solutions. The important feature of our scheme is that it conserves mass over the whole domain of domain decompositions for nonlinear multicomponent contamination flows, while the modified upwind technique ensures the second order accuracy to the multicomponent concentration equations. Our approach not only keeps the advantages of the non-overlapping domain decompositions and the splitting technique, while it satisfies the physical law of mass conservation, where we prove theoretically it conserves mass over the whole domain of multiple subdomains. It reduces computational complexities, large memory requirements and long compu-

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

3

tational durations, and is more suitable and powerful for parallel computing. Numerical experiments are first presented to illustrate the mass conservation, convergence order, stability and parallelism. We then simulate the multicomponent contaminant transport accompanied with precipitation reactions between contaminants Ca2+ , CO23− , Cl− and Na+ in groundwater. When the contaminant HCl flows into the porous media which includes CaCO3 and MgCO3 , there will produce dissolution reaction and give the new contaminants Ca2+ , Mg2+ and H2 CO3 . Numerical results show that our proposed mass-preserving S-DDM scheme simulates well multicomponent contamination flows in porous media. The developed algorithm of the conservative S-DDM works efficiently over multiple block-divided sub-domains, which can be applied in simulation of large scale multicomponent contamination flows in parallel computing. The rest of this paper is organized as follows. In Section 2, we describe mathematical models of the multicomponent contamination flows in porous media. In Section 3, we propose the conservative S-DDM scheme for solving the coupled nonlinear system of water head equation and multicomponent concentration equations. We prove mass conservation of sour scheme in Section 4 and a linearization iterative algorithm is proposed in Section 5. Numerical experiments are performed to illustrate the efficiency of our approach in Section 6. 2. Mathematical model of nonlinear contamination flow in porous media The mathematical model of multicomponent contamination flows in porous media can be described as the timedependent coupled nonlinear partial differential equations by two sets of governing equations of pressure equation and multicomponent contamination concentration equations. 2.1. Flow equation Based on the fluid mass conservation, the groundwater flow equation is

∂ (φ ρ ) + ∇ · (ρ u) = Q ρ , ∂t

(1)

where φ is the volumetric porosity, ρ is the fluid density, and Q is the flow volume rate of sources (sinks). Darcy’s law establishes the relationship between the velocity and the pressure gradient and can be described as

u=−

k

μ

(∇ p + ρ g ∇ z),

(2)

where k is the intrinsic permeability matrix, μ is the fluid dynamic viscosity, g is the magnitude of gravitational accelera∂(φ ρ ) tion, and z is the depth of a reference point to the datum. Let S p be the special storage coefficient, defined by S p = ∂ p , if we assume that the liquid and the porous media may both be compressible weakly, we can get that ∂p

S p ∂ t + ∇ · (ρ u) = Q ρ .

(3)

2.2. Multicomponent concentration equations The aqueous contaminants are transported by advection and diffusion in groundwater. Simultaneously, each contaminant is influenced by sources (or sinks) and geochemical reactions. The equations of multicomponent concentrations c α , α = 1, 2, · · · , N c , are described as ∂ ∂ t (φ c α ) + ∇

· (uc α − D∇ c α ) = Q c α∗ + R α (c ),

α = 1, 2, · · · , N c ,

(4)

∗ is the given special contaminant concentration at the inject well while where D is the diffusion coefficient matrix, and c α c ∗ = c at the production well.

α

α

The kinetically controlled reaction R α (c ), where  c = (c 1 , c 2 , · · · , c N c ), is described as the intra-aqueous reactions R aα (c ) s s and mineral dissolution-precipitation reactions R α (c ), i.e. R α (c ) = R aα (c ) + R α (c ) (see [20,24,29]). 2.3. The system of the coupled water-head and multicomponent concentrations In studying multicomponent contamination flows in groundwater, the water head h is used to replace the fluid pressure as

h=

p

ρg

+ z.

(5)

Then, from (2), (3) and (4), and letting the specific storage S s = S p g account for both porous media and fluid compresskρ g

ibility and K = μ be the hydraulic conductivity, we can have the coupled system of multicomponent contamination flow in two dimensions as follows

4

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

∂h − ∇ · (K∇ h) = Q , x ∈ , t ∈ (0, T ], ∂t x ∈ , t ∈ (0, T ], u = −K∇ h, ∂(φ c α ) + ∇ · (uc α − D∇ c α ) = Q c α∗ + R α (c ), α = 1, 2, · · · , N c , ∂t

Ss

(6) (7) x ∈ , t ∈ (0, T ],

(8)

where  = [ax , b x ] × [a z , b z ], u = (u x , u z )τ , K = diag ( K xx , K zz ) and D = diag ( D xx , D zz ). The boundary conditions can be given as Dirichlet boundary, Neumann boundary, or mixed boundary. We will consider the Neumann boundary conditions by

 = 0, K∇ h · n

(uc α − D∇ c α ) · n = 0,

α = 1, 2, · · · , N c ,

x ∈ ∂ ,

t ∈ (0, T ],

(9)

and the initial conditions are provided as

h(x, 0) = h0 (x),

c α (x, 0) = c α 0 (x),

x ∈ .

(10)

Besides of the Darcy’s velocity u = (u , u ), we define the concentration global flux w = { w α , w α } as x

x wα = u x c α − D xx

∂ cα , ∂x

z

x

z wα = u z c α − D zz

∂ cα , ∂z

z

α = 1, 2, · · · , N c .

(11)

By (11), Eqs. (6)–(8) are rewritten the following equivalent forms

∂ h ∂ ux ∂ uz + + = Q, x ∈ , t ∈ (0, T ], ∂t ∂x ∂z ∂h ∂h u x = − K xx , u z = − K zz , x ∈ , t ∈ (0, T ], ∂x ∂z ∂ w αz ∂(φ c α ) ∂ w αx + + = Q c α∗ + R α (c ), α = 1, 2, · · · , N c , x ∈ , t ∈ (0, T ], ∂t ∂x ∂z ∂ cα ∂ cα x wα = u x c α − D xx , w αz = u z c α − D zz , x ∈ , t ∈ (0, T ], ∂x ∂z

Ss

(12) (13) (14) (15)

with the boundary conditions

 = 0, u·n

 = 0, wα · n

x ∈ ∂ ,

t ∈ (0, T ],

(16)

and the initial conditions

h(x, 0) = h0 (x),

c α (x, 0) = c α 0 (x),

α = 1, 2, · · · , N c , x ∈ .

(17)

2.4. The general system of nonlinear multicomponent contamination flow For more complex multicomponent flows, the density of the liquid may depend on component concentrations while the porous media may be compressible (see [12,30,36,43], etc.).

.

ρ = ρ0 (1 +

Nc

(18)

α =1 ηα c α ),

where ρ0 is the reference density of the freshwater, and ηα = ερα is a coefficient describing the influence of concentration 0 of

α th component on the fluid density.

The general coupled nonlinear system of water-head equation and multi-component concentration equations is

Ss

∂ h ∂ ux ∂ uz + + = Q + f (h, c ), ∂t ∂x ∂z u x = − K xx

x ∈ , t ∈ (0, T ],

c ∂h  + ηα cα ), ∂z

(19)

N

∂h , ∂x

u z = − K zz (

x ∈ , t ∈ (0, T ],

(20)

α =1

∂ w αz ∂(φ c α ) ∂ w αx + + = Q c α∗ + R α (c ) + g α (h, c ), α = 1, · · · , N c , x ∈ , t ∈ (0, T ] ∂t ∂x ∂z ∂ cα ∂ cα x wα = u x c α − D xx , w αz = u z c α − D zz , x ∈ , t ∈ (0, T ], ∂x ∂z

(21) (22)

with the boundary conditions

 = 0, u·n

 = 0, wα · n

x ∈ ∂ ,

t ∈ (0, T ],

(23)

and the initial conditions

h(x, 0) = h0 (x),

c α (x, 0) = c α 0 (x),

α = 1, 2, · · · , N c , x ∈ .

(24)

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

5

Fig. 1. The domain  and its sub-domains  p ,q .

Remark 1. The coupled system (19)–(24) is derived based on the law of mass conservation in fluid flow in porous media. It is clear that the water-head and contamination component concentrations satisfy mass conservations as



t 

 S s h(x, t ) dx =







 φ c α (x, t ) dx =



[ Q + f (h, c )] dxdt ,

S s h0 (x) dx +

(25)

0 

t  φ c α 0 (x) dx +

t 

∗

0 



[ R α (c ) + g α (h, c )] dxdt ,

Q c α dxdt + 0 

α = 1, 2, · · · , N c .

(26)

In the next section, we will construct the mass conservative domain decomposition scheme for the coupled nonlinear system over non-overlapping domain decomposition.

3. The mass conservative S-DDM scheme for multicomponent contamination flows In this section, we will give the non-overlapping block-divided domain decomposition and its staggered meshes and propose the mass conservative splitting domain decomposition method (S-DDM) for solving the solution and its fluxes to the coupled nonlinear system of water head equation and multicomponent concentration equations. 3.1. Domain decomposition and staggered mesh The large domain  is divided into multiple non-overlapping block-divided sub-domains  p ,q , p = 1, 2, · · · , N p , q = 1, 2, · · · , N q , and let M = N p × N q be the number of sub-domains. Let γ px,q = ∂  p ,q ∩ ∂  p +1,q be the interior interface of  p ,q and  p +1,q , and γ pz,q = ∂  p ,q ∩ ∂  p ,q+1 be the interior interface of  p ,q and  p ,q+1 . The interface set  = x ∪ z consists of the x-directional interface x =

qz

=

N p

N p

x p =1  p

and xp =

Nq

q=1

γ px,q , and the z-directional interface z =

Nq

z q=1 q

with

z p ,q

(see, Fig. 1). Each sub-domain  p ,q is further divided into the staggered mesh (see, Fig. 2). p =1 γ For simplicity of description, we consider the staggered meshes over 2 × 2 block sub-domains (see, Fig. 2). Let the middle points of edges be (xi + 1 , z j ) and (xi , z j + 1 ) and let the center points be (xi , z j ) with 2

xi + 1 = i x, 2

z j + 1 = j z, 2

i = 0, 1 , · · · , I , j = 0, 1, · · · , J ,

2

xi = ( i −

1 2

zj = (j −

) x, 1 2

) z,

i = 1, 2, · · · , I , j = 1, 2, · · · , J ,

where x and z are the step sizes along x-direction and z-direction respectively, and x = b−I a , z = d−J c , and I , J

are integer numbers. The lines x : x = xi + 1 and z : z = z j + 1 are the interior interface lines of sub-domains. Over the 1 2 1 2 staggered mesh, we denote the solution on node (xi , z j ) and the x-direction and z-direction velocities (or fluxes) on nodes (xi + 1 , z j ) and (xi , z j + 1 ). 2

2

6

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

Fig. 2. The staggered meshes in sub-domains: ◦, , – the points of (i , j ), (i +

Let

1 , 2

j ), (i , j +

1 ). 2

ωi, j = ω(xi , z j ), ωi+ 1 , j = ω(xi+ 1 , z j ) and ωi, j+ 1 = ω(xi , z j+ 1 ). Define the spatial difference operators as 2

δ x ωi + 1 , j =

ω i +1 , j − ω i , j x

2

2

, δ x ωi , j =

2

2

ωi + 1 , j − ωi − 1 , j 2

2

x

, δ z ωi , j + 1 =

ω i , j +1 − ω i , j

2

z

, δ z ωi , j =

ωi , j + 1 − ωi , j − 1 2

2

z

.

3.2. The time step operator splitting to the nonlinear multicomponent contamination flow system For the general coupled nonlinear system of water-head equation and multi-component concentration equations  = ( S s h, φ c1 , φ c2 , · · · , (19)–(24), we will consider the time step operator splitting process. Let unknown variable vector y φ c N c )τ . Then, the system (19)–(24) can be rewritten as

∂ ∂ ∂ y = − ψ x − ψ z + F , ∂t ∂x ∂z

(27)

 x = (u x , w x , w x , · · · , w x )τ and ψ z = (u z , w z , w z , · · · , w z )τ are where the components of flux vectors ψ 1 1 Nc Nc 2 2 u x = − K xx

∂h , ∂x

c ∂h  + ηα c α ) ∂z

N

u z = − K zz (

(28)

α =1

x = u x c α − D xx wα

∂ cα , ∂x

z wα = u z c α − D zz

∂ cα , ∂z

α = 1, 2, · · · , N c

(29)

and the nonlinear source vector is

F = ( Q + f (h, c ), Q c 1∗ + R 1 ( c ) + g 1 (h, c ), · · · , Q c ∗N c + R N c ( c ) + g N c (h, c ))τ .

(30)

Let t be the time step size, and N = T / t be a total step number, denote t n = n t. In the time interval (t n , t n+1 ], a time step operator splitting ([19]) is used for solving the system (27)–(28) along x-direction and z-direction respectively. Let A1 represent the x-direction operator − ∂∂x and A2 represent the z-direction operator − ∂∂z + F , including the nonlinear source operator F . The time step operator splitting is as





(t + t ) = e t A1 e t A2 y(t ) + O ( t ). y

(31)

The time step splitting algorithm of the system of nonlinear multicomponent flow is described as that, for n = 0, 1, · · · , N − 1, do 1 n , it solves Stage 1: The first half time step [t n , t n+ 2 ] along x-direction, with y

∂ ∂ y = − ψ x ∂t ∂x

(32)

 x = (u x , w x , w x , · · · , w x )τ by and ψ 1 Nc 2 u x = − K xx

∂h , ∂x

x wα = u x c α − D xx

∂ cα , ∂x

α = 1, 2, · · · , N c

(33)

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974 1

7

1

n+ 2 , it solves Stage 2: The first half time step [t n+ 2 , t n+1 ] along z-direction, with y

∂ ∂ y = − ψ z + F , ∂t ∂z

(34)

 z = (u z , w z , w z , · · · , w z )τ by and ψ 1 Nc 2 c ∂h  + ηα cα ), ∂z

N

u z = − K zz (

z wα = u z c α − D zz

α =1

∂ cα , ∂z

α = 1, 2, · · · , N c

(35)

We will further discretize Stage 1 (32) (33) and Stage 2 (34) (35) on space domain decomposition and staggered mesh (see Fig. 1 and Fig. 2) defined in Section 3.1. The difficulty of construction is to ensure mass conservation of the numerical schemes to Stages 1 and 2 over the whole domain of domain decomposition, since many existing domain decomposition schemes break mass conservation over domain decompositions. In the next subsection, we will propose the mass preserving solution-flux coupled schemes to Stages 1 and 2 over domain decomposition of staggered meshes. 3.3. The mass conservative S-DDM scheme for the nonlinear coupled water-head and multicomponent concentration equations We take two half-time steps to compute the coupled nonlinear system of water-head equation and multicomponent con1

centration equations. In the first half time interval [t n , t n+ 2 ], we solve the intermediate water-head and Darcy’s velocity, and the intermediate concentrations and its global concentration flux by using the splitting equation along x-direction. The predict intermediate interface x-direction Darcy’s velocity and intermediate interface global concentration fluxes on interfaces of sub-domains are computed by the interface semi-implicit flux schemes, while we solve the intermediate water-head and the x-directional Darcy’s velocity and the intermediate multicomponent concentrations and the x-directional global concentration fluxes in the interiors of sub-domains by the x-direction splitting implicit solution-flux coupled schemes, and then we correct the interface Darcy’s velocity and the global concentration fluxes by the interior intermediate solutions. 1

In the second half time interval [t n+ 2 , t n+1 ], along z-direction, we find the predict intermediate z-direction Darcy’s velocity and the predict global concentration fluxes on the interfaces of sub-domains, and then compute the water-head and the z-direction Darcy’s velocity, multicomponent concentrations and the z-directional global concentration fluxes by the z-directional splitting implicit solution-flux coupled scheme, and finally, we correct the interface Darcy’s velocity and global concentration fluxes on interfaces of sub-domains by the interior solutions. There are two important and difficult problems that we need to solve in the construction of the mass conservative S-DDM scheme. Firstly, at Step 1, we require the predict intermediate interface x-direction flux (i.e. the predict intermediate interface x-direction Darcy’s velocity) and the intermediate concentration x-direction global fluxes on the interfaces of sub-domains. For doing this, we consider the x-direction flux equations by taking x-direction derivative to equation (32). Then, the semiimplicit flux scheme can be proposed by applying the spacial difference operators to the x-direction flux equations. Letting

{U˜

x,n+ 12 i p + 12 , j

} be the predict intermediate interface x-direction velocity and { W˜ α

x,n+ 12 i p + 12 , j

} be the predict intermediate x-direction

interface global concentration flux, the detailed interface semi-implicit flux schemes are defined as (39) (40) below. Similarly, z,n+1 at Step 2, the predict intermediate interface z-direction flux {U˜ 1 } and the predict intermediate interface z-direction i , jq + 2

z,n+1 global concentration flux { W˜ α i , j + 1 } on the interfaces of sub-domains are approximated by using the interface semi-implicit q 2

flux schemes defined as (47) (48). Secondly, in order to solve effectively the convention dominated multicomponent concentration equations, we use the modified upwind technique to approximate the convection and diffusion terms (see, for examples, [21,25]). Let λ(s) be the sign function given by



s ≥ 0, s < 0,

1, 0,

λ(s) =

(36)

and apply the second order modified upwind approximations to the global concentration fluxes

W αx i + 1 , j = U x 2

W αz i , j + 1 = 2

[λ(U x 1 )C α i , j i + 12 , j i+ 2 , j U z 1 [λ(U z 1 )C α i , j i, j+ 2 i, j+ 2

+ (1 − λ(U x

))C α i +1, j ] − i + 12 , j z + (1 − λ(U ))C α i , j +1 ] − i , j + 12

ˆ D xx,i + 1 , j δx C α i + 1 , j , 2

2

ˆ D zz,i , j + 1 δ z C α i , j + 1 , 2

(37)

2

where

ˆ D

xx,i + 12 , j

= D xx,i + 1 , j /(1 + 2

|U x

|

i+ 1 , j 2 2D xx,i + 1 , j 2

x),

ˆ D

zz,i , j + 12

= D zz,i , j + 1 /(1 + 2

|U z

i, j+ 1 2

2D

|

zz,i , j + 1 2

z).

(38)

8

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

Further, let {U

x,n+ 12 i + 12 , j

} and {U

head-water, and let { W α

x,n+ 12

z,n+ 12 i , j + 12

n+ 1

} be the intermediate numerical velocities and { H i , j 2 } be the intermediate numerical

} and { W α

z,n+ 12

n+ 1

} be the intermediate numerical global concentration fluxes and let {C α i , j 2 } be the intermediate numerical concentration. Similarly, other notations are defined for the n + 1 time level values. i + 12 , j

i , j + 12

Based on these notations, we can propose the mass conservative splitting domain decomposition method (S-DDM) as follows. For n = 0, 1, · · · , N − 1, do 1 Step 1: The first half time step [t n , t n+ 2 ] along x-direction, (1-a) The predict intermediate interface velocity U˜

˜ W

x,n+ 12 α i p + 12 , j

U˜

on the interface lines

x,n+ 12

φ

n+ 12

W˜ α

γ px,q are computed by the interface semi-implicit flux schemes as below,

+

i p + 12 , j

i p + 12 , j

and the predict intermediate interface global concentration flux

n

t K xx i p + 12 , j x,n x,n+ 12 U − 2U˜ + U x,n 3 , 1 1 2 S ip− 2 , j ip+ 2 , j ip+ 2 , j x sip+ 1 , j 2 x,n x,n x,n t = φ n 1 W α x,n 1 + U i , j λ(U i , j ) W α + (U ixp,n, j (1 − λ(U ixp,n, j )) x p p ip+ 2 , j ip+ 2 , j i p − 12 , j

x,n+ 1 −U ixp,n+1, j λ(U ixp,n+1, j )) W˜ α 12 − U ixp,n+1, j (1 − λ(U ixp,n+1, j )) W α x,n 3 ip+ 2 , j ip+ 2 , j

1 x , n + x,n + xt2 Dˆ n Wα − 2 W˜ α 12 + W α x,n 3 . 1 1

= U x,n

i p + 12 , j

x,n+ 12

i p + 12 , j

x,n+ 12 i p + 12 , j

xx,i p + 2 , j

ip− 2 , j

n+ 12

ip+ 2 , j

x,n+ 12 i + 12 , j

n+ 12

(39)

(40)

ip+ 2 , j

x,n+ 12 α i + 12 , j

(1-b) The intermediate variables H i , j , U and C α i , j , W in the interiors of sub-domains  p ,q are computed by the following x-directional splitting implicit solution-flux coupled scheme,

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

n+ 1 2

S s i, j U U

H i, j

x,n+ 12 i + 12 , j x,n+ 12 i p + 12 , j

n+ 1 φi , j 2

− H ni, j t

= − K xx = U˜

n+ 1 C α i, j 2

x,n+ 12

+ δx U i , j

n+ 1

i + 12 , j

i p + 12 , j

(xi , z j ) ∈  p ,q ,

n+ 1 2

H i +12, j − H i , j

n+ 12

x,n+ 12

= 0,

(41)

(xi + 1 , z j ) ∈  p ,q ,

,

x

2

(xi p + 1 , z j ) ∈ γ px,q ,

,

2

−φin, j C α ni, j

x,n+ 1

+ δx W α i , j 2 = 0, α = 1, 2, · · · , N c , (xi , z j ) ∈  p ,q ,

x,n+ 1 x,n+ 1 n+ 1 x,n+ 1 n+ 1 = U 1 2 λ(U 1 2 )C α i , j 2 + (1 − λ(U 1 2 ))C α i +12, j

t

Wα

x,n+ 12 i + 12 , j

i+ 2 , j

− Dˆ Wα

x,n+ 12 i p + 12 , j

= W˜ α

i+ 2 , j

i+ 2 , j

n+ 1 n+ 1 C α i +12, j −C α i , j 2

n+ 12 xx,i + 12 , j

x,n+ 12 i p + 12 , j

,

x

,

(42)

(xi + 1 , z j ) ∈  p ,q , 2

(xi p + 1 , z j ) ∈ γ px,q . 2

(1-c) The correct intermediate x-direction interface velocity U puted by the interior solutions n+ 1

x,n+ 12 i p + 12 , j

and W α

x,n+ 12 i p + 12 , j

on the interface lines

γ px,q are recom-

n+ 1

H i +21, j − H i , j2 p p U = − K xx 1 , 1 ip+ 2 , j ip+ 2 , j x x,n+ 12

n+ 12

(43) n+ 1

n+ 1



C α i +21, j − C α i , j2 x,n+ 12 n+ 12 x,n+ 12 n+ 12 n+ 12 p p ˆ λ( − = U U ) C + ( 1 − λ( U )) C D , Wα α ip, j α i p +1 , j i p + 12 , j i p + 12 , j i p + 12 , j i + 12 , j xx,i p + 12 , j x x,n+ 12

x,n+ 12

and the intermediate z-direction velocity U

U

z,n+ 12 i , j + 12

= − K zz

n+ 12 i , j + 12

n+ 1

z,n+ 12 i , j + 12

and W α

z,n+ 12 i , j + 12

(44)

on sub-domains  p ,q by

n+ 12

H i , j +21 − H i , j

z

,

(xi , z j + 1 ) ∈  p ,q , 2

(45)

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

Wα

z,n+ 12 i , j + 12

=U

z,n+ 12 i , j + 12

(λ(U

z,n+ 12

n+ 1 )C α i , j 2

i , j + 12

+ (1 − λ(U

z,n+ 12 i , j + 12

n+ 1 ))C α i , j +21 ) −

ˆ D

9

n+ 1

n+ 1

C α i , j2+1 − C α i , j2

n+ 12 zz,i , j + 12

z

(46)

.

1

Step 2: The second half time step [t n+ 2 , t n+1 ] along z-direction, z,n+1 (2-a) The predict z-direction interface Darcy’s velocities U˜ 1 and the predict z-direction interface global concentration

˜ fluxes W

U˜

z,n+1 α i , jq + 12

z,n+1 i , j q + 12

i , jq + 2

on the interface lines

=U

z,n+ 12 i , j q + 12

− t

n+ 1 K zz 2 1 i , jq + 2 Ss i , jq + 1 2

x,n+1

φ n +1

i , j q + 12

W˜ α i , j + 1 q 2

n+ 1 2 i , jq + 1 2 Ss i , jq + 1 2

K zz

t z2

+

=φ

γ pz,q are computed by the interface semi-implicit flux schemes as U

δ z ( Q n +1

i , j q + 12

n+ 12

W

i , j q + 12

z,n+ 12

+ f (H

x,n+ 12 α i, j + 1 q 2

z,n+ 12

1

i , j q − 12

+

z,n+ 2 − 2U˜ z,n+11 + U 3 i , jq + 2

n+ 12 i , j q + 12



, C

i , jq + 2

n+ 12

z,n+ 12

t z U i , jq

λ(U i , jq

z,n+ 12

(47)

)),

i , j q + 12

z,n+ 12



)W α

x,n+ 12 i , j q − 12

z,n+ 12

z,n+ 12

z,n+ 12

+ (U i , jq

(1 − λ(U i , jq 1

z,n+

z,n+ 12

))

z,n+1 2 −U i , jq +1 λ(U i , jq +1 )) W˜ α i , jq + 1 − U i , jq +1 (1 − λ(U i , jq +1 )) W α i , j q + 32 2

n+ 1 z,n+ 12 x,n+ 12 ˜ x,n+1 + zt2 Dˆ 2 1 W α 1 − 2W α i p + 1 , j + W α 3 zz,i , j q + 2

− t Dˆ

n+ 12 zz,i , j q + 12

i , jq − 2

δz

α

{ H ni ,+j 1 , U z,n+11 } i, j+

(2-b) The (n + 1)th level variables

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨

n+ 1 2

S s i, j

t

z,n+1 i , j + 12

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨

U

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

Wα

U

H ni ,+j 1 − H i , j

z,n+1 i , j q + 12

H n +1

2

= U˜ z,n+11 , i , jq + 2

− H ni ,+j 1 z

i , j +1

+

2

are computed by the z-directional splitting

2

(xi , z j ) ∈  p ,q ,

n +1

α = 1 ηα C α i , j + 1 , 2

(49)

(xi , z j + 1 ) ∈  p ,q , 2

(xi , z jq + 1 ) ∈ γ pz,q , 2

n+ 1 C α i, j 2

+ δz W α iz,,nj +1 = Q n+1 C α∗ ni ,+j 1 + R α (C in,+j 1 ) + g α ( H ni ,+j 1 , C in,+j 1 ),

t

α = 1, 2, · · · , N c , (xi , z j ) ∈  p,q ,

z,n+1 i , j + 12



= U z,n+11 λ(U z,n+11 )C α ni ,+j 1 + (1 − λ(U z,n+11 ))C α ni ,+j +1 1 i, j+ 2

i, j+ 2

zz,i , j + 12

z,n+1 Wα j q + 12

=

z

˜ z,n+1 1 , W α i , jq + 2

z,n+1 i , j q + 12

= − K zz n+1

,

(xi , z j + 1 ) ∈  p ,q 2

( xi , z j q + 1 ) ∈ γ 2

H ni ,+j 1+1 − H ni ,+j 1 q

q

i , j q + 12

z,n+1 i , j q + 12

= U z n +1

i , j q + 12 , j

z (λ(U z n+1

(50)

i, j+ 2

1 C α ni ,+j + −C α ni ,+j 1 1

− Dˆ n+1

the interior solutions

Wα

1 1 ) + g α ( H n+ 2 , C n+ 2 ))i , jq + 1 .

{C α ni ,+j 1 , W α z,n+11 } i, j+

Nc

z p ,q .

(2-c) The correct z-direction interface velocities {U

U

and

n+ 12

+ δz U iz,,jn+1 = Q in,+j 1 + f ( H ni ,+j 1 , C in,+j 1 ),

n +1 = − K zz i, j+ 1

n+ 1 φin,+j 1 C α ni ,+j 1 −φi , j 2

+ R α (C

2

implicit solution-flux coupled schemes as

i , jq + 2

2

∗n+ 12 ( Q n +1 C

(48)

+

z,n+1 } i , jq + 12

Nc

and { W α

z,n+1 } i , jq + 12

on the interface lines

γ pz,q are recomputed by



n +1

α = 1 ηα C α i , j q + 1 ,

(51)

2

n +1 n +1 n +1 ˆ n +1 1 )C α i , j q + (1 − λ( U z 1 ))C α i , j q +1 ) − D

i , jq + 2

i, j+ 2

C α ni ,+j 1+1 − C α ni ,+j 1 q

zz,i , j q + 12

q

z

,

(52)

and the x-direction numerical velocities and global concentration fluxes on sub-domains  p ,q are computed by

U

x,n+1 i + 12 , j

Wα

= − K xx n+11

x,n+1 i + 12 , j

i+ 2 , j

H ni ++11, j − H ni ,+j 1

x

,

(xi + 1 , z j ) ∈  p ,q ,

(53)

2

= U x n+11 (λ(U x n+11 )C α ni ,+j 1 + (1 − λ(U x n+11 ))C α ni ++11, j ) − Dˆ n+1 1 i+ , j i+ , j i+ , j xx,i + , j 2

2

2

2

C α ni ++11, j

−

x

C α ni ,+j 1

.

(54)

The important feature is that our scheme (39)–(54) is globally conservative over the whole domain since it is conservative both in multiple block-divided sub-domains and across inner interfaces, which is proved in the next section.

10

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

4. Mass conservation of our scheme In this section, we will prove our proposed scheme to satisfy mass conservation over the whole domain of domain decompositions. Theorem 1. (Mass conservation) The scheme (39)–(54) satisfies mass conservation in the discrete form over the whole domain of domain decompositions, i.e., for n ≥ 0



S s i , j H ni ,+j 1 x z =



i, j



S s i , j H 0 i , j x z +



i, j

φi0, j C α0 i , j x z +

(55)

n +1  

Q il, j C α∗ i , j x z t l

l =1 i , j

i, j n +1  

+

( Q il, j + f ( H li , j , C αl ,i , j )) x z t ,

l =1 i , j

i, j

φin,+j 1 C αn+i ,1j x z =

n +1  

( R α (C il , j ) + g α ( H li , j , C il , j )) x z t , α = 1, 2, · · · , N c .

(56)

l =1 i , j

Proof. We first rewrite Step 1. From (41), we have that for i = i p and i = i p + 1 n+ 12

S s i, j

H i, j

− H ni, j

t

+

Uα

x,n+ 12

− Uα

i + 12 , j

x,n+ 12 i − 12 , j

x

= 0,

i = i p , i p + 1.

(57)

When i = i p , using (41)c and (43), we have n+ 12

S s i, j

H i, j

− H nip , j

t

U

+

x,n+ 12

−U

i p + 12 , j

x,n+ 12 i p − 12 , j

x

U

=

x,n+ 12

1

x,n+ 2 − U˜ 1

i p + 12 , j

ip+ 2 , j

x

i = ip.

,

(58)

When i = i p + 1, using (41)c and (43), we have x,n+ 1

1

x,n+ 1

x,n+ 1

x,n+ 1

n+ 2 2 2 U − U 12 U − U˜ H i +21, j − H ni +1, j ip+ 2 , j i p + 32 , j i p + 12 , j i p + 12 , j p p S s i p +1 , j + =− , t x x

i = i p + 1.

(59)

Multiplying x z to (57)–(59), and summing them together, it gets that



n+ 1 S s i , j H i , j 2 x z

+

x,n+ 1

x,n+ 1

2

2

2 2  U i+ 1 , j − U i− 1 , j

i, j

x

i, j x,n+ 12

Noting the boundary conditions U 1

,j 2



n+ 1

S s i , j H i , j 2 x z =



i, j

x z =



S s i , j H ni, j x z.

(60)

i, j

= 0, U

x,n+ 12 I + 12 , j

= 0, we have

S s ni, j H α ni, j x z.

(61)

i, j

From (42), we have that for i = i p and i = i p + 1 n+ 1

n+ 1

φi , j 2 C α i , j 2 − φin, j C α ni, j

t

+

Wα

x,n+ 12

− Wα

i + 12 , j

x,n+ 12 i − 12 , j

x

= 0,

i = i p , i p + 1.

(62)

When i = i p , using (42)c and (44), we have n+ 1

n+ 1

φi p , j2 C α i , j 2 − φinp , j C α nip , j

t

+

Wα

x,n+ 12 i p + 12 , j

− Wα

x,n+ 12 i p − 12 , j

x

W

=

x,n+ 12 i p + 12 , j

1

˜ x,n+1 2 −W

ip+ 2 , j

x

,

i = ip.

(63)

When i = i p + 1, using (42)c and (46), we have n+ 1

φinp++11, j C α i p +21, j − φinp +1, j C α nip +1, j

t

+

Wα

x,n+ 12 i p + 32 , j

− Wα

x

x,n+ 12 i p + 12 , j

W

=−

x,n+ 12 i p + 12 , j

1

˜ x,n+1 2 −W

ip+ 2 , j

x

,

i = i p + 1.

(64)

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

11

Multiplying x z to (62)–(64), and summing them together, it gets that



n+ 12

φin,+j 1 C α i , j x z +

x,n+ 1

x,n+ 1

2

2

2 2  W α i+ 1 , j − W α i− 1 , j

i, j

x

i, j x,n+ 12

Noting the boundary conditions W α 1 ,j 2



n+ 12

φin,+j 1 C α i , j x z =



i, j

x z =



φin, j C α ni, j x z.

(65)

i, j x,n+ 1

= 0, W α 1 2 = 0, we have I+ , j 2

φin, j C α ni, j x z.

(66)

i, j

We then rewrite Step 2. From (49), we have that for j = jq and j = jq + 1 n+ 12

S s i, j

H ni ,+j 1 − H i , j

U

+

t

z,n+1 i , j + 12

− U z,n+11 i, j− 2

z

= Q in,+j 1 + f ( H ni ,+j 1 , C in,+j 1 ),

j = j q , j q + 1.

(67)

When j = jq , using (49)c and (51), we have n+ 12

S s i , jq

H ni ,+j 1 − H i , j q

U

q

+

t

z,n+1 i , j q + 12

− U z,n+11

i , jq − 2

z

=

1 Q in,+ jq

f ( H ni ,+j 1 , C in,+j 1 ) +

+

U

z,n+1 i , j q + 12

− U˜ z,n+11

i , jq + 2

z

,

j = jq .

(68)

When j = jq + 1, using (49)c and (51), we have z,n+1

1

z,n+1

z,n+1

z,n+1

n+ U −U U − U˜ H ni ,+j 1+1 − H i , j 2+1 i , j q + 12 i , j q + 32 i , j q + 12 i , j q + 12 q q S s i , j q +1 + = Q in,+jq1+1 + f ( H ni ,+j 1 , C in,+j 1 ) − , j = jq + 1. (69) t z z

Multiplying x z to (67)–(69), and summing them together, it gets that



S s i , j H ni ,+j 1 x z

i, j

+

z,n+1

z,n+1

2

2

 U i+ 1 , j − U i− 1 , j

x z =



n+ 1

S s H i , j 2 x z

z i, j i, j  n +1 n +1  n +1 + ( Q i , j + f ( H i , j , C i , j )) x z.

(70)

i, j

Noting the boundary conditions U



S s i , j H ni ,+j 1 x z =

i, j



z,n+1 i , 12

= 0, U z,n+11 = 0, we have i, J + 2

n+ 1

S s i , j H i , j 2 x z +

i, j

 ( Q in,+j 1 + f ( H ni ,+j 1 , C in,+j 1 )) x z.

(71)

i, j

Adding (61) (71), and summing all equations from 0 to n + 1, we obtain that



S s i , j H ni ,+j 1 x z =

i, j



S s i , j C α 0i , j x z +

n   1  l +1 ( Q il,+j 1 + f ( H li+ , j , C i , j )) x z t .

(72)

l =0 i , j

i, j

From (50), we have that for j = jq and j = jq + 1 z,n+1 z,n+1 n+ n+ Wα 1 − Wα 1 φin,+j 1 C α ni ,+j 1 − φi , j 2 C α i , j 2 i+ 2 , j i− 2 , j + = Q in,+j 1 C α∗ ni ,+j 1 + R α (C in,+j 1 ) t z + g α ( H ni ,+j 1 , C in,+j 1 ), j = jq , jq + 1. 1

1

(73)

When j = jq , using (42)c and (44), we have n+ 1

n+ 1

φin,+jq1 C α ni ,+jq1 − φi , jq 2 C α i , jq2

t

+

Wα

z,n+1 i , j q + 12

− W α z,n+11

i , jq − 2

z

+ g α ( H ni ,+j 1 , C in,+j 1 ) + When j = jq + 1, using (42)c and (46), we have

Wα

= Q in,+jq1 C α∗ ni ,+jq1 + R α (C in,+jq1 ) z,n+1 i , j q + 12

z,n+1 − W˜ α i , jq + 1

2

z

,

j = jq .

(74)

12

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

n+ 1

n+ 1

φin,+jq1+1 C α ni ,+jq1+1 − φi , jq 2+1 C α i , jq2+1

t +

+

Wα

z,n+1 i , j q + 32

− W α z,n+11

i , jq + 2

= Q in,+jq1+1 C α∗ ni ,+j 1 + R α (C in,+jq1+1 )

z Wα

g α ( H ni ,+j 1 , C in,+j 1 ) −

z,n+1 i , j q + 12

z,n+1

− W˜ α i , jq + 1

2

j = j q + 1.

,

z

(75)

Multiplying x z to (73)–(75), and summing them together, it gets that



φin,+j 1 C α ni ,+j 1 x z +

z,n+1 z,n+1  W α i, j+ 1 − W α i, j− 1 2

2

x z =



n+ 1

n+ 1

φi , j 2 C α i , j 2 x z

z i, j i, j  + ( Q in,+j 1 C α∗ ni ,+j 1 + R α (C in,+j 1 ) + g α ( H ni ,+j 1 , C in,+j 1 )) x z.

i, j

(76)

i, j z,n+1

Noting the boundary conditions W α 1 i, 2



φin,+j 1 C α ni ,+j 1 x z =

i, j



n+ 1

= 0, W α z,n+11 = 0, we have i, J + 2

n+ 1

φi , j 2 C α i , j 2 x z +

i, j

 ( Q in,+j 1 C α∗ ni ,+j 1 + R α (C in,+j 1 ) + g α ( H ni ,+j 1 , C in,+j 1 )) x z.

(77)

i, j

Adding (66) (77), and summing all equations from 0 to n + 1, we obtain that



φin,+j 1 C α ni ,+j 1 x z =

i, j



φi0, j C α 0i , j x z +

(78)

l =0 i , j

i, j

This ends the proof.

n   1 l +1  l +1  l +1 ( Q il,+j 1 C α∗ li+ , j + R α (C i , j ) + g α ( H i , j , C i , j )) x z t .

2

5. The linearization iterative algorithm of the mass conservative S-DDM scheme In the section, we develop the linearization iterative algorithm of the mass conservative S-DDM scheme (39)–(54) for solving the general coupled system of multicomponent contamination flows (19)–(24). Because of the coupling and nonlinearity of equations (41)–(42) and (49)–(50), that can be solved by Newton’s method at each sub-step, we can also propose a linearization iterative scheme. Considering nonlinear functions K xx (C ), K zz (C ), U x ( H , C ), U z ( H , C ), φ( H ), D xx (U x , U z ), D zz (U x , U z ), f ( H , C ), R (C ) and g α ( H , C ) to be required in a linearization in computation. We describe the linearization algorithm as below A. Initialization: Input the water head h and the multicomponent concentration c α at initial time t = t 0 and set H i0, j = h(t 0 , xi , z j ) and 0 C α i , j = c α (t 0 , xi , z j ), for each sub-domain  p ,q . B. Time Stepping Procedure: for the time steps t n , n = 0, 1, · · · , N − 1 do 1

1

1

1

B0. Set H n+ 2 ,0 = H n , C n+ 2 ,0 = C n , U x,n+ 2 ,0 = U x,n and W α x,n+ 2 ,0 = W α x,n for k = 0, and set the linearization index k = 1. B1. Along x-direction, 1

x,n+ 2 ,k (1-a) Compute the predict intermediate interface Darcy’s velocity U˜ and the predict intermediate interface global 1

flux W˜ α

x,n+ 12 ,k i p + 12 , j

ip+ 2 , j

on the interface lines

γ px,q by the linearized interface semi-implicit flux schemes n+ 1 ,k−1

2 K

1 t xx i p + 12 , j x,n x,n ˜ ˜ x,n+1 2 ,k + U x,n 3 , U =U 1 + U − 2 U 1 1 ip+ 2 , j ip+ 2 , j ip+ 2 , j ip− 2 , j ip+ 2 , j x2 S s i p + 1 , j 2 t x,n n+ 1 ,k−1 x,n+ 12 ,k φ 21 W˜ α = φin + 1 , j W α x,n 1 + U λ(U ixp,n, j ) W α x,n 1 1 ip+ 2 , j ip+ 2 , j ip+ 2 , j ip− 2 , j p 2 x i p , j

x,n+ 12 ,k

1

x,n+ ,k + (U ixp,n, j (1 − λ(U ixp,n, j )) − U ixp,n+1, j λ(U ixp,n+1, j )) W˜ α 12

− U ixp,n+1, j (1 − λ(U ixp,n+1, j )) W α x,n

i p + 32 , j



ip+ 2 , j

+

t ˆ n+ 12 ,k−1 x,n+ 1 ,k x,n D Wα − 2 W˜ α 12 + W α x,n 3 . 1 2 xx,i p + 1 , j i i − , j + , j i + , j p 2 p 2 p 2 2 x

(79)

(80)

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974 n+ 12 ,k

(1-b) Compute the intermediate interior solutions and fluxes { H i , j

 p ,q

x,n+ 12 ,k

13 n+ 1 ,k

x,n+ 1 ,k

} and {C α i , j 2 , W α 1 2 } in sub-domain i+ 2 , j by the linearized x-directional splitting solution-flux coupled scheme as, for (xi , z j ) ∈  p ,q n+ 1 ,k S s i, j H i, j 2

U U

x,n+ 12 ,k i + 12 , j

x,n+ 12 ,k i p + 12 , j

U

+ t

x,n+ 1 ,k x,n+ 1 ,k 2 −U 1 2 i+ 1 , j i− , j 2 2

x

i + 12 , j

(81)

= S s i , j H ni, j ,

n+ 1 ,k

= − K xx

,U

n+ 1 ,k 2

2 n+ 12 ,k−1 H i +1, j − H i , j x i+ 1 , j

(82)

,

2

1

x,n+ 2 ,k = U˜ , 1

(83)

ip+ 2 , j

and n+ 1 ,k n+ 1 ,k φi , j 2 C α i , j 2

Wα

x,n+ 12 ,k

=U

i + 12 , j

Wα

x,n+ 12 ,k

+ t

Wα

x,n+ 1 ,k x,n+ 1 ,k 2 −W α 1 2 i+ 1 , j i− , j 2 2

x

x,n+ 12 ,k i + 12 , j

(λ(U

x,n+ 12 ,k i + 12 , j

= φin, j C α ni, j ,

n+ 1 ,k )C α i , j 2

+ (1 − λ(U

x,n+ 12 ,k i + 12 , j

U

x,n+ 12 ,k i + 12 , j

ˆ D

n+ 1 ,k

n+ 12 ,k

C α i +12, j − C α i , j

n+ 12 ,k xx,i + 12 , j

x

,

(85)

1

(86)

ip+ 2 , j

(1-c) Recompute the correct intermediate Darcy’s velocity U

γ

n+ 1 ,k ))C α i +12, j ) −

x,n+ ,k = W˜ α 12 .

i p + 12 , j

interface lines

(84)

α = 1, 2, · · · , N c ,

x p ,q

x,n+ 1 ,k

and global concentration flux W α 1 2 i+ 2 , j

i + 12 , j

on both

and subdomains  p ,q by interior solutions n+ 1 ,k

=

x,n+ 12 ,k

n+ 1 ,k 2

2 n+ 1 ,k H i +1, j − H i , j − K xx 12 x i+ 2 , j

(87)

,

and

Wα

x,n+ 12 ,k i + 12 , j

=U

x,n+ 12 ,k i + 12 , j

(λ(U

x,n+ 12 ,k i + 12 , j

n+ 1 ,k )C α i p , j2

+ (1 − λ(U

x,n+ 12 ,k i + 12 , j

n+ 1 ,k ))C α i +12, j ) −

1

ˆ D

n+ 1 ,k

xx,i + 12 , j

1

n+ 12 ,k

C α i +12, j − C α i , j

n+ 12 ,k

1

x

.

(88) 1

B2. Along z-direction. With setting of H n+1,0 = H n+ 2 , C n+1,0 = C n+ 2 , U z,n+1,0 = U z,n+ 2 and W α z,n+1,0 = W α z,n+ 2 for k = 0. z,n+1,k (2-a) Compute the predict intermediate interface Darcy’s velocity U˜ and the predict intermediate interface global 1 i , jq + 2

z,n+1,k flux W˜ α i , j + 1 on the interface lines γ pz,q by the linearized interface semi-implicit flux schemes q 2

U˜

z,n+1,k i , j q + 12

=U

z,n+ 12 i , j q + 12

− t

+

K zz

n+1,k−1 1

i , jq + t 2 z2 S s i , jq + 1

n+1,k−1 K zz i , jq + 1 2 Ss i , jq + 1 2

2

U

z,n+ 12 i , j q − 12

− 2U˜ z,n+11,k + U i , jq + 2

z,n+ 12



i , j q + 32

(89)

δz ( Q n+1 + f ( H n+1,k−1 , C n+1,k−1 ))i , jq + 1 , 2

and 1

1

n+ z,n+ 2 z,n+1,k φ n+1,k1−1 W˜ α i , jq + 1 = φ 2 1 W α 1 + i , jq + 2

2

i , jq + 2

i , jq + 2



z,n+ 12 z,n+ 12 z,n+ 12 t U λ( U ) W α i , jq i , jq z i , jq − 1

2

z,n+ 1 z,n+ 1 +(U i , jq 2 (1 − λ(U i , jq 2 )) −

z,n+ 1 z,n+ 1 z,n+1,k U i , j +21 λ(U i , j +21 )) W˜ α i , j + 1 q q q 2

z,n+ 1 z,n+ 1 z,n+ 12 −U i , jq +21 (1 − λ(U i , jq +21 )) W α 3 i , jq + 2

1 z , n + z,n+ 12 n + 1 , k − 1 2 ˜ z,n+1,k + zt2 Dˆ 1 Wα 1 − 2 W α i , jq + 1 + W α 3 zz,i , j q + 2 i , jq − 2 i , jq + 2 2 n+1,k−1 n+1 ∗ n+1,k−1 n + 1 , k − 1 (Q C α + R α (C ) − t δz D zz

+ g α ( H n+1,k−1 , C n+1,k−1 )) . 1 i , jq + 2

(90)

14

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974 n+1,k

(2-b) Compute the interior solutions and fluxes { H i , j

, U z,n+11,k } and {C α ni ,+j 1,k , W α z,n+11,k } in sub-domain  p ,q by the i, j+ 2

i, j+ 2

linearized z-directional splitting solution-flux coupled scheme as, for (xi , z j ) ∈  p ,q n+1,k

S s i, j H i, j

U

+ t

z,n+1,k

i, j+ 1 2

−U z,n+11,k i, j−

2

z

U

z,n+1,k i , j + 12

= − K zz n+1,1k−1

U

z,n+1,k i , j q + 12

= U˜ z,n+11,k ,

n+ 12 ,k

H n+1,k − H n+1,k i , j +1

i, j+ 2

+ t ( Q in,+j 1 + f ( H ni ,+j 1,k−1 , C in,+j 1,k−1 )),

 c n+1,k−1 + αN= , 1 1 ηα C α

= S s i, j H i, j i, j

(91) (92)

i, j+ 2

z

(93)

i , jq + 2

and

φin,+j 1,k C α ni ,+j 1,k + t

Wα

z,n+1,k

i, j+ 1 2

− W α z,n+11,k i, j−

2

z

n+ 12 ,k

= φi , j

n+ 12 ,k

C α i, j

+ t Q in,+j 1 C α∗ ni ,+j 1,k−1

+ R α (C n+1,k−1 ) + g α ( H n+1,k , C n+1,k−1 ) , i, j

i, j

α = 1, 2, · · · , N c ,

Wα

z,n+1,k i , j + 12



= U z,n+11,k λ(U z,n+11,k )C α ni ,+j 1,k + (1 − λ(U z,n+11,k ))C α ni ,+j +1,1k − Dˆ n+1,k

Wα

z,n+1,k i , j q + 12

z,n+1,k = W˜ α i , jq + 1 .

i, j+ 2

i, j+ 2

γ

z p ,q

U

n+1,k

(95)

,

z

(96)

2

z,n+1,k i , j + 12

and global multicomponent concentration flux W α

z,n+1,k i , j + 12

on both interface

and subdomains  p ,q by the interior solutions

z,n+1,k i , j + 12

Wα

n+1,k

C α i , j +1 − C α i , j

zz,i , j + 12

i, j+ 2

(2-c) Recompute the Darcy’s velocity U lines

i, j

(94)

= − K zz n+1,1k

z,n+1,k i , j + 12

i, j+ 2

H n+1,k − H n+1,k i , j +1

i, j

z

+



Nc

n+1,k α = 1 ηα C α i , j + 1 ,

(97)

2

= U z,n+11,k (λ(U z,n+11,k )C α ni ,+j 1,k + (1 − λ(U z,n+11,k ))C α ni ,+j +1,1k ) − Dˆ n+1,k i, j+ 2

i, j+ 2

zz,i , j + 12

i, j+ 2

n+1,k

n+1,k

C α i , j +1 − C α i , j

z

.

(98)

B3. Set the iteration index k = k + 1. B4. Iterations. Repeat B1-B3 at current time step until Error < tolerance or k > k0 . Here we let Error = max{ Err1 , Err2 }, where

Err1 = H n+1,k − H n+1,k−1 ,

Err2 =

max

α =1,2··· , N c

C α n+1,k − C α n+1,k−1 . n+1,k

n+1 B5. Let H n+1 = H n+1,k , (U x , U z )n+1 = (U x , U z )n+1,k , and C α = Cα end do

, ( W αx , W αz )n+1 = ( W αx , W αz )n+1,k .

6. Numerical experiments In the section, we will show the performance of our mass conservative S-DDM scheme for solving the coupled nonlinear system of contamination flow in porous media. In Experiment 6.1, numerical results will first illustrate our scheme to be mass conservative, stable and second-order accuracy in space and first-order in time respectively, and we also show its parallelism. In Experiment 6.2, we will use our scheme to simulate the flow of four contaminants (Cl− , Na+ , Ca2+ and CO23− ) accompanied with precipitation reaction in porous media. In Experiment 6.3, we will use our scheme to simulate that the contaminant HCl flows through different porous media block of CaCO3 and porous media block of MgCO3 along with the dissolution reaction. 6.1. Experiment 1 Consider the following coupled system as

⎧ ∂h ⎪ ⎪ − ∇ · ( K (∇ h + ηc ∇ z)) = f , ⎪ ⎪ ⎪ ⎨ ∂t u = − K (∇ h + ηc ∇ z), ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ c + ∇ · (uc − D ∇ c ) = g , ∂t

(99)

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

15

Table 1 The space convergence orders for the water head H and the concentration C at t = 0.1, 0.2, 0.3, 0.5. h

1/50

1/100

1/200

1/300

H

enH Order

4.1665E-5 –

1.0443E-5 1.9963

2.6217E-6 1.9940

1.2098E-6 1.9074

C

enC Order

1.6328E-4 –

4.2946E-5 1.9268

1.0753E-5 1.9978

4.6120E-6 2.0878

H

enH Order

5.5678E-5 –

1.3971E-5 1.9947

3.4839E-6 2.0037

1.5615E-6 1.9792

C

enC Order

2.4465E-4 –

6.4080E-5 1.9328

1.5943E-5 2.0070

6.7884E-6 2.1057

H

enH Order

6.1786E-5 –

1.5509E-5 1.9942

3.8561E-6 2.0079

1.7104E-6 2.0049

C

enC Order

2.7979E-4 –

7.2996E-5 1.9385

1.8086E-5 2.0129

7.6782E-6 2.1130

H

enH Order

6.8742E-5 –

1.7239E-5 1.9955

4.2803E-6 2.0099

1.8867E-6 2.0204

C

enC Order

2.2816E-4 –

7.3006E-5 1.6440

1.8021E-5 2.0183

7.6693E-6 2.1070

t = 0.1

t = 0.2

t = 0.3

t = 0.5

Table 2 The time convergence orders for water head H and concentration C at t = 0.1, 0.2, 0.5, 1.0.

t

1/500

1/1000

1/2000

1/3000

H

enH Ratio

9.0889E-5 –

4.3556E-5 1.0612

2.1770E-5 1.0005

1.4559E-5 0.9923

C

enC Ratio

1.3664E-4 –

5.7621E-5 1.2457

3.1106E-5 0.8894

2.1797E-5 0.8771

H

enH Ratio

9.0168E-5 –

4.0815E-5 1.1435

2.0192E-5 1.0153

1.3505E-5 0.9920

C

enC Ratio

2.2173E-4 –

1.0245E-4 1.1139

5.5813E-5 0.8762

3.8916E-5 0.8894

H

enH Ratio

8.4306E-5 –

3.4272E-5 1.2986

1.6570E-5 1.0485

1.1089E-5 0.9906

C

enC Ratio

3.4462E-4 –

1.6945E-4 1.0241

8.9837E-5 0.9155

6.1623E-5 0.9297

H

enH Ratio

9.6723E-5 –

3.2115E-5 1.5906

1.3571E-5 1.2427

8.7688E-6 1.0771

C

enC Ratio

3.9531E-4 –

1.8911E-4 1.0638

9.4846E-5 0.9956

6.3569E-5 0.9868

t = 0.1

t = 0.2

t = 0.5

t = 1.0

with boundary conditions u · n |∂  = 0, ( D ∇ c ) · n |∂  = 0, where the domain is  = [0, 1] × [0, 1] and the physical parameters are K = 1, D = 0.02 and η = 0.5. The initial conditions are h0 (x, z) = 5x2 (1 − x)2 z2 (1 − z)2 and c 0 (x, z) = 200x2 (1 − x)2 z2 (1 − z)2 . The exact solutions of the water head and the concentration are h(x, z, t ) = 5e −t x2 (1 − x)2 z2 (1 − z)2 and c (x, z, t ) = 200e −t x2 (1 − x)2 z2 (1 − z)2 respectively. The right side functions f and g are given by substituting the exact solutions into the equations. Letting H ni, j and C in, j be the numerical solutions, we define numerical errors in discrete L 2 norm as



eH = (

i , j (h (xi , z j , t ) −

1

H i , j )2 h x h z ) 2 ,



eC = (

i , j (c (xi , z j , t ) −

1

C i , j )2 h x h z ) 2 ,

and errors of mass as

H Masserror = | C Masserror =



H ni, j h x h z − t

n 



0 l i, j f i, j hxh z − i , j H i , j h x h z |, l =1     | i , j C in, j h x h z − t ln=1 i , j g li , j h x h z − i , j C i0, j h x h z |. i, j

The domain  is divided into 2 × 2 block sub-domains. Take the uniform space steps h = h x = h z and the uniform time step t respectively. The convergence orders of water head H and concentration C in space step and in time step are listed in Tables 1 and 2. Numerical results at t = 0.1, 0.2, 0.3, and 0.5 are in Table 1, where the spatial step h is varied from 1/50 to 1/400 while the time step is taken as t = 1/100000. From Table 1, we can see clearly that our scheme is of second-order convergence in spatial step in L 2 -norm for the water head H , and the convergence order for the concentration C tends to two in L 2 -norm in the space step h by the second-order modified upwind technique. Numerical results of time convergence orders at time t = 0.1, 0.2, 0.5, 1.0 are in Table 2, where

16

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

Table 3 Mass errors of the water head H and the concentration C at t = 0.1, 0.2, 0.3, 0.5.

t

1/800

1/2000

1/3000

1/4000

t = 0.1

H C

Masserror Masserror

9.8879E-17 4.3091E-15

2.9195E-15 2.1407E-15

3.8216E-15 2.8345E-15

9.5497E-15 1.5335E-15

t = 0.2

H C

Masserror Masserror

2.6021E-18 1.6653E-16

5.3768E-15 9.0206E-16

7.3552E-15 1.7347E-16

1.8697E-14 1.3808E-15

t = 0.3

H C

Masserror Masserror

2.6108E-16 2.0678E-15

7.6631E-15 1.0270E-15

1.0184E-14 2.7339E-15

2.6611E-14 3.6982E-16

t = 0.5

H C

Masserror Masserror

1.2143E-16 4.5658E-15

1.1367E-14 1.6376E-15

1.5764E-14 4.8711E-15

4.0582E-14 1.5266E-16

Table 4 The stability for water head H and concentration C at t = 0.1, 0.2, 0.3, 0.5.

t

1/800

1/2000

1/3000

1/4000

t = 0.1

H C

enH enC

4.5752E-3 2.0453E-2

2.1814E-5 3.5309E-5

1.4576E-5 2.2938E-5

1.0943E-5 1.6817E-5

t = 0.2

H C

enH enC

3.6662E-3 1.8238E-2

2.0316E-5 6.2678E-5

1.3543E-5 4.0625E-5

1.0158E-5 2.9917E-5

t = 0.3

H C

enH enC

3.6242E-3 1.8064E-2

1.8973E-5 7.9246E-5

1.2645E-5 5.2014E-5

9.4730E-6 3.8413E-5

t = 0.5

H C

enH enC

2.3230E-3 1.0508E-2

1.6782E-5 9.5911E-5

1.1167E-5 6.3272E-5

8.3572E-6 4.6948E-5

Table 5 Errors and mass errors for the water head H and the concentration C over multi-block sub-domains. Sub-domains

2×2

3×3

4×4

5×5

H

enH Masserror

4.3943E-6 2.4980E-16

4.9078E-6 9.5844E-17

5.6916E-6 6.2667E-16

6.5841E-6 8.1196E-16

C

enc Masserror

6.0505E-6 4.1633E-15

3.1572E-5 5.6552E-16

3.3275E-5 1.2490E-15

3.7650E-5 2.8380E-15

H

enH Masserror

4.0788E-6 9.4282E-16

4.5746E-6 1.5179E-15

5.2815E-6 1.0517E-15

6.1027E-6 1.8232E-15

C

enc Masserror

1.0894E-5 1.6098E-15

3.3641E-5 3.5874E-15

3.5016E-5 2.7756E-17

3.8131E-5 1.0131E-15

H

enH Masserror

3.7940E-6 1.3652E-15

4.2663E-6 2.0383E-16

4.8820E-6 1.3884E-15

5.6197E-6 2.5435E-15

C

enc Masserror

1.4161E-5 2.8172E-15

3.4269E-5 1.5266E-16

3.4874E-5 3.2405E-15

3.7156E-5 4.2952E-15

H

enH Masserror

3.3530E-6 1.9555E-15

3.7767E-6 1.7044E-16

4.2373E-6 2.2994E-15

4.8327E-6 4.0198E-15

C

enc Masserror

1.7734E-5 6.0646E-15

3.4072E-5 2.5258E-15

3.3510E-5 3.4278E-15

3.5913E-5 5.8426E-15

t = 0.1

t = 0.2

t = 0.3

t = 0.5

h = 10 t, which shows clearly that our scheme is of first-order convergence for both water head H and concentration C in time step. In Table 3, the errors of mass for the water head H and the concentration C at time t = 0.1, 0.2, 0.3, and 0.5. Let h = 1/400. We can see clearly from Table 3 that the errors of mass for the water head H and concentration C reach accuracy 10−14 , which is machine precision. Thus, it shows that our scheme preserves mass for both water head H and concentration C over the whole domain of domain decompositions. The stability is tested in Table 4, where we let t from 1/200 to 1/4000 and take h = 1/400. From Table 4, we can 1 t see that even when t = 800 , i.e. r = = 200, our scheme is stable where the water head and concentration have good h2 accurate results. The performance of our scheme over more block-divided sub-domains is shown in Table 5, which shows mass conservation and numerical errors of our scheme over multiple sub-domains at t = 0.1, 0.2, 0.3 and 0.5, where we take t = 1/10000 and h = 1/400. Finally, in Table 6, the efficiency of our scheme over multi-block sub-domains is presented at t = 0.3. Define the speedup by S d = TT 1 , where T 1 is the CPU time, computed by the mass-preserving and second order splitting scheme over the domain d  without domain decomposition, T d is the CPU time, simulated by our scheme over d block-divided sub-domains. The

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

17

Table 6 Efficiency of our scheme over multi-block sub-domains.

t = 0.3

enH Masserror enC Masserror T p (sec) Sp E f f (%)

1

2×2

3×3

4×4

5×5

1.2306E-6 2.1669E-14 4.1202E-6 8.0005E-15 3140.6633 1.0000 100.00

1.2397E-6 1.6577E-14 4.1180E-6 1.0270E-15 802.4049 3.9141 97.85

1.2378E-6 1.3641E-14 4.8159E-6 9.0206E-16 384.7370 8.1631 90.70

1.2371E-6 1.0852E-14 4.5275E-6 8.4655E-16 217.2240 14.4582 90.36

1.2455E-6 8.0864E-15 4.6728E-6 1.3531E-15 137.1170 22.9050 91.62

Fig. 3. Geometry and boundary conditions for the vertical aquifer model. S

parallelism efficiency is defined by E f f = dd . In computation, we take the time step t = 1/40000 and the space step h = 1/400. From Table 6, it clearly shows that the speed up S p is smaller when the domain is divided into more block-divided sub-domains. Certainly, because of the communication between the sub-domains, the efficiency E f f will decrease when the more block-divided sub-domains are used. The efficiency E f f is 90% when the domain is divided into 5 × 5 sub-domains. The algorithm is implemented on the Server, which contained 14 nodes and twenty-eight cores is included in each node. 6.2. Experiment 2 In the multicomponent contamination flow, complex chemical reactions between multicomponent contaminants not only reduce the concentrations of the contaminant components, or even consume the contaminants (such as precipitation reactions and complexation reactions, etc.), but also increase the concentrations of the contaminants and produce the new contaminants (such as dissolution reactions, hydrolysis and oxidation-reduction reactions, etc.). In this experiment, we simulate the multicomponent contamination flow accompanied with precipitation reactions and dissolution reactions by two kinds of experiments. The reaction stoichiometry of the calcite crystal growth (see [4,13,15,22]) is

k Ca2+ + CO23− GGGGGGGGGGGGGGGGGA CaCO3 . precipitation

(100)

The solubility constant is K sp = 5 × 10−9 at 25 ◦ C and the rate of growth of calcium carbonate is r = kS e C Ca2+ C CO2− . The kinetic controlled reaction R for Ca2+ and CO23− can be described as

R Ca2+ = − v Ca2+ r ,

R CO2+ = − v CO2− r , 3

3

(101)

3

where the stoichiometric coefficients are v Ca2+ = v CO2− = 1. 3

The initial water flow contains the contaminant CO23− and the geometry of the vertical two dimensional model and the boundary conditions are listed in Fig. 3. The horizontal and vertical dimensions of the aquifer are L = 1000 m and Z = 20 m respectively. An initial steady horizontal flow along x-direction is created by imposing a hydraulic head difference between the left boundaries hl = 20 m and right boundaries hr = 10 m. The Neumann boundary conditions ∂∂hz = 0

are assumed on the top and bottom boundaries. The injected multicomponent contaminants (Cl− and Ca2+ ) are located ∗ , and the Neumann boundary between x1 = 50 m and x2 = 75 m by imposing a Dirichlet boundary condition on c α = c α conditions are specified at the other boundaries. The main parameters in simulation are the porosity φ = 0.35, the density of groundwater ρ0 = 1 × 103 kg/m3 , the x-directional permeability kxx = 8 × 10−11 m2 and z-directional k zz = 2 × 10−11 m2 , the storage coefficient S p = 10−5 , the fluid viscosity μ = 1 × 10−3 kg m−1 s−1 , the averaged free-liquid diffusion coefficient D ∗ = 5.0 × 10−3 m2 /day, the transverse dispersivity dl = 1 m, and the longitudinal dispersivity dt = 0.001 m. We choose ηCl− = 1 × 10−3 L/mol, ηCa2+ = 1 × 10−3 L/mol and ηCO2− = 1 × 10−3 L/mol in the test. The injected concentration of 3

the contaminants Cl− and Ca2+ are c ∗ − = 0.4 mol/L and c ∗ Cl

Ca2+

= 0.03 mol/L, respectively. The initial CO23− concentration is

c CO2− = 0.005 mol/L. The reaction rate is set as k = 300 mol m−2 day−1 . The effective surface area is taken as S e = 100 m−1 . 3

18

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

Fig. 4. Concentration distribution of contaminant Cl− after t = 2, 5, 8, 11 months. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

We assume that the region is divided into 2 × 2 sub-domains. The domain is further discretized by the space step h x = 10 m and h z = 0.4 m, and the time step is set t = 1 day. In the test, we consider that the contaminants Ca2+ and Cl− are with intermittent injection on [x1 , x2 ], by imposing a Dirichlet boundary condition on [x1 , x2 ], i.e., the contaminants are injected every 3 months. From Fig. 4(a), the contaminant Cl− transports along the flow at t = 2 months. After t = 3 months, the injected contaminant Cl− concentration becomes zero and transports along the flow under velocity and dispersion (see Fig. 4(b)). When at t = 6 months, due to the refilled contaminant Cl− , there will produce the new concentration distribution near the injected point (see Fig. 4(c)). After t = 9 months, the two contaminant cells will transport along the flow, respectively (see Fig. 4(d)). In Fig. 5(b), due to the contaminant Ca2+ intermittent injecting at t = 3 months, the injected contaminant Ca2+ is depleted in presence of the contaminant CO23− due to calcite precipitation. Besides, the Ca2+ concentration is higher than

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

Fig. 5. Concentration distribution of contaminant Ca2+ after t = 2, 5, 8, 11 months.

Table 7 Parameter values for numerical simulations. Property

Value

Porosity (φ ) Permeability (k xx ) Permeability (k zz ) Longitudinal dispersivity (dl ) Transverse dispersivity (dt ) Diffusion (D ∗ ) Domain length Domain width

0.3 8 × 10−11 m2 4 × 10−11 m2 3.0 × 10−4 m 0 2 × 10−3 m2 /day 2000 m 25 m

19

20

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

Fig. 6. Concentration distribution of contaminant CO23− after t = 2, 5, 8, 11 months.

Fig. 7. Geometry and boundary conditions.

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

21

Fig. 8. The concentration distribution of contaminant Ca2+ after t = 6, 12, 18, 24 months.

CO23− , which leads to the excess Ca2+ in groundwater and the excess contaminant Ca2+ will transport along flow and will

be consumed by the contaminant CO23− , which leads to Ca2+ decrease. After t = 6 months, the refilled contaminant Ca2+ increases the Ca2+ concentration in groundwater (see Fig. 5(c)). At t = 9 months, the contaminant repeats the process as Fig. 5(b) (see Fig. 5(d)). In Fig. 6, because the contaminant CO23− is consumed by the contaminant Ca2+ which causes CO23− depletes in groundwater where the contaminant Ca2+ flows. The results indicate that the multicomponent contamination transport along the flow is accompanied by calcite precipitation reaction. Comparing Fig. 5 with Fig. 6, we can see that the contaminant CO23− will deplete if and only the contaminant Ca2+ exists.

22

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

Fig. 9. The concentration distribution of contaminant Mg2+ after t = 12, 18, 24 months.

6.3. Experiment 3 The simulations are presented to show the effects of acid on the advancement of carbonate mineral dissolution front. As known, the kinetics of heterogeneous carbonate mineral dissolution on pH < 7 is described as

k MCO3 + 2H+ GGGGGGGGGGGGGGGA M 2+ + H2 CO3 , dissolution

(102)

and the rate of reaction depends on the concentration of hydrogen ions, or the activity of H+ (aH+ ) (see [28]), i.e.

r = k± S e aH+ .

(103)

In the following, the investigated numerical trials take place in a homogeneous, saturated carbonate minerals medium, which includes 10% calcium carbonate (CaCO3 ), 10% magnesium carbonate (MgCO3 ) and other insoluble porous media in acid. The CaCO3 is located between x1 = 400 m and x2 = 600 m, and MgCO3 is located between x3 = 1200 m and x4 = 1400 m (see Fig. 7). The source contaminant HCl measures 3 m and is located 5 m from the top of the domain on the left hand side of the porous media tank. The boundaries at the top and bottom of the tank are impermeable. The left and right boundaries are modeled as specified equivalent freshwater heads (hl = 25 m, hr = 5 m). Boundary conditions for concentration are assigned as Dirichlet boundary at the left side of the porous media tank, and all other concentration boundaries are assigned as zero-dispersive flux. The physical parameters used for the simulations are shown in Table 7. Take ηCa2+ = 8 × 10−3 L/mol, ηMg2+ = 4 × 10−3 L/mol, ηCl− = 8 × 10−3 L/mol, ηH2 CO3 = 1.2 × 10−2 L/mol and ηH+ =

2 × 10−4 L/mol. The reaction rates are kCaCO3 = 5 mol m−2 day−1 for CaCO3 and kMgCO3 = 20 mol m−2 day−1 for MgCO3 ,

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

23

Fig. 10. The concentration distribution of contaminant H2 CO3 after t = 6, 12, 18, 24 months.

respectively. The reactive surface areas for CaCO3 and MgCO3 are S e = 1 m−1 . In simulation, we assume that the porosity doesn’t vary with dissolution reaction and the weak acid H2 CO3 is soluble in water. The rectangular domain is divided into 2 × 2 sub-domains, and the space steps are h x = 20 m and h z = 0.5 m and the time step is t = 1 day. In Fig. 8, when the contaminant H+ arrives at [x1 , x2 ], i. e., flows into the CaCO3 media, the contaminant H+ reacts with CaCO3 which produces the new contaminants Ca2+ . The component Ca2+ will increase with dissolution reaction and will move along the flow. In Fig. 9, when the contaminant H+ transports through [x3 , x4 ] (MgCO3 media), the contaminant H+ will react with MgCO3 and produce the contaminant Mg2+ and H2 CO3 . The concentration of Mg2+ will then increase with time due to the dissolution reaction between the MgCO3 media with the contaminant H+ and will move along the flow. In Fig. 10, when the contaminant H+ arrives at [x1 , x2 ], i.e., flows into the CaCO3 media, the contaminant H+ reacts with CaCO3 which produces the new contaminants Ca2+ and H2 CO3 . The component H2 CO3 will increase with dissolution

24

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

Fig. 11. The concentration distribution of contaminant Ca2+ after t = 12, 18, 24 months.

Fig. 12. The concentration distribution of contaminant Mg2+ after t = 18, 24 months.

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

25

Fig. 13. The concentration distribution of contaminant H2 CO3 after t = 12, 18, 24 months.

reaction and will move along the flow (Fig. 10(a)(b)). When the contaminant H+ transports through [x3 , x4 ] (MgCO3 media), the contaminant H+ will react with MgCO3 and produce the new contaminant Mg2+ and the new H2 CO3 (Fig. 10(c)). The concentration of H2 CO3 will increase bigger due to the moving of the concentration of H2 CO3 produced in the dissolution reaction between H+ with CaCO3 in the first part [x1 , x2 ]. After t = 12 months, the concentration of H2 CO3 accumulate and move along the flow (Fig. 10(d)). Finally, we consider the porous media is not homogeneous, i.e., where x ≤ 14 m, the permeability in horizontal direction is kxx = 1.5 × 10−10 m2 , when x > 14 m, it is kxx = 5 × 10−11 m2 , while the permeability in vertical direction is k zz = 1.5 × 10−10 m2 . Numerical results of the concentrations of multicomponent contaminants are shown after t = 6, 12, 18, 24 months in Figs. 11–13. From Fig. 11, due to the bigger permeability horizontal direction in bottom region, the contaminant H+ will first react with CaCO3 and MgCO3 , it leads to the concentration of Ca2+ is bigger in bottom region than in upper region. In Fig. 12, the contaminant H+ first reacts with MgCO3 and produces Mg2+ in bottom region. From Fig. 13, there is a superposition phenomenon about H2 CO3 in bottom region, that is because there exists dissolution reaction between H+ and MgCO3 , which leads Mg2+ increase in bottom region. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements D. Liang’s work was supported partially by Natural Sciences and Engineering Research Council of Canada. Z. Zhou’s work was supported partially by Natural Science Foundation of China (under grant 61703250), Natural Science Foundation of Shandong Province (under grant ZR2017BF002) and Shandong Agricultural University (under grant xxxy201704). The authors would thank referees for their comments which have helped to improve the paper.

26

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

References [1] S. Bea, J. Carrera, C. Ayora, F. Batlle, M. Saaltink, CHEPROO: a Fortran 90 object-oriented module to solve chemical processes in Earth science models, Comput. Geosci. 36 (2009) 1098–1112. [2] B. Bekins, B. Rittmann, J. MaDonald, Natural attenuation strategy for ground water cleanup focuses on demonstrating cause and effect, Eos Trans. AGU 82 (2001) 53–58. [3] Y. Bonbendir, X. Antonie, C. Geuzaine, A quasi-optimal non-overlapping domain decomposition algorithm for Helmholtz equations, J. Comput. Phys. 231 (2012) 262–280. [4] L. Chou, M. Garrels, R. Wollast, Comparative study of the kinetics and mechanisms of dissolution of carbonate minerals, Chem. Geol. 78 (1989) 269–282. [5] C.N. Dawson, T.F. Dupont, Explicit/implicit conservative domain decomposition procedures for parabolic problems based on block centered finite differences, SIAM J. Numer. Anal. 31 (1994) 1045–1061. [6] C.N. Dawson, Q. Du, T.F. Dupont, A finite difference domain decomposition algorithm for numerical solution of heat equations, Math. Comput. 57 (1991) 63–71. [7] M. Dryja, X. Tu, A domain decomposition discretization of parabolic problems, Numer. Math. 107 (2007) 625–640. [8] C. Du, D. Liang, An efficient S-DDM iterative approach for compressible contamination fluid flows in porous media, J. Comput. Phys. 229 (2010) 4501–4521. [9] Q. Du, M. Mu, Z. Wu, Efficient parallel algorithms for parabolic problems, SIAM J. Numer. Anal. 39 (2001) 1469–1487. [10] V. Freedman, M. Ibaraki, Effects of chemical reactions on density-dependent fluid flow: on the numerical formulation and the development of instabilities, Adv. Water Resour. 25 (2002) 439–453. [11] T. Graf, R. Therrien, Variable-density groundwater flow and solute transport in porous media containing nonuniform discrete fractures, Adv. Water Resour. 28 (2005) 1351–1367. [12] T. Graf, R. Therrien, Coupled thermohaline groundwater flow and single-species reactive solute transport in porous media, Adv. Water Resour. 30 (2007) 742–771. [13] M. Hayek, G. Kosakowski, A. Jakob, V. Churakov, A class of analytical solutions for multidimensional multispecies diffusive transport coupled with precipitation-dissolution reactions and porosity changes, Water Resour. Res. 48 (2012) W03525. [14] H. Henderson, K. Mayer, L. Parker, A. Al, Three-dimensional density-dependent flow and multicomponent reactive transport modeling of chlorinated solvent oxidation by potassium permanganate, J. Contam. Hydrol. 106 (2009) 195–211. [15] W. Inskeep, P. Bloom, An evaluation of rate equation for calcite precipitation kinetics at P CO2 less than 0.01 atm and pH greater than 8, Geochim. Cosmochim. Acta 49 (1985) 2165–2180. [16] H. Jansen, R. Zeebe, D. Wolf-Gladrow, Modeling the dissolution of settling CaCO3 in the ocean, Glob. Biogeochem. Cycles 16 (2002). [17] J. Kanney, C. Miller, D. Barry, Comparison of fully coupled approaches for approximating nonlinear transport and reaction problems, Adv. Water Resour. 26 (2003) 353–372. [18] D. Kirkner, H. Reeves, Multicomponent mass transport with homogeneous and heterogeneous chemical reactions: effect of the chemistry on the choice of numerical algorithm 1. Theory, Water Resour. Res. 24 (1988) 1719–1729. [19] D. Lanser, J.G. Verwer, Analysis of operator splitting for advection–diffusion–reaction problems from air pollution modelling, J. Comput. Appl. Math. 111 (1999) 201–216. [20] A. Lasaga, Kinetic Theory in the Earth Sciences, Princeton Univ. Press, Princeton, N.J., 1998. [21] R. Lazarov, I. Mishev, P. Vassilevski, Finite volume methods for convection-diffusion problems, SIAM J. Numer. Anal. 33 (1996) 31–55. [22] L. Li, C. Steefel, L. Yang, Scale dependence of mineral dissolution rates within single pores and fractures, Geochim. Cosmochim. Acta 72 (2008) 360–377. [23] D. Liang, C. Du, An efficient S-DDM scheme and its analysis for solving parabolic equations, J. Comput. Phys. 272 (2014) 46–69. [24] P. Lichtner, Modeling of Reactive Flow and Transport in Natural Systems, Paper presented at the Rome Seminar on Environmental Geochemistry, Castelnuovo di Porto, Rome, 1996. [25] D. Liang, W. Zhao, An optimal weighted upwind covolume method on non-standard grids for convection-diffusion problems in 2D, Int. J. Numer. Methods Eng. 67 (2006) 553–577. [26] P. Lions, On the Schwarz alternating method I, in: R. Glowinski, G. Golub, G. Meurant, J. Periaux (Eds.), Proceeding of the First International Symposium on Domain Decomposition Methods for Partial Differential Equations, Paris, 1987, SIAM, Philadelphia, USA, 1988, pp. 1–42. [27] P. Lions, On the Schwarz alternating method II, in: T. Chan, R. Glowinski, J. Periaux, O. Wildlund (Eds.), Proceeding of the Second International Symposium on Domain Decomposition Methods for Partial Differential Equations, Los Angeles, 1988, SIAM, Philadelphia, USA, 1988, pp. 47–70. [28] C. Liu, T. Narasimhan, Redox-controlled multiple-species reactive chemical transport 1. Model development, Water Resour. Res. 25 (1989) 869–882. [29] K. Mayer, O. Emil, W. David, Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation kinetically controlled reactions, Water Resour. Res. 38 (2002) 1174. [30] X. Mao, H. Prommer, D. Barry, C. Langevin, B. Panteleit, L. Li, Three-dimensional model for multi-component reactive transport with variable density groundwater flow, Environ. Model. Softw. 21 (2006) 615–628. [31] W. Rivera, J. Zhu, D. Huddleston, An efficient parallel algorithm with application to computational fluid dynamics, J. Comput. Math. Appl. 45 (2003) 165–188. [32] M. Saaltink, C. Ayora, J. Carrera, On the behavior of approaches to simulate reactive transport, J. Contam. Hydrol. 48 (2001) 213–235. [33] X. Sanchez-Vila, L. Donado, A. Guadagnini, J. Carrera, A solution for multicomponent reactive transport under equilibrium and kinetic reactions, Water Resour. Res. 46 (2010) W07359. [34] Z. Sheng, G. Yuan, X. Hang, Unconditional stability of parallel difference schemes with second order accuracy for parabolic equation, Appl. Math. Comput. 184 (2007) 1015–1031. [35] H. Shi, H. Liao, Unconditional stability of corrected explicit/implicit domain decomposition algorithms for parallel approximation of heat equations, SIAM J. Numer. Anal. 44 (2006) 1584–1611. [36] M. Simpson, T. Clement, Theoretical analysis of the worthiness of Henry and Elder problems as benchmarks of density dependent groundwater flow models, Adv. Water Resour. 26 (1995) 17–31. [37] A. Toselli, O. Widlund, Domain Decomposition Methods-Algorithms and Theory, Springer, 2005. [38] M. Tran, Overlapping optimized Schwarz methods for parabolic equations in n-dimensions, Proc. Am. Math. Soc. 141 (2013) 1627–1640. [39] P. Vabishchevich, Domain decomposition methods with overlapping subdomains for the time-dependent problems of mathematical physics, Comput. Methods Appl. Math. 8 (2008) 393–405. [40] A. Walter, E. Frind, D. Blowes, C. Ptacek, J. Molson, Modeling of multicomponent reactive transport in groundwater 1. Model development and evaluation, Water Resour. Res. 30 (1994) 3137–3148. [41] A. Walter, E. Frind, D. Blowes, C. Ptacek, J. Molson, Modeling of multicomponent reactive transport in groundwater 2. Metal mobility in aquifers impacted by acidic mine tailings discharge, Water Resour. Res. 30 (1994) 3149–3158. [42] T. Xu, E. Sonnenthal, N. Spycher, K. Pruess, TOUGHREACT user’s guide: a simulation program for non-isothermal multiphase reactive geochemcial transport in variably saturated geologic media: applications to geothermal injectivity and CO2 geological sequestration, Comput. Geosci. 32 (2006) 145–165.

D. Liang, Z. Zhou / Journal of Computational Physics 400 (2020) 108974

27

[43] H. Zhang, F. Schwartz, W. Wood, S. Garabedian, D. LeBlanc, Simulation of variable-density flow and transport of reactive and non-reactive solutes during a tracer test at Cape Cod, Massachusetts, Water Resour. Res. 34 (1998) 67–82. [44] Z. Zheng, B. Simenon, L. Petzold, A stabilized explicit Lagrange multiplier based domain decomposition method for parabolic problems, J. Comput. Phys. 227 (2008) 5272–5285. [45] Z. Zhou, D. Liang, The mass-preserving S-DDM scheme for two-dimensional parabolic equations, Commun. Comput. Phys. 19 (2016) 411–441. [46] S. Zhu, Conservative domain decomposition procedure with unconditional stability and second-order accuracy, Appl. Math. Comput. 216 (2010) 3275–3282. [47] L. Zhu, G. Yuan, Q. Du, An explicit-implicit predictor-corrector domain decomposition method for time dependent multi-dimensional convection diffusion equations, Numer. Math. Theor. Meth. Appl. 2 (2009) 1–25. [48] Y. Zhuang, X. Sun, Stabilitized explicit-implicit domain decomposition methods for the numerical solution of parabolic equations, SIAM J. Sci. Comput. 24 (2002) 335–358.