Simulation of reactive transport in porous media using radial point collocation method

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Simulation of reactive transport in porous media using radial point collocation method Aatish Anshuman a, T.I. Eldho a,∗, Laishram Guneshwor Singh b a b

Department of Civil Engineering, Indian Institute of Technology (IIT) Bombay, Mumbai 400076, India Health Physics Division, Bhabha Atomic Research Centre, Trombay, Mumbai, India

a r t i c l e

i n f o

Keywords: Point collocation method Radioactive decay Reactive transport Radial basis function Linear adsorption

a b s t r a c t In this paper, a meshfree Radial Point Collocation method (RPCM) for modelling of reactive transport involving first-order decay and adsorption is presented. Unlike conventional grid-based methods such as Finite Element (FEM) and Finite Difference Methods (FDM), it does not require nodal connectivity information for domain representation. The proposed method does not necessitate operator splitting for accommodating reactions as it is done in FEM and FDM. The model domain is divided into intersecting circular shaped local support domains. This technique solves the ill-conditioning problem associated with globally supported system and allows use of more number of nodes which increases the accuracy of models. The multi-quadrics radial basis function (MQ-RBF) is used for the approximation/interpolation of solute concentration values in the local support domain. 1D and 2D models are developed and verified against benchmark analytical solutions. Sensitivity analyses are performed for shape parameters of the MQ-RBF, support domain size and nodal density. Two case studies are presented with regular and irregular shaped domains and results are compared with FEM. The results demonstrate advantages of the proposed model over grid-based methods while handling adsorption and first-order decay reactions. This study shows that the proposed model is effective for simulating reactive transport problems in porous media.

1. Introduction Groundwater is increasingly being used as the primary source of water for drinking, agriculture and industrial uses in many parts of the world. However, the quality of groundwater is being compromised due to various natural and anthropogenic pollution sources. Some of the factors include deliberate or, inadvertent spills from pollution sources such as landfills, industrial facilities, septic tanks, agricultural lands and nuclear waste storage facilities [1]. After entering the groundwater, the transport or migration of these pollutants are governed by processes such as advection and dispersion. In addition, reactive contaminants undergo through various other processes such as degradation, adsorption, production, etc. [2]. In particular, contaminants involving radioactive elements undergo first-order decay which emits harmful radioactive emissions. Such emissions include alpha, beta and gamma rays and also neutrons which ionize the medium through which they pass. In the absence of proper barriers, these radioactive elements can find their way to groundwater from waste disposal sites [3]. If the contaminated groundwater is ingested or, comes into physical contact, it may cause radiological damage to body organs and may eventually lead to cancer with prolonged exposure. Other organic pollutants such as BTEX compounds and chlorinated hydrocarbons may enter groundwater through ∗

contaminated lakes/rivers, leakages from petroleum industries, widespread use of insecticides and from sewage [4]. These contaminants are soluble in water [5,6] and may also get adsorbed in the soil mass [7]. Chlorinated hydrocarbons are known to affect the nervous system [8] and remain in the food chain for a long duration as they are highly soluble in fats [9]. The fate of such contaminants and their movement in the aquifer needs to be studied/tracked to protect the population from their harmful effects. This will enable us to assess the impact of contamination, planning of such industries or facilities, proper planning of aquifer remediation, etc. The finite difference and finite element based models are the conventional techniques for the modelling of flow and transport problems in porous media [2]. These methods are grid or, mesh based. For many problems, mesh constructions are very complex and also computationally expensive. Performing adaptive analysis, which requires repeated re-meshing of the entire domain, with such mesh-based methods may become prohibitively difficult [10]. Numerical dispersion and instability of the solution while handling highly advective flows [11,12] are also well-known issues with these conventional methods. This issue can be alleviated by creating very fine meshes, but at the cost of making the problem computationally very expensive [13]. For these reasons, many Eulerian, particle-based and mixed mesh-free

Corresponding author. E-mail address: [email protected] (T.I. Eldho).

https://doi.org/10.1016/j.enganabound.2019.03.016 Received 25 September 2018; Received in revised form 6 February 2019; Accepted 13 March 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

numerical methods are actively being developed by researchers [14–18]. In meshfree methods, shape functions are used to approximate the state variable across the domain. Some examples of these shape functions include Weighted Least Squares (WLS), Moving Least Squares (MLS), Polynomial Point Interpolation, Radial basis functions (RBF) etc. [10]. MQ-RBFs were first introduced by Hardy [19] for data interpolation in surveying. Kansa [20,21] introduced the applicability of MQRBFs to solve partial differential equations (PDEs) using point collocation methods. Kansa’s method involved the formation of a single global support domain which leads to a densely populated system matrix. This formulation leads to an increase in the rank of the matrix and yields an ill-conditioned matrix. Subsequently, Kansa and Hon [22] suggested some ways to alleviate the issue of globally supported domains. Another proposed way to deal with such problems is the formation of several overlapping local support domains [23,24]. The use of local support domains ensures a sparse system matrix and curbs the ill-conditioning problems efficiently. The support domain shape and size can be fixed or, varying for a given problem [10]. RBF based approximations are used for developing many models to simulate flow and contaminant transport in groundwater using point collocation method. Meenal and Eldho [17,18] used a rectangular shaped local support domain for modelling two-dimensional groundwater flow and transport problems. They used dummy nodes around the boundary to keep the support domain size uniform. Singh et al. [25] used a circular shaped local support domain to model a coupled groundwater flow and transport problem. The developed model was successfully used in solving inverse problems with simulation-optimization approach viz., contaminant source identification [26] and aquifer parameter estimation [27] by linking it to swarm intelligence optimization techniques. In above-mentioned studies, the solute transport being modelled is assumed to be conservative in nature without degradation. Apart from these, Kumar and Dodagoudar [28] developed models based on element free Galerkin method using MLS approximations to simulate reactive transport in porous media. However, the model uses back ground mesh to compute concentrations using MLS approximations. In this study, we developed a model based on RPCM to simulate transport in porous media involving reactions such as first-order decay and linear adsorption. The model is truly meshfree which does not use any background mesh for approximation. MQ-RBF is used to approximate shape functions and its spatial derivatives. The Neumann boundary condition is implemented using the direct collocation technique [10]. The model is compared against benchmark/standard analytical solutions and with FEM for short-lived solutes having fast decay rates i.e. in order of days and found to be effective.

point sources or sinks. ne is the effective porosity of the media and n is the direction perpendicular to the plane considered. Dxx and Dyy are dispersion coefficients in x- and y-direction, respectively, computed as [4]: 𝐷𝑥𝑥 =

𝛼𝐿 𝑉𝑥 2 + 𝛼𝑇 𝑉𝑦 2 𝑉̄

(2a)

𝐷𝑦𝑦 =

𝛼𝐿 𝑉𝑦 2 + 𝛼𝑇 𝑉𝑥 2 𝑉̄

(2b)

where 𝛼 L and 𝛼 T are longitudinal and transverse dispersivities, respec√ tively, and 𝑉̄ = 𝑉𝑥 2 + 𝑉𝑦 2 . ΣRxn is the reaction term in Eq. (1a) which represent different reactions that the dissolved solute undergoes. In this study, we considered first-order decay and linear adsorption. For solute undergoing first-order decay, the reaction rate is given by Batu [1]: 𝐾𝑑𝑒𝑐𝑎𝑦 = −𝜆𝐶𝑓 𝑜

(3)

Where Cfo is the concentration and 𝜆 is the rate of decay of the solute undergoing first-order decay. The half-life (t1/2 ) is the time elapsed for the activity or concentration of element to reduce to half of its initial value. The decay rate 𝜆 is calculated from the half-life as [1]: 𝜆=

0.693 𝑡1∕2

(4)

In the case of linear equilibrium adsorption, the adsorbed mass (S) is linearly proportional to the solute concentration (Cla ) and it is computed as the fraction of solute concentration as 𝑆 = 𝑘𝑑 𝐶𝑙𝑎

(5)

where kd is known as the distribution coefficient which is measured experimentally for different solute species [4]. The reaction rate for the solute undergoing linear equilibrium adsorption is given by Rastogi [4] 𝜌 𝜕𝑆 𝐾𝑠𝑜𝑟𝑝𝑡𝑖𝑜𝑛 = − 𝑏 (6) 𝑛𝑒 𝜕𝑡 Considering both reactions, the ΣRx n term for solute species C is derived using Eqs. (3), (5) and (6) as ∑ 𝜌 𝑘 𝜕𝐶 𝑅𝑥𝑛 = − 𝑏 𝑑 − 𝜆𝐶 (7) 𝑛𝑒 𝜕𝑡 By substituting above term in Eq. (1a) and rearranging, the governing equation for reactive transport in groundwater is given by ( ) ( ) [ ] [ ] 𝜕 𝑉𝑦 𝐶 𝜕 𝑉𝑥 𝐶 𝜕𝐶 𝜕 𝜕𝐶 𝜕 𝜕𝐶 𝐷 + 𝑅 = 𝐷 − − 𝜕𝑡 𝜕𝑥 𝑥𝑥 𝜕𝑥 𝜕𝑦 𝑦𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦 )( ) 𝑞𝑤 𝐶 ′ ( (8) − 𝜆𝐶 − 𝛿 𝑥 − 𝑥𝑖 𝑦 − 𝑦𝑖 𝑛𝑒

2. Governing equations and boundary conditions

Where 𝑅 = 1 +

The governing equation for transient solute transport in groundwater in two dimensions is given by [2,4]. ( ) ( ) [ ] [ ] 𝜕 𝑉𝑦 𝐶 𝜕 𝑉𝑥 𝐶 𝜕𝐶 𝜕 𝜕𝐶 𝜕 𝜕𝐶 = 𝐷 − − 𝐷 + 𝜕𝑡 𝜕𝑥 𝑥𝑥 𝜕𝑥 𝜕𝑦 𝑦𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦 ′ ∑ )( ) 𝑞𝑤 𝐶 ( (1a) 𝑅𝑥𝑛 + + 𝛿 𝑥 − 𝑥𝑖 𝑦 − 𝑦𝑖 𝑛𝑒

𝜌𝑏 𝐾𝑑 𝑛𝑒

, which is known as the retardation factor. It may

be noted that when only linear adsorption is considered, the governing equation can be solved alternatively by using modified velocities and dispersion coefficients which are obtained by dividing original velocities and dispersion coefficients by a factor R. 3. RPCM formulation

Initial condition for the transport problem is given by 𝐶 (𝑥, 𝑦, 0) = 𝐶0 𝑥, 𝑦 ∈ Ω

3.1. Construction of radial basis function and derivatives

(1b)

The general form of boundary conditions is: Dirichlet Boundary ∶ Neumann Boundary ∶

𝐶 (𝑥, 𝑦, 𝑡) = 𝑔 (𝑥, 𝑦, 𝑡) 𝜕𝐶 𝜕𝑛

= 𝑓 (𝑥, 𝑦, 𝑡)

𝑥, 𝑦 ∈ 𝜕 Ω1

𝑛 ∈ 𝜕 Ω2

In the radial point collocation method (RPCM), the problem domain is represented as a set of nodes scattered throughout the domain. Radial basis functions are employed for function interpolation while the governing equations and boundary conditions (Eq. (1)) are discretized at each of the nodes using point collocation method [18,25]. For interpolation of scalar variable C, a combination of ‘n’ radial basis functions (𝜙) and ‘m’ polynomial basis functions (p) can be used as follows [10]:

(1c) (1d)

where Vx and Vy are seepage velocities in x- and y-directions, C is the concentration of the solute in the groundwater and Ć is the concentration of the solute species injected or extracted at the rate of qw by the 9

A. Anshuman, T.I. Eldho and L.G. Singh

𝐶 (𝑥, 𝑦) =

𝑛 ∑ 𝑖=1

𝜙𝑖 (𝑥, 𝑦)𝑎𝑖 +

𝑚 ∑ 𝑗=1

𝑝𝑗 (𝑥, 𝑦)𝑏𝑗

Engineering Analysis with Boundary Elements 104 (2019) 8–25

This formulation is used in this study. By substituting Eqs. (16) and (17) into Eq. (8) we get following equation for node (xr , yr ). ( [ )] 𝑅𝜙𝑇 𝜕 2 𝜙𝑇 𝜕 2 𝜙𝑇 𝜕 𝜙𝑇 𝜕 𝜙𝑇 − 𝜃 𝐷𝑥𝑥𝑟 + 𝐷 𝑦𝑦𝑟 − 𝑉𝑥𝑟 − 𝑉𝑦𝑟 − 𝜆𝜙𝑇 C𝑠 𝑡+Δ𝑡 2 2 Δ𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 [ ] 2 𝜙𝑇 2 𝜙𝑇 𝑇 𝜕 𝜕 𝜕 𝜙 𝜕 𝜙𝑇 = 𝑅𝐶𝑟 𝑡 + (1 − 𝜃) 𝐷𝑥𝑥𝑟 + 𝐷 𝑦𝑦𝑟 − 𝑉𝑥𝑟 − 𝑉𝑦𝑟 − 𝜆𝜙𝑇 𝐶𝑠 𝑡 2 2 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 )( ) 𝑞𝑤 𝐶 ′ ( (18) + 𝛿 𝑥 − 𝑥𝑖 𝑦 − 𝑦𝑖 𝑛𝑒

(9)

The coefficients a and b in Eq. (9) are calculated by forcing C(x,y) to pass through the nodal values at all n points in the local support domain. This yields n number of equations with n + m variables (unknowns) as follows. 𝑛 𝑚 ( ) ∑ ( ) ∑ 𝐶 𝑥𝑘 , 𝑦𝑘 = 𝑎𝑖 𝜙𝑖 𝑥𝑘 , 𝑦𝑘 + 𝑝𝑗 𝑏𝑗 , 𝑘 = 1, … , 𝑛 𝑖=1

(10)

𝑗=1

Using the nodal indices for the ith node at the left-hand side of Eq. (18), it is modified for each node as [ )] ( 𝑅𝜙𝑇 𝜕 2 𝜙𝑇 𝜕 2 𝜙𝑇 𝜕 𝜙𝑇 𝜕 𝜙𝑇 − 𝜃 𝐷𝑥𝑥𝑖 + 𝐷 𝑦𝑦𝑖 − 𝑉𝑥𝑖 − 𝑉𝑦𝑖 − 𝜆𝜙𝑇 𝐶𝑖𝑡+Δ𝑡 Δ𝑡 𝜕𝑥 𝜕𝑦 𝜕 𝑥2 𝜕 𝑦2 [ ] 𝜕 2 𝜙𝑇 𝜕 2 𝜙𝑇 𝜕 𝜙𝑇 𝜕 𝜙𝑇 = 𝑅𝐶𝑟 𝑡 + (1 − 𝜃) 𝐷𝑥𝑥𝑟 + 𝐷 𝑦𝑦𝑟 − 𝑉𝑥𝑟 − 𝑉𝑦𝑟 − 𝜆𝜙𝑇 𝐶𝑠𝑡 𝜕𝑥 𝜕𝑦 𝜕 𝑥2 𝜕 𝑦2 )( ) 𝑞𝑤 𝐶 ′ ( + 𝛿 𝑥 − 𝑥𝑖 𝑦 − 𝑦𝑖 (19) 𝑛𝑒

Additional m equations are generated by employing the following m constraint conditions: 𝑚 ∑ 𝑖=1

( ) 𝑝𝑗 𝑥𝑖 𝑎𝑖 = 0, 𝑗 = 1, … , 𝑚

(11)

In this study, multi-quadrics radial shape functions (MQ-RBF) as defined in Eq. (12) are used. 𝜙𝑖 (𝑥, 𝑦) = (𝑟𝑖 2 + 𝑐𝑠2 )𝑞

(12)

Here q and cs are the shape parameters of the MQ-RBF. The shape parameter q is important for the accuracy of the interpolation. The standard multi-quadrics RBF has a q value of 0.5 [20,21]. The shape parameter cs is often expressed as 𝛼 × dc [10]. Here 𝛼 is factor of dc which is the average spacing between the nodes. The method for calculating dc for regularly and irregularly distributed nodes is given in Liu and Gu [10]. The radial distance ri is calculated using the point of consideration (xi , yi ) and other nodes in its local support domain. In case of 2D problems, the radial distance is computed as √ ( )2 ( )2 𝑟𝑖 = (13) 𝑥 − 𝑥𝑖 + 𝑦 − 𝑦𝑖

The above system of equations is for the internal nodes in the domain Ω. For the nodes on the Dirichlet boundary 𝜕Ω1 , the discretization is implemented as ( ) 𝜙𝑇 𝐶𝑠 = 𝐶0 𝑥𝑖 , 𝑦𝑖 , 𝑡 (20) Where C0 is the specified concentration at time t for location (xi ,yi ). Neumann boundary conditions are known to cause instability in the solution when using strong form methods [10]. Various approaches such as direct collocation, Hermite type collocation, fictitious point method etc. to deal with such problem are presented in Liu ad Gu [10]. In this study, a direct collocation approach is implemented for the nodes on the Neumann type boundary (𝜕Ω2 ).

The details for solving the system of Eqs. (10) and (11) are explained in Liu and Gu [10]. In this study polynomial basis functions are not used as accurate results have been reported in previous studies by using MQRBFs only [17,18,27]. Upon solving, the shape functions at n nodes are derived as { } 𝜙 = 𝜙1 , 𝜙2 , … … … … 𝜙𝑛 (14)

𝜕 𝜙𝑇 𝐶 = 𝑓𝑖 (𝑥, 𝑦, 𝑡) (21) 𝜕𝑛 𝑖 Where, fi is the known solute flux at location (x, y) and node i in n direction during time t. For impervious boundaries, the fluxes are set to zero. Singular sources e.g. a pump, are known to cause instability [29] which can be dealt with a Gaussian regularization technique [30]. Combining Eqs. (19), (20) and (21), for N boundary and internal nodes, the following system of equations will be obtained:

The variable C at each node is approximated by using the 𝜙 values for nodes within its support domain. For n nodes in the support domain s, the variable C at location (x,y) is computed as the following equation: 𝐶 (𝑥, 𝑦) = 𝜙𝑇 𝐶𝑠 =

𝑛 ∑ 𝑖=1

𝜙𝑖 𝐶𝑖

+Δ𝑡 = [𝐹 ]𝑁×𝑁 {𝐶 }𝑡𝑁×1 + {𝑄}𝑁×1 [𝐾 ]𝑁×𝑁 {𝐶 }𝑡𝑁×1

The system matrix formed due to the use of local support domains is sparse as the values of the shape function and its directional derivatives are considered zero outside the support domain of a given node. The sparsity of the matrix depends on the number of nodes in the support domain.

(15)

Similarly, the spatial derivative terms are also computed for all nodes using the shape function derivatives at their respective support domains [10].

4. Model development and verification

3.2. Discretization of governing equations using RPCM

4.1. Model development

For transport problems, the scalar variable considered in Eq. (9) is concentration of the dissolved solute. The governing Eq. (8) is discretized by collocation at the nodes and using the approximation of concentration from Eq. (9). The time discretization at any given node ‘i’ is implemented by the following forward difference scheme [2]: 𝐶 𝑡+Δ𝑡 − 𝐶𝑖𝑡 𝜕𝐶 = 𝑖 𝜕𝑡 Δ𝑡

Based on the formulation discussed in the previous section, ReactiveRPCM transport models named RT-RPCM-1D and RT-RPCM-2D are developed, respectively for one and two dimensions and implemented numerically in Python using Numpy package [31]. The important steps in this implementation are listed in the following steps and shown in the flowchart in Fig. 1.

(16)

The concentration C in the governing Eq. (8) is replaced by the semiimplicit formulation as follows. 𝐶 = 𝜃 𝐶 𝑡+Δ𝑡 + (1 − 𝜃)𝐶 𝑡

(22)

Step 1. In this step, the parameters of the model are defined. The variables such as velocity, dispersion coefficient, the rate of decay, retardation factor and initial concentrations are input into the model. Step 2. Domain discretization: The domain of the model is represented using uniformly/non-uniformly distributed nodes. The radius of the circular support domain is defined in terms of multiples of nodal distance dc .

(17)

Here 𝜃 is a weight coefficient and 𝜃 ∈ [0, 1] [25]. The numerical scheme is unconditionally stable if the value of 𝜃 is higher than 0.5. When the value of 𝜃 is 0.5, the formulation is known as the Crank– Nicolson scheme which is a fairly stable and widely used method [4]. 10

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 1. Flowchart for model development.

Step 3. The shape functions and their derivatives for each node are then evaluated using nodes within its local support domain. The shape function values for all locations are calculated using matrix G (n × n) and vector R (n × 1) as in Eq. (23), where n is the number of nodes in the support domain. Each row of G matrix is computed for a concerned node by estimating radial basis term (see Eq. (12)) with respect to other nodes in its support domain. The diagonal terms of the G matrix constitute R vector as shown in Eq. (23). The terms 𝜕 𝜙/𝜕 x, 𝜕 𝜙/𝜕 y, 𝜕 2 𝜙/𝜕 x2 and 𝜕 2 𝜙/𝜕 y2 are calculated easily by replacing the R vector with its derivatives keeping the G matrix unaltered. ⎡𝑔11 [𝐺]{𝜙} = {𝑅} ⇒ ⎢ ⋮ ⎢ ⎣𝑔 1 𝑛

⋯ ⋱ ⋯

𝑔𝑛1 ⎤⎧𝜙1 ⎫ ⎧𝑔11 ⎫ ⎪ ⎪ ⎪ ⎪ ⋮ ⎥⎨ ⋮ ⎬ = ⎨ ⋮ ⎬ ⎥⎪ ⎪ ⎪ ⎪ 𝑔𝑛𝑛 ⎦⎩𝜙𝑛 ⎭ ⎩𝑔𝑛𝑛 ⎭

tion values are then used as the initial values to compute matrices K and F for the next time step. Step 6. The steps 4 to 5 above are repeated for all the time steps. 4.2. Model verification 4.2.1. 1D transport for continuous source 4.2.1a. Single species decay. The governing equation for 1D advective– dispersive transport with single species decay is,

(23)

𝜕𝐶 𝜕2 𝐶 𝜕𝑐 =𝐷 −𝑣 − 𝜆𝐶 𝜕𝑡 𝜕𝑥 𝜕 𝑥2

(24)

The initial condition used is 𝐶 (𝑥, 0) = 𝐶𝑖

(25)

The boundary conditions used are

Step 4. The governing equations and boundary conditions are then discretized at each node using point collocation method. The resulting discretized nodal equations are assembled into a global stiffness matrix K and load diagonal matrix F (see Eq. (22)). Step 5. The concentration values for the current time step are then computed by solving above system of equations. These concentra-

𝐶 (0, 𝑡) = 𝐶0 𝜕𝐶 (∞, 𝑡) = 0 𝜕𝑥

(26)

The analytical solution for simultaneous linear equilibrium adsorption and first-order decay is given in [32]. 11

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

simulation closely matches the analytical solution which demonstrates the accuracy of the proposed RT-RPCM-1D model for first-order decay. 4.2.1b. Linear equilibrium adsorption. The governing equation for 1D advective–dispersive transport with linear equilibrium adsorption is, 𝑅𝜕𝐶 𝜕2 𝐶 𝜕𝑐 =𝐷 −𝑣 𝜕𝑡 𝜕𝑥 𝜕 𝑥2

(28)

The initial and boundary conditions for this case are kept same as in the previous case. Consider the same domain discussed in the previous subsection (Fig. 2(a)). In this case, the contaminant is taken as Hexachloroethane (HCE). It is produced as a by-product of chlorination processes in the industry. It affects the central nervous system as depressant [34]. The distribution coefficient (kd ) for HCE is 0.31 cm3 /g [1]. The porosity and bulk density of the material is assumed to be 0.3 and 1.59 g/cm3 , respectively. The retardation factor (R) is calculated to be 2.643 as in Eq. (8). Rest of the parameters are kept same as in Section 4.2.1a. The concentration profile obtained after a simulation period of 5 years and 10 years are shown in Fig. 4. The model simulations show good agreement with the analytical solution.

Fig. 2. Representation of (a) study domain (b) local support domain for case 4.2.1.

𝐶 (𝑥, 𝑡) { ( [ ] [ ]} ) ( ) 𝐶 (𝑣 − 𝑢)𝑥 (𝑣 + 𝑢)𝑥 𝑅𝑥 − 𝑢𝑡 𝑅𝑥 + 𝑢𝑡 = 𝑜 exp 𝑒𝑟𝑓 𝑐 √ + exp 𝑒𝑟𝑓 𝑐 √ 2 𝐷𝑥 𝐷𝑥 2 𝐷𝑥 𝑅𝑡 2 𝐷𝑥 𝑅𝑡 [ ] [ ]} { ( ) ( ) 𝐶 𝜆𝑡 𝑅𝑥 − 𝑣𝑡 𝑣𝑥 𝑅𝑥 + 𝑣𝑡 − 𝑖 exp − 𝑒𝑟𝑓 𝑐 √ + exp 𝑒𝑟𝑓 𝑐 √ 2 𝑅 𝐷𝑥 2 𝐷𝑥 𝑅𝑡 2 𝐷𝑥 𝑅𝑡 ( ) 𝜆𝑡 + 𝐶𝑖 exp − (27) 𝑅

4.2.2. 1D transport for injection type source In this example, a radioactive isotope of Iodine (131 I) is considered as the contaminant. It is widely used in the field of medicine and industry as tracer [35]. It is a fast decaying solute having a half-life of 8.02 days. The decay rate is computed using Eq. (4). A 200 m length domain is taken for the study and the source of the contaminant is located at the left boundary. The boundary conditions used are

where R is the retardation factor, v is the velocity, Dx is the dispersion √ 4𝜆𝐷 coefficient and 𝑢 = 𝑣 1 + 𝑣2 𝑥 . This test case represents a one-dimensional domain of 700 m length (Fig. 2(a)) which is initially in pristine condition (Ci = 0). The boundaries are far enough for the second boundary condition to be satisfied. A continuous source of Tritium (3 H2 ) with concentration C0 = 50 ppm is assumed to be present at the left boundary, while the right end boundary is impervious. Tritium is a radioactive isotope of hydrogen (H) which is used in radioactive tracing [33]. It has a half-life of 12.32 years and its decay rate (𝜆) can be calculated by Eq. (4) to be 0.000154 day−1 . The velocity v and dispersion coefficient D are 0.2 m/d and 0.2 m2 /d, respectively. 101 nodes are used to represent the domain with a regular nodal distance of 7 m as shown in Fig. 2(b) and radius of circular local support domain is kept as 4 × dc . Shape parameters q and cs are taken as 0.98 and 14 × dc respectively. The model is run for a simulation period of 15 years with a time-step of 10 days. The concentration profile after 5 and 10 years are plotted in Fig. 3. It is observed that the RPCM model

{ 𝐶 (0, 𝑡) =

10 ppm, 0,

𝑖𝑓 0 ≤ 𝑡 ≤ 10 days 𝑖𝑓 𝑡 > 10 days

𝐶 (200m, 𝑡) = 0

(29)

The initial condition used is 𝐶 (𝑥, 0) = 0

(30)

The flow velocity and dispersion coefficients are taken to be 0.2 m/day and 0.35 m2 /day, respectively. The domain is represented using 201 nodes with a support domain radius of 4 × dc . The shape parameters q and cs are the same as in case 4.2.1a. The model is run for a

Fig. 3. Concentration profiles at different periods (case 4.2.1a). 12

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 4. Concentration profiles at different periods (case 4.2.1b).

Fig. 5. Concentration profiles at different time periods (case 4.2.2).

simulation period of 1 year with a time step size of 1 day. For comparison of model accuracy, a semi-analytical solution using Anamodel tool [36] is used. The concentration profiles plotted in Fig. 5 show that the semi-analytical solution and RT-RPCM-1D simulation are in good agreement. These results show the applicability of the model for an injection type source and a rapidly decaying solute.

days. Remaining boundaries are constant concentration boundaries with C = 0 (Fig. 6). The initial condition and the source term are given by 𝑐 (𝑥, 𝑦, 0) = 0

(31)

𝑐 (0, 10 ≤ 𝑦 ≤ 40, 0 ≤ 𝑡 ≤ 30) = 5 ppm

(32)

The problem is solved using RPCM models having uniform and nonuniform nodal distributions (see Fig. 7) named as RT-RPCM-2D and RTRPCM-2D-NU respectively. In both models the domain is represented using 1681 nodes. In RT-RPCM-2D, uniform nodal distance (dc ) is 1.25 m. In RT-RPCM-2D-NU nodal distance is computed using an estimated area of local support domain and number of nodes [10]. The radius of local support domain for both models is 4 × dc . Similarly, MQ-RBF parameters q and cs are taken as 0.98 and 40 × dc for both models. The models are run for 100 days with a time step size of 1 day. In order to compare the results obtained by the models, Mean absolute error (MAE) and root mean square (RMSE) error are used as error indicators/metrics and are

4.2.3. 2D transport for injection type source The RPCM models developed in this study is verified against a semianalytical solution [36]. We consider 57 Co as the contaminant here. It is used in heavy metal industries [37] and also in medical applications [38]. Exposure to 57 Co for a prolonged period of time can be carcinogenic. An aquifer of dimension 50 m × 50 m is considered having a groundwater velocity field in the x-direction as 0.4 m/day. The dispersion coefficients are Dxx = 0.1 m2 /day and Dyy = 0.1 m2 /day. The contaminant taken is 57 Co which is injected at the left boundary for 30 13

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 6. Representation of study domain for (case 4.2.3).

Fig. 7. Nodal arrangements used for RT-RPCM-2D and RT-RPCM-2D-NU.

14

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 8. Breakthrough curves at points p1, p2, p3 and p4 for case 4.2.3. Table 1 RMSE and MAE values with varying nodal arrangements for case 4.2.3. Observation point

p1 p2 p3 p4

RMSE

MAE

Uniform nodes

Non uniform nodes

Uniform nodes

Non uniform nodes

0.073 0.037 0.012 0.005

0.075 0.040 0.010 0.006

0.041 0.024 0.005 0.002

0.041 0.026 0.004 0.002

defined as below: 𝑁 ∑

Mean absolute error, 𝑀𝐴𝐸 =

| 𝑒𝑥𝑎𝑐𝑡 | − 𝐶𝑖𝑛𝑢𝑚 | |𝐶 𝑖 |

𝑖=1 |

Root mean square error, 𝑅𝑀𝑆𝐸 =

𝑁

Table 1). However, the magnitudes of errors are very close. The results show that the RPCM models accurately predict the concentration for the concerned problem. In order to test the efficiency of the models further, more complicated problems are solved in the case studies.

(33)

𝑁 √ ∑ (𝐶𝑖𝑒𝑥𝑎𝑐𝑡 − 𝐶𝑖𝑛𝑢𝑚 ) 2

𝑖=1

𝑁

5. Sensitivity studies Parameter sensitivity studies are required to establish the optimal values of shape parameters, local support domain size etc. [17,18]. MQRBF produces accurate solutions for both fluid flow and solid problems when the value of shape parameter q is 0.98 or, 1.03 [10]. The value of shape parameter q is fixed at 0.98 for this study. In this section, we test the sensitivity of the proposed model with respect to shape parameter cs , local support domain radius and nodal distance. MAE and RMSE are used are error indicators in order to check the model accuracy.

(34)

where Ci exact is the semi-analytical solution and Ci num is the model simulation at node i. Concentration profiles at four randomly chosen locations are plotted in Fig. 8. It is observed that the simulations obtained by RT-RPCM-2D and RT-RPCM-2D-NU show little difference. RT-RPCM-2D model gives better results compared to RT-RPCM-NU at all points except p3 (see 15

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 9. Concentration profile after 5 years (case 4.2.1a) with varying cs . Table 2 Errors in RT-RPCM-1D model simulation (t = 5 years, case 4.2.1a) and the condition number of stiffness matrix for different cs values.

MAE RMSE Condition number

cs = 0.01 × dc

cs = 0.1 × dc

cs = 5 × dc

cs = 10 × dc

cs = 46 × dc

cs = 48 × dc

cs = 100 × dc

0.52 0.08 11.92

0.29 0.05 10.63

0.12 0.02 8.41

0.04 0.01 6.69

0.05 0.01 11.72

– – 86.34

– – 640.48

5.1. Sensitivity with respect to cs

tration gradient is the reason for this. Further, the error values decrease with respect to time. They become almost negligible beyond t > 4000 days when then the concentration finally reaches the far boundary. From Fig. 11 it is observed for support domain size 2 × dc the model shows instability where concentration changes sharply and upon increasing the support domain radius the solution is stabilized. The results are almost identical for support domain size 4 × dc and 7 × dc . The times taken for simulations by varying support domain radius are presented in Table 3. The ratio between the time taken by different support domain radii and that of support domain radius 2 × dc are also presented as relative time. The computation is done on a PC having specifications Core i7-6770, 3.40 GHz CPU and 16 GB RAM. From the Table 3, we observe that the time taken for model simulation increases with local support domain size which is expected as the size of the G matrix increases (see Eq. (23)). However, the increase is much higher in the case of the 2D model. It is due to more number of nodes are used in the 2D model for the same support domain size (see Figs. 2(b) and 6(b)) which increases the size of the G matrix even more. It is desirable for a model to produce accurate results using minimum computational effort. Considering the error values and relative time of computation from this test, it can be inferred that model performance is optimum when the support domain radius is 4 × dc.

This parameter of the MQ-RBF function must be properly chosen to avoid ill-conditioning of the system matrix and thus to provide stable solutions [10]. Ill-conditioning is indicated by high values of the condition number of the matrix. To examine its effect on model simulation, the simulations are run repeatedly with different cs values. The analytical solution from Section 4.2.1a is considered here for the tests keeping rest of the variables unchanged. Concentration profile after 5 years is used for comparison to track the concentration front. From Fig. 9 it is observed that the accuracy of model simulation declines when values of cs decreases to 0.1. Results match closely to the analytical solution when cs = 10 × dc and 46 × dc . However, when the parameter is increased to 48 × dc or more, the model performance becomes erratic and unreliable. It is due to the high value of condition number which denotes that the stiffness matrix is ill-conditioned. The error indicators in Table 2 show that model simulations are relatively accurate when cs = 10 × dc to 46 × dc . 5.2. Sensitivity with respect to local support domain radius The number of neighbouring nodes considered for the approximation of state variable for each node is determined by support domain radius. It is varied from 2 × dc to 7 × dc for problem discussed in Section 4.2.1a using 701 nodes. The error indicators are calculated for all time steps and plotted in Fig. 10. The results indicate that the magnitudes of errors are highest when t = 190 days. Presence of steep concen-

5.3. Sensitivity with respect to nodal density The size of the globally assembled stiffness matrix has the biggest impact on the computation time. In order to reduce the computational 16

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 10. MAE and RSME variation with respect to time of simulation case 4.2.1a.

plotted in Fig. 12. It shows that model simulations are in good agreement with analytical solutions when Pe = 10. Instability is observed near the steep concentration gradient region when Pe = 35 which increases further when Pe = 50. This instability can be alleviated by a number of methods as suggested by Liu and Gu [10]. Decreasing the Pe value by decreasing the nodal distance is one of them which is shown in Fig. 12 using 71 nodes (i.e. Δx = 10 m). Other methods include increasing the support domain size near boundary nodes and use of total or, adaptive upwind support domains [10]. 6. Case studies In this section, two case studies with reactive transport in twodimension for a regular rectangular and an irregular shaped domain are presented. The RPCM models are developed using uniform (RT-RPCM2D) and non uniform (RT-RPCM-2D-NU) nodal distributions. The model simulations are compared with FEM based models which are developed in COMSOL Multiphysics [39]. In both case studies, the porous medium is considered to be anisotropic. In order to compare model accuracy, Normalized MAE (NMAE) values are computed which are given by the following equation: ∑ ⎡ 𝑁 ||𝐶 𝑅𝑃 𝐶𝑀 − 𝐶 𝐹 𝐸𝑀 || ⎤ 𝑖=1 | 𝑖 𝑖 1 |⎥ ⎢ 𝑁 𝑀 𝐴𝐸 = (35) ( ) 𝐹 𝐸𝑀 ⎢ ⎥ 𝑁 max 𝐶1𝐹 𝐸𝑀 , … , 𝐶𝑁 ⎣ ⎦

Fig. 11. Comparison of concentration profiles at t = 190 days for different support domain radii.

effort, it is desirable to use less number of nodes i.e. increase the nodal distance. However, this increases the Peclet number (𝑃 𝑒 = 𝑉 𝐷Δ𝑥 ) and may cause instability in the model simulation. Here Δx is the nodal distance, V is the average velocity of flow and D is the dispersion coefficient. To study the sensitivity of the model with respect to nodal density, the nodal distance is gradually increased for the problem (in Section 4.2.1a) keeping the other parameters constant. The results are

6.1. Case study 1 Here a rectangular homogenous anisotropic aquifer of dimension 200 m × 120 m is considered (see Fig. 13). The source of pollution (C0 = 10 ppm) is located between 0–50 m in the left boundary. The values of velocities and dispersion coefficients are vx = 1 m/d, vy = 0.2 m/d, 17

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Table 3 Computational time variation with support domain radius. Support Domain Radius

2 × dc 3 × dc 4 × dc 5 × dc 6 × dc 7 × dc

1D problem 701 nodes

2D problem 1681 nodes

547 simulations (Case 4.2.1a)

100 simulations (Case 4.2.3)

Simulation time (s)

Relative time

Total time (s)

Relative time

14.51 17.89 22.43 26.53 29.64 33.42

1.00 1.23 1.55 1.83 2.04 2.30

11.96 25.39 40.88 102.37 126.91 133.86

1.00 2.12 3.42 8.56 10.61 11.19

Fig. 12. Concentration profile after 5 years (case 4.2.1a) for different nodal densities.

Dxx = 1 m2 /d and Dyy = 0.5 m2 /d. In first case, the solute is assumed to undergo linear adsorption while in the second simultaneous linear adsorption and decay is considered. The retardation factor (R) for linear adsorption is 3. The decay rate for first order decay is 0.05 d−1 which corresponds to half-life of approximately 13.86 days (see Eq. (4)). Two observation wells OB1 and OB2 are located at (20 m, 20 m) and (50 m, 50 m) respectively to measure the concentration profiles with respect to time. Both cases are solved with 1025 uniformly (RT-RPCM-2D) and non uniformly distributed nodes (RT-RPCM-2D-NU). The nodal arrangements are given in Fig. 14. The average nodal distance (dc ) for uniformly distributed nodes is 5 m in both x- and y- directions. The radius of local support domain is 4 × dc for both types of nodal arrangements. The MQRBF shape parameters q and cs are fixed as 0.98 and 5 × dc , respectively. The FEM model for the case consists of 2328 nodes and 2204 triangular elements. The average nodal distance in FEM model is 5.2 m. The Pe for both the cases is approximately 5. The simulation time for all the models is 150 days with time step of 5 days. The concentration profiles at OB1 and OB2 are presented in Fig. 15. The NMAE values are computed with respect to FEM results and presented in Table 4. It is observed that the NMAE values at the

Fig. 13. Representation of the study domain for case-study 1.

18

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 14. Nodal arrangements for case study 1: uniform nodal distribution (left) and non uniform nodal distribution (right).

Table 4 NMAE values for RT-RPCM-2D and RT-RPCM-2D-NU with respect to FEM for case study 1. Observation well

OB1 OB2

Reactions Nodal distribution

Linear adsorption

Linear adsorption and first order decay

Uniform Non uniform Uniform Non uniform

0.033 0.034 0.023 0.006

0.036 0.037 0.024 0.029

Table 5 Percentage difference of RT-RPCM-2D-NU results with respect to that of RT-RPCM-2D for case study 1. Reaction type

Adsorption Adsorption and decay

Observation well

OB1 OB2 OB1 OB2

t = 50 days

t = 150 days

RT-RPCM-2D

RT-RPCM-2D-NU

% Difference

RT-RPCM-2D

RT-RPCM-2D-NU

% Difference

2.728 0.000 1.426 0.000

2.698 0.000 1.408 0.000

1.01 0 1.26 0

9.989 4.803 3.857 0.604

9.990 5.371 3.853 0.699

0.10 11.83 0.10 15.73

observation wells nearer to the source OB1 are higher compared to that of OB2 for both reaction types. The concentration profile obtained by FEM for the linear adsorption case is fairly stable as seen in Fig. 15(a) and (b). However, in case of simultaneous adsorption and decay FEM results shows oscillations at both observation wells. The RPCM models give almost similar results at OB1 for both reaction types. The RT-RPCM-2D-NU model gives higher concentration values at OB2 for both reaction types when compared to RT-RPCM-2D (see Fig. 15(b) and (d)). At t = 50 days the percentage difference is almost negligible at OB1 and OB2 for both reaction types (see Table 5). At OB2 the results obtained by RT-RPCM-2D-NU model are comparable to that of FEM as seen in Fig. 15. It is to be noticed that OB1 is located near to the boundary. Since boundary nodes are kept same for both type of nodal distribution, most of the nodes in the support domain of OB1 are same for RT-RPCM-2D and RT-RPCM-2D-NU. On the other hand, OB2 is located far from the boundary. The locations of the nodes in the support domain of OB2 are different for different nodal arrangements. That is why the difference in RPCM model results is higher at OB2 compared to OB1. Nonetheless, the concentration profiles obtained by the proposed models don’t exhibit oscillations unlike FEM.

6.2. Case study 2 In this case, an aquifer having irregular boundary is considered (see Fig. 16). The study domain is more complicated than that of case study 1 and similar to the aquifers encountered in field case studies. The medium of aquifer is assumed to be anisotropic with transport parameters vx = 1.25 m/d, vy = −0.2 m/d, Dxx = 5 m2 /d and Dyy = 5 m2 /d. An aerial source of dimensions 50 m × 50 m injects contaminants at the rate of 1 ppm/d. The concentrations of contaminants at the entire boundary are maintained at 0 ppm. To measure the variation of concentration, two observation wells are located at OB1 (900 m, 500 m) and OB2 (1100 m, 400 m). Two different cases are considered in which the contaminant is assumed to undergo linear adsorption and simultaneous adsorption and decay, respectively. For the linear adsorption case, the retardation factor (R) is 1.8. The decay rate (𝜆) for the simultaneous adsorption and decay case is 0.0086 day−1 which corresponds to half-life of 80.58 days (see Eq. (4)). Two RPCM based models are developed with uniformly (RT-RPCM2D) and non-uniformly (RT-RPCM-2D-NU) distributed nodes as shown in Fig. 17. For the uniform nodal distribution, the nodal distance (dc ) 19

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 15. Comparison of concentration profiles obtained at observation wells OB1 and OB2 for linear adsorption case (a and b) and for simultaneous linear adsorption and decay case (c and d) for case study 1.

is fixed as 25 m. The MQ-RBF shape parameters q and cs are taken as 0.98 and 5 × dc , respectively. For both types of nodal arrangements, the support domain radius is fixed at 4 times the nodal distance (dc ). The models are run for 2 years of simulation period with time step of 5 days. The model simulations are compared with FEM model with 5142 nodes and triangular 4868 elements. The Pe for the RPCM models and FEM model are 6.25 and 5.94, respectively. The concentration contours after 2 years of simulation period for RT-RPCM-2D and RT-RPCM-2D-NU models are presented in Fig. 18. In order to analyze the variation of concentration with respect to time, the concentration profiles at observation wells OB1 and OB2 are compared in Fig. 19 for all models. It is observed that, the FEM model simulations are observed to be stable for the linear adsorption case (see Fig. 19(a) and (b)). Nevertheless, noticeable oscillations in the concentration pro-

file are present for simultaneous adsorption and decay case at OB1 and OB2 (see Fig. 19(c) and (d)). These oscillations are more prominent at OB1 which is nearer to the source. The concentration profiles at OB1 and OB2 obtained by RPCM models are noticeably different. The NMAE values computed with respect to FEM show higher for RT-RPCM-2D compared to RT-RPCM-2D-NU for most of the cases (Table 6). The RT-RPCM-2D model over predicts compared to RT-RPCM-2D-NU and FEM the concentration at OB1 and OB2 for both reaction types. The simulations obtained by RT-RPCM-2DNU are more comparable to that of FEM. The percentage absolute difference between RPCM model results is observed to increase with time for most of the cases (see Table 7). Compared to case study 1, the difference in RPCM model results is more distinguishable as both observation points are located far from the boundary. Therefore, the locations of 20

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 16. Representation of the study domain for case-study 2.

Fig. 17. The nodal arrangement for uniformly distributed nodes (left) and non uniformly distributed nodes (right).

the nodes in the local support domain of OB1 and OB2 are different for RT-RPCM-2D and RT-RPCM-2D-NU models. In Section 4.2.3, the RPCM models give almost similar results in the regularly shaped homogenous and isotropic porous medium. Nonetheless, the model simulations with uniformly and non-uniformly distributed nodes vary when more complicated problems are solved including anisotropy and irregular boundaries as in the case studies. This variation is more pronounced when the observation points are located far from the boundary. However, the concentration profiles obtained by RPCM models are stable for the problems considered whereas FEM models show prominent oscillations are observed when simultaneous adsorption and first order decay are considered.

that the concentration profiles for linear adsorption case are stable and similar for both radial support domain sizes. However, oscillation near the concentration front is noticeable for simultaneous linear adsorption and decay case when radius of local support domain reduced to 2 × dc (see Fig. 20(b)). Hence, increasing the radius of local support domain improves the model simulation by eliminating the oscillations near the concentration front. The advantage of the proposed method is that this increase can be easily implemented while it is difficult to do so in FEM and FDM [10]. 7. Discussion As observed in this study, the proposed RT-RPCM 1D and 2D models are effective in simulating reactive transport in porous media. Collocation using local support domains solves the ill-conditioning problem associated with Kansa approach [20,21] by yielding a sparse matrix system. From the sensitivity results, we have found that a local support domain radius (size) of 4 × dc gives the optimal balance in terms of computational effort and accuracy. The accuracy and stability of model simulations are also affected by shape parameters q and cs of multiquadrics RBF (MQ-RBF). In this study, the value of parameter q is fixed to 0.98 as suggested by Liu and Gu [10]. The shape parameter cs was

6.3. Stability analysis for linear adsorption and first order decay case In the case studies presented in previous subsections, it is observed that the FEM models show oscillations upon introduction of first order decay term. Codina [40] showed that the instability caused by the first order decay in FEM solution is similar to that of the advective term. However, the RPCM model simulations are free from such instability. The concentration profile for OB1 from case study 1 is analyzed by decreasing the radius of local support domain to 2 × dc . It is observed 21

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 18. Contour plots at t = 730 days for concentration(units in ppm) for linear adsorption case: (a) RT-RPCM-2D (b) RT-RPCM-2D-NU and simultaneous linear adsorption and first order decay case: (c) RT-RPCM-2D (d) RT-RPCM-2D-NU.

Table 6 NMAE values for RT-RPCM-2D and RT-RPCM-2D-NU with respect to FEM for case study 2. Observation well

OB1 OB2

Nodal distribution

Uniform Non uniform Uniform Non uniform

Reactions Linear adsorption

Linear adsorption and first order decay

0.025 0.009 0.013 0.004

0.034 0.019 0.022 0.011

varied from 0.01 × dc to 100 × dc. It was noted that the condition number for the system matrix becomes very high when cs = 48 × dc or, more rendering unreliable simulation. Besides, error in the simulation was observed to increase gradually with a decrease in cs value when cs ≤ 5 × dc . Fixing it within 10 × dc to 46 × dc yielded good results for the problem considered. Further, the model sensitivity was tested with different nodal densities. When the nodal density is decreased, the Peclet number increases and the problem becomes advection dominant which destabilizes the solution. The standard grid based methods such as FEM and FDM is known to show oscillations when Pe > 2 [41,42]. Special treatments such as upwinding can be employed to simulate transport for higher Pe values [41,42]. However, the proposed model simulations were satisfactory even when Pe = 10 without adopting any such stabilization techniques. Two case studies are presented involving linear adsorption and simultaneous adsorption and first order decay. The problems are solved in

regular and irregular domains with uniform and non uniform nodal distributions. In the first case study, continuous source is considered at the boundary while areal injection type source is considered for the second case study. The nodal arrangement did not influence the results considerably for a regular boundary problem with homogenous and isotropic porous media (see Section 4.2.3). However, for more complicated problems considered in the case study, noticeable variation is observed for different nodal arrangement especially for observation points located far from boundary. The FEM models produce oscillating concentration profiles at the observation points in case of simultaneous linear adsorption and first order decay. These oscillations may be countered by implementing special treatments/ techniques such as the introduction of artificial diffusion, crosswind diffusion and upwinding. Apart from that, operator splitting is another technique that is used in grid-based models for handling reaction terms [43,44]. However, in the proposed models reactions are solved simultaneously without adopting any spe22

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 19. Comparison of concentration profiles obtained at observation wells OB1 and OB2 for linear adsorption case (a and b) and for simultaneous linear adsorption and first order decay case (c and d) for case study 2. Table 7 Percentage absolute difference in RT-RPCM-2D-NU with respect to RT-RPCM-2D for case study 2. Reaction type

Observation well

Adsorption

OB1 OB2 OB1 OB2

Adsorption and decay

t = 365 days

t = 730 days

RT-RPCM-2D

RT-RPCM-2D-NU

% Difference

RT-RPCM-2D

RT-RPCM-2D-NU

% Difference

0.099 0.000 0.017 0.000

0.100 0.000 0.019 0.000

1.01 0 11.76 0

12.447 1.341 1.032 0.052

11.727 1.186 0.976 0.046

6.01 11.56 5.43 11.54

decaying solute such as 131 I. Since the model is meshfree, it is flexible in adaptive analysis and representation of any irregularly shaped domain. The use of local support domain yields sparse system matrices which remove the ill-conditioning associated with global support domain. The proposed model can be used to simulate the transport of solutes undergoing first order decay and linear adsorption in groundwater.

cial treatments. It is also observed that the role of radius of local support domain is important for producing stable solutions for simultaneous adsorption and decay cases. Use of larger support domain eliminates the oscillations near the concentration front. However, care should be taken while enlarging the radius of local support domain as it increases the computational effort tremendously especially for two-dimensional problems (see Table 2).The model solutions are also stable for very fast 23

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

Fig. 20. Comparison of concentration profiles at OB1 obtained by RT-RPCM-2D for (a) linear adsorption case and (b) simultaneous adsorption and decay case.

8. Conclusions

References

In this paper, a meshfree RPCM model is developed for simulation of reactive solute transport in porous media which takes into account of first-order decay and linear adsorption. Shape function evaluation for function interpolation is performed using circular shaped local support domains. The use of local support domain yield sparse system matrix. Hence it is free from the ill-conditioning problem associated with using a globally supported domain. The one- and two-dimensional models are verified against analytical and semi-analytical solutions for transport problems involving continuous and injection type sources. Sensitivity analysis is performed with respect to model parameters in order to determine their optimum ranges. Case studies are presented to compare the model performance against FEM. In the case studies, two dimensional RPCM models with uniform and non uniform nodal distributions are compared with FEM based models for rectangular and irregular boundary problems. The proposed model produces oscillation-free solutions for the given high Peclet number problems without requiring any special stabilization techniques such as upwinding or, operator splitting. Besides that, since the method is truly meshfree, it reduces the computational cost required in meshing and is suitable for adaptive analysis. As demonstrated, the RPCM model can be effectively used for irregular boundary problems. Hence the proposed model can be a good alternative to traditional mesh-based methods for reactive transport problems in porous media.

[1] Batu V. Applied flow and solute transport modeling in aquifers: fundamental principles and analytical and numerical methods. CRC Press; 2005. [2] Bear J, Cheng AH. Modeling groundwater flow and contaminant transport. Springer Science & Business Media; 2010. [3] Means JL, Crerar DA, Duguid JO. Migration of radioactive wastes: radionuclide mobilization by complexing agents. Science 1978;200(4349):1477–81. [4] Rastogi AK. Numerical groundwater hydrology. India: Penram International Publication; 2007. [5] Essaid HI, Cozzarelli IM, Eganhouse RP, Herkelrath WN, Bekins BA, Delin GN. Inverse modeling of BTEX dissolution and biodegradation at the Bemidji, MN crude-oil spill site. J Contam Hydrol 2003;67(1–4):269–99. [6] Poulsen M, Lemon L, Barker JF. Dissolution of monoaromatic hydrocarbons into groundwater from gasoline-oxygenate mixtures. Environ Sci Technol 1992;26(12):2483–9. [7] Calvet R. Adsorption of organic chemicals in soils. Environ Health Persp 1989;83:145. [8] Gupta RC, Milatovic D. Insecticides. Biomarkers in toxicology. Cambridge, Massachusetts: Academic Press; 2014. p. 389–407. [9] Bolla KI, Cadet JL. Exogenous acquired metabolic disorders of the nervous system: toxins and illicit drugs. Textbook of clinical neurology. Third Ed. Elsevier Inc; 2007. [10] Liu GR, Gu YT. An introduction to meshfree methods and their programming. Springer Science & Business Media; 2005. [11] Hirsch C. Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics. Elsevier; 2007. [12] Patankar S. Numerical heat transfer and fluid flow. CRC Press; 1980. [13] Oñate E, Manzan M. Stabilization techniques for finite element analysis of convection-diffusion problems. Devel Heat Transf 2000;7:71–118. [14] Boddula S, Eldho TI. A moving least squares based meshless local Petrov–Galerkin method for the simulation of contaminant transport in porous media. Eng Anal Bound Elem 2017;78:8–19. [15] Herrera PA, Massabó M, Beckie RD. A meshless method to simulate solute transport in heterogeneous porous media. Adv Water Resour 2009;32(3):413–29. [16] Majumder P, Eldho TI. Vectorized simulation of groundwater flow and contaminant transport using analytic element method and random walk particle tracking. Hydrol Process 2017;31(5):1144–60. [17] Meenal M, Eldho TI. Simulation of groundwater flow in unconfined aquifer using meshfree point collocation method. Eng Anal Bound Elem 2011;35(4):700–7. [18] Meenal M, Eldho TI. Two-dimensional contaminant transport modeling using meshfree point collocation method (PCM). Eng Anal Bound Elem 2012;36(4):551–61. [19] Hardy RL. Multiquadric equations of topography and other irregular surfaces. J Geophys Resear 1971;76(8):1905–15. [20] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics—I surface approximations and partial derivative estimates. Comput Math Appl 1990;19(8–9):127–45. [21] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics—II solutions to parabolic, hyperbolic and elliptic partial differential equations. Comput Math Appl 1990;19(8–9):147–61.

Acknowledgments The work reported in this paper is supported by a project sanctioned by Board of Research in Nuclear Sciences (BRNS) titled “Modeling of reactive transport in groundwater using meshfree based numerical methods” (Project no. 16BRNS002). Authors are thankful to BRNS for their support. Authors would also like to thank the anonymous reviewers for their comments and thereby improving the quality of the manuscript to a great extent. Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.enganabound.2019.03.016. 24

A. Anshuman, T.I. Eldho and L.G. Singh

Engineering Analysis with Boundary Elements 104 (2019) 8–25

[22] Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations. Comput Math Appl 2000;39(7–8):123–38. [23] Chen CS, Ganesh M, Golberg MA, Cheng AD. Multilevel compact radial functions based computational schemes for some elliptic problems. Comput Math Appl 2002;43(3–5):359–78. [24] Šarler B, Vertnik R. Meshfree explicit local radial basis function collocation method for diffusion problems. Comput Math Appl 2006;51(8):1269–82. [25] Singh LG, Eldho TI, Kumar AV. Coupled groundwater flow and contaminant transport simulation in a confined aquifer using meshfree radial point collocation method (RPCM). Eng Anal Bound Elem 2016;66:20–33. [26] Singh LG, Eldho TI, Kumar AV. Identification of groundwater contamination sources using meshfree RPCM simulation and particle swarm optimization. Water Resour Manag 2018;32(4):1517–38. [27] Thomas A, Majumdar P, Eldho TI, Rastogi AK. Simulation optimization model for aquifer parameter estimation using coupled meshfree point collocation method and cat swarm optimization. Eng Anal Bound Elem 2018;91:60–72. [28] Kumar RP, Dodagoudar GR. Meshfree analysis of two‐dimensional contaminant transport through unsaturated porous media using EFGM. Int J Numer Methods Biomed Eng 2010;26(12):1797–816. [29] Hesthaven JS, Gottlieb S, Gottlieb D. Spectral methods for time-dependent problems. Cambridge University Press; 2007. [30] Jung JH. A note on the spectral collocation approximation of some differential equations with singular source terms. J Sci Comput 2009;39(1):49–66. [31] Oliphant TE. A guide to NumPy. USA: Trelgol Publishing; 2006. [32] Van Genuchten MT, Alves WJ. Analytical solutions of the one-dimensional convective-dispersive solute transport equation. Economic Research Service; Technical Bulletin, No 1661. Washington DC: US Department of Agriculture; 1982. p. 8–56. [33] Jacobs DG. Sources of tritium and its behavior upon release to the environment. AEC critical review series. Oak Ridge National Lab., Tenn; 1968.

[34] Ruder AM. Potential health effects of occupational chlorinated solvent exposure. Annal New York Acad Sci 2006;1076(1):207–27. [35] Hoefnagel CA, Voute PA, Marcuse HR. Radionuclide diagnosis and therapy of neural crest tumors using iodine-131 metaiodobenzylguanidine. J Nuclear Med Offic Publ Soc Nuclear Med 1987;28(3):308–14. [36] Goltz M, Huang J. Analytical modeling of solute transport in groundwater: Using models to understand the effect of natural processes on contaminant fate and transport. John Wiley & Sons; 2017. [37] Goldoni M, Catalani S, De Palma G, Manini P, Acampa O, Corradi M, Bergonzi R, Apostoli P, Mutti A. Exhaled breath condensate as a suitable matrix to assess lung dose and effects in workers exposed to cobalt and tungsten. Environ Health Persp 2004;112(13):1293. [38] Versluis A, Rasker JJ, Jurjens H, Woldring MG. Labelling of bleomycin with cobalt-57, indium-111, technetium-99m, mercury-197, lead-203, and copper-67. Nuklearmedizin. Nucl Med 1976;15(2):86–90. [39] Comsol AB. Comsol Multiphysics Reference Manual. Stockholm: COMSOL AB; 2007. [40] Codina R. Comparison of some finite element methods for solving the diffusion-convection-reaction equation. Comput Methods Appl Mech Eng 1998;156(1–4):185–210. [41] Spalding DB. A novel finite difference formulation for differential expressions involving both first and second derivatives. Int J Numer Methods Eng 1972;4(4):551–9. [42] Heinrich JC, Huyakorn PS, Zienkiewicz OC, Mitchell AR. An ‘upwind’finite element scheme for two‐dimensional convective transport equation. Int J Numer Methods Eng 1977;11(1):131–43. [43] Zienkiewicz OC, Taylor RL, Nithiarasu P. The finite element method for fluid dynamics, volume 3. Elsiever; 2005. [44] Zheng C, Wang PP. MT3DMS: a modular three-dimensional multispecies transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems; documentation and user’s guide. Alabama University; 1999.

25