Peridynamics modelling of coupled water flow and chemical transport in unsaturated porous media

Abstract A wide range of phenomena in natural or engineered systems emerge from strongly coupled hydraulic, chemical and mechanical processes involving a mix of clearly discrete and quasi-continuous mechanisms. Non-local formulations of these processes, e.g. based on Peridynamics, offer significant advantages compared to classical local formulations. Existing Peridynamics models for water flow and chemical transport are applicable only to saturated systems and use loosely coupling schemes, such as explicit time stepping approaches. This work advances the non-local approach by developing a bond-based formulation for coupled water flow and chemical transport in partially saturated porous media. An implicit solution is proposed for coupling the PD formulation of chemical transport with water flow formulation. Firstly, the proposed formulation is verified against results from finite element/finite difference transient solutions for 1-D and 2-D coupled problems. The agreement between results demonstrates the accuracy of the proposed methodology. Secondly, a series of case studies are presented to illustrate the model’s capability to capture discontinuities and heterogeneities, including stationary cracks, propagating cracks, and randomly distributed permeable and impermeable inclusions. The results show that the multi-physics Peridynamics-based formulations and computational model (Pyramid) provide clear advantages to classical local formulations for analyses of reactive transport in partially saturated porous media with physically realistic microstructures.


Introduction
Modelling approaches based on classic formulations of transport phenomena in porous material utilise the basic physics of fluid flow and chemical transport where the fluid and chemical fluxes are proportional to the gradient of a potential which are applied alongside the constitutive equations, describing the material's properties. In such formulations, the conservations of mass, momentum and energy are locally imposed via partial differential equations. Examples include Fick's law, Darcy's law and Fourier's law. For problems where processes are strongly coupled, additional complexity arises in the modelling approaches from the non-linearity and cross-coupling effects. Approximate solutions to coupled transport processes, described by the local theories, can be obtained by various numerical techniques (e.g. Finite Element Method (FEM), Finite Difference Method (FDM), and Boundary Element Method (BEM)).
The local formulations, however, are less efficient for analysing heterogeneous media or media with discontinuities (Jafarzadeh et al., 2019b). Examples include porous systems, whose fundamentally discrete structures span several length scales with associated discretely changing transport properties (e.g., permeability and diffusivity), and systems with discontinuous transport across phase boundaries (Carciopolo et al., 2020). Such underlying structures result in fast and slow transport processes existing in the same computational domain. It becomes extremely challenging to solve problems with heterogeneities and discontinuities without extra requirements such as multi-scale algorithms and adaptive multilevel schemes (Helmig et al., 2010). Moreover, in a natural heterogeneous porous medium, the moisture flux at a given location depends on the pressure distribution in a finite (not infinitesimal) size neighbourhood, which requires 'connections' in discretised domain (Delgoshaie et al., 2015;Jenny and Meyer, 2017). In the presence of discontinuities, such as cracks or shocks, additional requirements such as dependence of density on concentration, need to be imposed to represent the flow realistically (Koldoba and Koldoba, 2003). Without imposing these requirements, the solutions of the partial differential equations may not converge. Furthermore, the classical local formulations of chemical transport in porous media have limitations to capture important phenomena, such as anomalous transport observed in heterogeneous or discontinuous environments, where nonlocal effects are important (Liu et al., 2019(Liu et al., , 2017Zhao et al., 2018).
Nonlocal formulations, e.g. weekly nonlocal theories, nonlocal elasticity, nonlocal strain-softening models and Peridynamics, have been developed in recent decades to address problems with discontinuous systems, large deformations and heterogeneities (Bažant and Chang, 1984;Eringen and Edelen, 1972;Rogula, 1982;Silling, 2000). Among those, the Peridynamics (PD) theory is receiving increasing interest among the continuum mechanics community. The theory of peridynamics (Silling, 2000) presents a set of equations of continuum mechanics to overcome the limitations of the classical local theory in solving problems with discontinuities. A major development in this theoretical framework is that integral-differential equations replace the partial differential equations of the classical local theory. This results in an obtaining consistent formulation from mathematical point of view where discontinuities are developed in the porous media (e.g. fracture, fragmentation). The PD theory has successfully been applied in different areas, including fluid driven fractures (Oterkus et al., 2017;Ouchi et al., 2015), thermal fracture Wang and Zhou, 2019), membranes and fibers (Silling and Bobaru, 2005), dynamic fracture of nanofiber (Bobaru, 2007), pitting corrosion damage (Jafarzadeh et al., 2019a) and electro-migration (Gerstle et al., 2008).
The PD theory bypass the calculations of spatial derivatives, which can be mathematically undefined in the presence of sharp discontinuity in the field variables (such as a crack or a shock) . This feature allows for developing PD formulations for physical phenomena beyond the original applications to solids behaviour proposed by Silling (2000). For example, Bobaru and Duangpanya (2010), Bobaru and Duangpanya (2012) presented a PD model for thermal diffusion which adopted the bond-based version of PD. The generic PD formulations have also been used for studying mass diffusion (Du et al., 2012;Tian et al., 2017;Tian and Du, 2013). Katiyar et al. (2014) proposed a PD formulation for steady state flow of water in saturated porous media; addressing discontinuities and heterogeneities in the modelling of fluid flow. PD models for coupled flow and deformation have also been presented to study fracture initiation and propagation by a cumulative damage concept, which tracks broken mechanical bonds for material points Ouchi et al., 2015). A PD model for transient moisture flow in unsaturated porous media was also proposed by Jabakhanji and Mohtar (2015); which provides a solution to the classic Richard's equation of water flow in unsaturated porous media. The PD formulations presented for water flow and chemical transport problems are limited to those solving either the water flow (Jabakhanji and Mohtar, 2015;Katiyar et al., 2014) or the chemical transport (Gerstle et al., 2008) separately, or offering a weak coupling scheme to deal with the two processes only in saturated systems (Zhao et al., 2018).
The coupling methods used in multiphysics problems include the strong coupling approach (i.e., direct coupling), and the weak coupling approach (i.e., iterative coupling) (Zhang and Cen, 2015). A weakly coupled method has been introduced in PD to solve coupled problems such as thermo-mechanical (D'Antuono and Morandini, 2017;Wang and Zhou, 2019;Wang et al., 2018). In a weakly coupled method, different time stepping strategies are generally employed for different physics/processes. PD formulations of coupled multiphysics problems have been developed which are based on sequential numerical approaches to solve the coupled problem (Agwai, 2011;Madenci and Oterkus, 2014;Oterkus et al., 2014Oterkus et al., , 2017Ouchi et al., 2015). In these models the processes are sequentially addressed within a time step which can fail to satisfy convergence in the solutions. This is often tacked at the expense of adopting a small-time step .
To the authors' knowledge, there is no PD formulation which can describe an implicit coupled water flow and chemical transport in partially saturated porous media. Such coupled hydro-chemical formulation is an essential step to develop PD-based fully coupled hydrochemical-mechanical models for studying the chemo-mechanical problems such as erosion assisted degradation of materials or hydraulic-chemo-mechanical behaviour of fractured systems where coupling effects are inherently combined with modelling challenges associated with discontinuity and heterogeneity of the system. By using the non-local theory of peridynamics, transport in systems with spatially varying geometrical and physical properties can be solved without resort to special multi-scale or adaptive schemes. This paper presents the development of a PD formulation for coupled water flow and chemical transport in discontinuous and/or heterogeneous unsaturated porous media. Bond-based Peridynamics formulations for water flow and chemical transport, considering advection, dispersion and diffusion processes in partially saturated porous media, are derived and presented. The proposed coupling of water flow and chemical transport processes is implemented based on an implicit numerical scheme. The new coupled hydraulic and chemical Peridynamics model (Pyramid) is tested against several verification bench marks to ensure accurate derivation and implementation of the formulations. A series of case studies are presented of transient water flow and chemical transport in fractured and heterogeneous unsaturated porous media, by which the effects of fractures and randomly distributed permeable and impermeable regions are investigated by the new capabilities of the proposed model.

Bond-based peridynamics formulations for water flow and chemical transport
This section describes a bond-based PD formulation for the transport equations. The classical continuum equations for water flow and chemical transport are provided in Appendix A. The PD considers a body occupying a region (R) as a collection of interacting material points x. The points are identified with their coordinates: (x) in 1-D, (x, y) in 2-D, and (x, y, z) in 3-D problems. Specifically for transport problems, each point is assigned physical variables: volume and concentration, from where its mass is calculated; or volume and mass, from where its concentration is calculated. In the following description, we will assume that the mass is the primary variable for transport problems, instead of the alternative concentration. A material point x interacts with (is connected to) all materials points x ′ within a certain finite region H x , called the horizon of x (see Fig. 1). The radius of the horizon is denoted by δ. The neighbourhood of x is divided into a region, where the points have larger mass than x, denoted by H x+ , and a region, where the points have smaller mass than x, denoted by H x-. S denotes the equivalent mass at the surface of the horizon of x. The conductor of mass between two interacting material points x and x′ is called 'mass bond' or 'm-bond' (Bobaru andDuangpanya, 2010, 2012). The Peridynamics mass flux per unit volume along PD bond is assumed to be dependent on the distance between the material points x and x′ (Bobaru and Duangpanya, 2012). The pairwise mass exchange between material points x and x′ is denoted by f h .
The reactive chemical transport and water flow PD formulations will be presented in Sections 2.1 and 2.2. There are four steps for the construction of PD transport formulations.
Step I: we construct the peridynamic flux along the bond between material points x and x′ (Fig. 1). Step II: the conservation of mass law is used to develop the peridynamic transport formulations along a single m-bond between points x and x′.
Step III: We integrate peridynamic flux over the horizon of material point x to obtain the peridynamic transport formulations. It is noted that x ′ − x or ‖x ′ − x‖ is treated as a finite distance, which is different from the way that the mass conservation equation is obtained in the classical transport model where a limit of x ′ →x is taken (Bobaru and Duangpanya, 2010).
Step IV: We connect peridynamics' parameters to the measurable material properties to obtain the PD micro-diffusivity and micro-hydraulic conductivity by using different influence functions, which will be presented in Section 2.3. To clarify further these points, we provide detailed derivations of 1-D pure chemical diffusion PD formulations in Appendix B.

PD formulations for chemical transport
The vector of Peridynamics' chemical mass flux along the 'm bond' connecting points x and x ′ is given by Zhao et al. (2018): is the diffusivity of the 'm-bond' between the material points x and x ′ , C(x ′ , t) and C(x, t) are chemical concentrations at material points x ′ andx, respectively; ξ represents the distance vector between the two material points x and x ′ , i.e. ξ = x ′ − x (Silling, 2000), ξ ‖ξ‖ is the unit vector along the mass bond; and q(x, x ′ , t) is the average water flux between the two material points x and x ′ .
In previous studies where saturated systems were of interest (e.g. (Zhao et al., 2018)) the chemical advective flux induced by the water flux are assumed to be constant. However, considering a transient advection flux (Eq. A-2 in Appendix A) and variable unsaturated hydraulic conductivity (Eq. A-4 in Appendix A) is necessary for modelling the chemical advection in partially saturated system which creates further nonlinearity. In order to couple the transient advection induced by the water flow, an average water flux q(x, x ′ , t) is defined (Oterkus et al., 2017): q(x, t) and q(x ′ , t) are the liquid velocity vectors of the material points x and x ′ . We note that for 2-D problems, q(x, t) has components in x and y direction: q x (x, t) and q y (x, t). The velocity/advection vector can be obtained by solving the water transport equations, which will be given in Section 2.2. The conservation of mass law led to the peridynamic chemical transport in a single m-bond (x, x ′ ) in the direction ξ ‖ξ‖ as follow (Bobaru and Duangpanya, 2010): Combining Eqs. (1) and (3), the mass conservation equation for the bond along the two material points x and x ′ can be presented as: In order to obtain the mass conservation at the material point x due to a chemical flux from all the bonds attached to material point x within its horizon (H x ), we integrate Eq. (4) over the horizon H x of material pointx: ∫ The uniform grid (spacing Δx) is used for desretization in this study. The volume of x ′ is therefore represented by Δx 3 (Δx 2 for 2D, Δx for 1D). It is worth noting that we have implemented successfully non-uniform grids, based on Voronoi tesselations within peridynamics models. However, for the purposes of the present work, it is sufficient to use uniform grids, and we have used only those. The relationship between the concentration at point x and time t and the average concentration in all the m-bonds connected to material point x is given by Bobaru and Duangpanya (2010): where, V Hx is the horizon volume for the node centred at x. Combining Eqs. (5) and (6) leads to the following equation: The Peridynamics microscopic diffusivity and the Peridynamics microscopic liquid flux are defined respectively as: v The pairwise mass exchange f h is also defined by: By substituting Eqs. (8) and (9) into Eq. (7) and combining with Eq. (10), the bond-based Peridynamics formulation for the coupled water flow and chemical transport in partially porous media at the material point at xcan be expressed in the following form: The micro-advection term, v(⋅, ⋅, ⋅) in Eq. (11), is induced by the water flow due to the differences in water contents of interacting points lead the water flux into their connecting points. In order to calculate the micro-advection, a Peridynamics formulation for water flow is required to enable a coupled chemical transport and water flow is expanded. The following section describes the PD formulation for the water flow.

PD formulations for water flow
Similarly to Eq. (1), the Peridynamics water flux between the two interacting material points x and x ′ per unit volume per unit volume of x ′ is described as: Eq. (13) is integrated over the horizon of x to obtain the mass conservation for point x due to flow by all the bonds attached to the material point x within its horizon H Similar to Eq. (6), we assume that the relationship between the water content at point x and time t and the average water content in all the bonds connected to x is given by: From Eqs. (13)-(15), we obtain the mass conservation equation for liquid water flow for the material point x as: is the Peridynamics' micro-hydraulic conductivity function. The newly derived Peridynamics based formulations for the coupled water flow and chemical transport, Eqs. (11) and (16), do not involve any partial derivatives. It is also indicated that the proposed PD formulations can be applied everywhere in the domain including in cases of discontinuities of the parameters (i.e. d(⋅, ⋅, ⋅) and v(⋅, ⋅, ⋅)) and discontinuities of the domain.

Peridynamics' parameters
The Peridynamics' parameters, including micro-hydraulic conductivity κ(x ′ , x, t) and micro-diffusivity d(x ′ , x, t), are non-measurable parameters. However, the classical (local) material properties including diffusivity D and hydraulic conductivity K are experimentally measurable parameters. We present κ(x ′ , x, t) and d(x ′ , x, t) as functions of the hydraulic conductivity and diffusivity, respectively to develop the connection between the Peridynamics and the classical descriptions of the material constitutive relationships. Inside a nonlocal Peridynamics region, the micro-hydraulic κ(x ′ , x, t) conductivity and micro-diffusivity d(x ′ , x, t) depend on the distance between the interating points (Chen and Bobaru, 2015). In order to assess the level (significance) of contributions of all the bonds participating in calculating the volume-dependent properties, the influence function is introduced (Seleson and Parks, 2011;Silling et al., 2007). The most common of influence functions are the shape of "uniform" and "triangular" (see Fig. 2). The shape of the different influence function is important in the numerical approximations. However, investigating the specific effects that the shape of the influence functions on the overall behaviour of a system are yet to be fully investigated (Jabakhanji and Mohtar, 2015). Different influence functions have been tested by Jabakhanji and Mohtar (2015) which indicate that the solutions obtained by utilising "uniform" and "triangular" influence functions correlate closely with the solutions obtained by FEM. A uniform influence function is considered in this study for both micro-hydraulic conductivity and micro-diffusivity. By using a uniform influence function (d(x ′ , x, t) = d and as shown in Fig. 2a and b) the Peridynamics hydraulic conductivity and Peridynamics diffusivity will be independent to the distance between the points in the horizon which provides a simplified form of the solution.
The Peridynamics pure diffusive flux and the classic pure diffusive flux under steady-state conditions are considered to be equal (Bobaru and Duangpanya, 2010) which enables determining the Peridynamics micro-diffusivity. Based on Fick's law of diffusion, the steady-state concentration and flux for pure diffusion at position x can be given by: where, m and n are constants, J diff (x) is the classic chemical flux for pure diffusion. The peridynamic flux for pure diffusion along the 'c bond' connecting points x and x ′ can be obtained by reducing the advection term in Eq. (1): In order to derive the Peridynamics chemical flux in all the bonds attached to x in its horizon H x , Eq. (19) is integrated over the horizon of material point x. For 1D and 2D diffusion cases, the Peridynamics' chemical flux can be expressed as: By substituting the uniform influence function (see Fig. 2) and the steady-state concentration equation, Eq. (17), into Eqs. (20) and (21), we obtain: If the classic chemical flux, Eq. (18), and Peridynamics chemical flux, Eq. (22), at material point x are equal, the micro-diffusivity for 1D and 2D transport cases are obtained to be: Similar process is also applied to obtain the micro-hydraulic conductivity. Based on Darcy's law, the steady-state water flux due to a linear hydraulic potential field (H(x) = ax + b) at position x is: where, a is constant and q c (x) is the classic Darcy's flux. The Peridynamics' water flux (Eq. (12)) in all the bonds connected to material point x within its horizon H x , for 1D and 2D transport cases can be expressed as: Equivalence between the classic Darcy's flux, Eq. (24), and the Peridynamics flux, Eqs. (25) and (26), at material point x, provide the micro-conductivity relations:

Numerical solution
The numerical solution for the PD equations of water flow (Eq. (16)) and chemical transport (Eq. (11)) are developed for a finite number of PD particles. The domain is discretized using a uniform grid with spacing Δx. Each node has a ''volume" (length in 1D and area in 2D). Fig. 3 shows the 1D discretization around a node x. In this illustration, the points locate in yellow region are the points that their chemical concentrations and water contents are larger than that of the material point x. The points located in red region points are those points that their chemical concentrations and water contents are smaller than that of the material point x. N is the total number of particles in the whole domain. For a given PD particle x, the set of particles x ′ within H x is called P x and defined as: The spatial discretization of Peridynamics moisture flow, Eq. (16), and Peridynamics chemical transport in partially saturated soil, Eq. (11), for the given PD particle can be approximated as: The implicitly coupled method is used for solving the coupled problems. In this method, equations for water flow and chemical transport are solved simultaneously. The spatially discretized equations for water flow and chemical transport can be combined in a matrix form. The matrix equation can conveniently be expressed as: where A, B and C are the matrices of coefficients and ϕ is the vector of variables.
We use the forward Euler method for the time integration of Eq. (31): Rearranging Eq. (32) gives: For the numerical treatment of the implicitly coupled water flow and chemical transport Peridynamics system of equations, the same time stepping is adopted. The equation of water flow is solved for the fluid field and the chemical transport equation is solved for the concentration field. In order to approximate the solutions to both equations, we adopt an explicit time stepping scheme.
In order to illustrate the numerical implementation, the solution applied to the 1D Peridynamics coupled equations are discussed. Similar procedure has been applied to the 2D problem. In this study, the horizon size of δ = 3Δ has been selected which is based on the efficient numerical convergence study confirmed by various researchers (Agwai, 2011). Δ is the spacing between the material points. As shown in Fig. 3, the material point of interest x is denoted by x i which interacts with the three points to its left and right. Hence, the set point x′ from left to right within the horizon of x i are x i-3 , x i-2 , x i-1 , x i+1 , x i+2 , x i+3 , as shown in Fig. 3. The Peridynamics' water flow and chemical transport formulations, Eqs. (29) and (30), at the nth time step (i.e., the current time step), can be presented in a matrix form as: The details of matrix [K] and [M] are provided in Appendix C. The approach proposed by Bobaru and Duangpanya (2012) is adopted for prescribing the boundary conditions. The nonlocal Dirichlet boundary condition is applied directly at the boundary nodes. The nonlocal Dirichlet boundary condition converges to the local one as the size of the horizon becomes smaller and approaches zero. The Neumann boundary condition is applied by equating the known volumetric fluxes with their peridynamic equivalents.
The peridynamics bonds caused by either concentration or hydraulic gradient that cross a pre-existing crack are assigned different micro properties, specifically zero micro diffusivity and zero micro hydraulic conductivity. Consider an insulating crack defined by the segment AB and peridynamics bonds between material points x and x′ denoted by f x x′. If the line segment AB intersects the line segment xx′, the peridynamics bonds f x x′ is set to zero; mathematically if (xA × xx For modelling heterogeneous inclusions in the system (e.g. oil and rock segments), a sub-domain is generated to represent the heterogeneity. The particles located inside the sub-domain are assigned with different properties (e.g. micro diffusivity and hydraulic conductivity). Therefore, different peridynamics bonds can be defined to describe the interaction between points of different phases. For example, the peridynamics bonds between the material points x and x′ are defined as • f rp , the material points x and x′ are located in rock and porous media, respectively. • f rr , the material points x and x′ are both located in rock domain.
• f pr the material points x and x′ are located in porous media and rock, respectively • and f pp , the material points x and x′ are both located in porous media domain.

Verifications
In this section, there are two verification examples including i) chemical transport in a saturated porous medium and ii) coupled water flow and chemical transport in unsaturated porous system. Chemical transport in a saturated porous medium will be first presented. The proposed model will be verified by an existing analytical solution for 1D problem. The problem is then extended to a coupled water flow and chemical transport in unsaturated porous media. The proposed PD solutions will be compared with alternative finite element analysis of the same problem.

Chemical transport in a saturated porous medium
The verification tests presented in this section focused on testing the PD model for pure chemical transport in a fully saturated system. Two problems for which analytical solution exists are considered: Case 1, where neither diffusion nor advection is dominant; and Case 2, where the transport is advection dominated.
Diffusion and advection mechanisms for chemical flow are considered. A constant velocity of liquid is considered in the simulations. The length of domain is 10 cm. The initial concentration of chemical in the domain is considered to be zero. The boundary condition at x = 10 cm is assumed to be zero flux boundary conditions (i.e. Neumann boundary condition). The boundary condition at the injection side (x = 0) is maintained at C = C 0 . In peridynamics approach, Δ is set to be 0.01 cm and the horizon sizes of 0.03 cm is considered. The classical solution for this problem can be found in (Cleary and Ungs, 1978;Ogata and Banks, 1961) as: For case 1, the diffusion coefficient and advection are assumed to be D = 1.0 cm 2 /s and V = 1.0 cm/s, respectively. Fig. 4 shows the chemical concentration profiles at 1 s, 2 s and 3 s, obtained by the proposed formulation (solid lines) and the results by the analytical solution (symbols). The simulation results with the PD formulation are clearly in excellent agreement with the analytical solution.
For case 2, the diffusion coefficient and advection are assumed to be D = 0.01 cm 2 /s and V = 1.0 cm/s, respectively. The Peclet number, Pe = VL D = 1000≫1. This smaller dispersivity or advection-dominated transport case is in principle difficult to solve as numerical oscillations develop when advection dominates. Nevertheless, the PD formulation is able to capture the transport successfully. The chemical concentration distributions at 1 s, 2 s and 3 s obtained by the PD model and those obtained by the analytical solution are shown in Fig. 5. It can be seen that the results are again in very good agreement with the analytical solution.
The verification examples demonstrate that a robust implementation of the theoretical formulation of chemical transport in the computational model has been achieved.

1D water flow and chemical transport in unsaturated porous media
The problem considered is the water flow and chemical transport in a one-dimensional domain of partially saturated medium with a length of 10 cm (L = 10 cm). The parameters used in the simulation of water flow and chemical transport are listed in Table 1.
The initial water content is assumed to θ r and the initial chemical concentration in the domain is considered to be zero. The boundary conditioned applied at x = 10 cm for chemical transport and water flow analysis are C = 0 and θ l = θ s , respectively. The boundary condition at injection side (x = 0) is maintained at C = C 0 and θ l = θ s . Δ is 0.05 cm and a horizon size of 0.15 cm was used for the PD simulations. The PD simulation results (Fig. 6) are compared with analysis of the same problem solved by Finite Element method which is implemented in Comsol Multiphysics with the element size of 0.05 cm. It can be seen from Fig. 6 that there is a good agreement between PD solutions and FE solution indicating that the proposed PD formulations can be used to simulate the water flow in unsaturated porous medium.

2D water flow and chemical transport in unsaturated porous system
For the 2D simulation case, the domain is 10 cm long and 10 cm wide. The parameters of simulation used are listed in Table 1. The initial chemical concentration and water content in the domain are considered to be zero and θ r , respectively. The boundary condition at x = 10 cm is assumed to be zero flux (impermeable boundary). The boundary condition at the injection side (x = 0) is divided into two areas. Area One is 4 cm wide at x = 0 cm (left boundary). The area is between coordinates (0, 3) and (0, 7). Zone Two covers the excluding area one at the left boundary of domain (Fig. 7). At the injection boundary in area one the chemical concentration and moisture content are maintained at C = C 0 and θ l = θ s , respectively.
The chemical concentration and moisture content in zone two are assumed to be C = 0 and θ l = θ s , respectively (Fig. 7). Δ is 0.1 cm and a horizon size of 0.3 cm is assumed for the PD simulations. Similar to previous verification case, the PD results are compared with Finite Element analysis of a same problem. The same time step of Δt = 10 − 4 d is used for PD and FE.
The results of the saturation profile and chemical concentration profile along the centre line (x = 5 cm) from the PD simulation are compared with those from FE solution at t = 5, 10 and 15 d in Fig. 8. Fig. 9 shows the comparison of peridynamics solution against the numerical solution obtained by FE analysis along x = 0.02 m. It can be seen that the results obtained by the proposed PD are in excellent agreement with the FE solution for 2D cases (Figs. 8 and 9). A difference of saturation and concentration profile can be found in Figs. 6 and 8. The   Fig. 4. Comparison between analytical solution and PD analysis (this study) for simulation case 1.     results in Fig. 6 show that chemical migration is faster than water transport; however, the results in Fig. 8 show a contradicting trend. These phenomena may be explained by dilution effects and transversal dispersion. Moreover, the results in Fig. 8 indicated a faster migration of chemical and water compared to the results observed in Fig. 6. This is because the source of chemical and water for 1D problems is a point injection; however, the injection position for 2D problems ranges from (0, 3) and (0, 7). The computational times required for the PD and FE analysis for this simulation case study using a I7-6200 CPUs is 16 min 38 sec and 8 min 9 sec, respectively. The computational time for the PD analysis compared to the FE is approximately twice (for this case) and this is due to the fact that the PD calculations involves analysis of interactions between a large number of collocation points (Fig. 3). However, it is noted that the PD can provide solutions for discontinuous or strongly heterogenous problems without any extra requirements. This point will be further clarified in Section 6. Moreover, the PD solutions can readily be implemented in a parallel programming environment by which a significant reduction in run times can be expected.

Model application
In this section, we present three examples to show the capacity of PD models to capture the transport features in discontinuous and heterogeneous unsaturated porous media. A simulation case for coupled water flow and chemical transport in unsaturated porous media with multiple insulating cracks will be presented first. The case is then extended to consider propagating cracks. Finally, two case studied of analysis of the coupled water flow and chemical transport model in heterogeneous unsaturated porous media will be shown.
The domain is 10 cm long and 10 cm wide. The initial chemical concentration and moisture content in the domain are considered to be zero and θ r , respectively. The boundary at x = 5 cm is assumed to be impermeable, i.e. natural (Neumann) boundary condition of zero fluxes. The boundary at the injection side (x = − 5 cm) is maintained at C = C 0 and θ l = θ s , i.e. essential (Dirichlet) boundary condition with constant concentration and moisture for chemical transport and water flow, respectively. The properties of the porous medium and the parameters used for the simulations are listed in Table 1. The problem is analysed with horizon sizes δ = 0.15 cm and Δ = 0.05 cm.

Chemical transport in unsaturated porous media with insulating cracks
Solving the mathematical equation of transport processes with cracks and fractures using the classical local formulations is computationally challenging due to the difficulties in solving derivatives over discontinuities. The PD eliminates such difficulties by its non-local and integrating approach. In this section, water flow and chemical transport in a system with pre-existing cracks (discontinuity) are studied. The case study was designed to demonstrate the capability of proposed PD model in solving the transport processes in a discontinuous system with insulated cracks. A 2D domain with two single insulated vertical cracks of length 2 cm is considered, as illustrated in Fig. 10. The implementation of peridynamics bonds that cross the crack is described in Section 4. The concentration distribution and saturation profile are shown in Fig. 10.
It can be seen that PD models capture successfully the shielding effects of the cracks: both the water flow and the chemical concentration profile are strongly affected by the presence of cracks. Sharp jumps of concentration and saturation near the cracks can be observed. For example, the saturation at the points (0.1, 0.2− ) and (0.1, 0.2+) is 0.94 and 0.62, respectively. The jumps are natural consequence of the zero peridynamics bonds caused by concentration gradient or hydraulic across the cracks. This saturation behaviour in the presence of cracks/ fractures is challenging for the local methods including FEM and FDM, where the cracks must be explicitly modelled as narrow slits and meshed. This is because in the local formulations the constitutive parameters, such as diffusivity and conductivity, are volume-based and hence cannot be assigned to interfaces. For example, the concentration distribution and saturation profile shown in Fig. 11 is obtained by FEM. It can be found that geometry of the domain is divided into two sections. One is the porous medium and another one is the insulating cracks which are separated and marked as an extra boundary. The insulating cracks are considered as a 0.3 cm × 2.5 cm slits in FEM. A slight difference can be found between Figs. 10 and 11. This may be due to the different approaches adopted to represent the insulating cracks by PD and FEM. The crack is represented naturally with the broken bonds in PD; however, extra boundaries are introduced to form the cracks in FEM.

Chemical transport in unsaturated porous media with propagating cracks
In this section, the problem described in Section 5.1 is extended to study a dynamic discontinuous system for which the length of discontinuities in the domain is enlarged with time. We did not aim to describe the mechanics of crack propagation in this work, but to explore the The concentration distribution and saturation profiles at 5 d, 10 d, 15 d and 20 d are shown in Fig. 12. It can be seen that the effect of cracks increases from being practically insignificant at 5 d to very significant at 15 d, where strong discontinuous profiles of saturation and concentration near the crack at x = 0.1 cm are observed. The distribution profiles of saturation and concentration are orthogonal to the crack surfaces. The shielding effect induced by the growing cracks is correctly captured by the peridynamic solution. A nonnegligible difference can be noticed between Figs. 10 and 12 at 20 d. A faster transport phenomenon is observed in Fig. 10. A sharper jump of concentration and saturation near the second crack is observed in Fig. 12 compared to the results in Fig. 10. For example, the saturation difference along the cracks can reach 0.32 in Fig. 10; however, this value increases to 0.68 in Fig. 12. The is due to the fact that a larger crack with a length of 6 cm is formed in Fig. 12 compared to the length of 2.5 cm in Fig. 10. The larger cracks show stronger shielding effects. Notably, water flow and chemical transport through a system of propagating cracks is even more challenging for local methods, where adaptive mesh refinement may be required to track the crack extension.
The examples provided in Sections 5.1 and 5.2 demonstrate the ability of the proposed PD formulation to handle complex static and dynamic discontinuous systems. The PD models deal with the problem naturally. PD-based extension to include the mechanical and geochemical reaction processes for a fully coupled hydro-chemical-mechanical model is currently under development.

Coupled water flow and chemical transport in heterogeneous system
The effects of heterogeneous system on water flow and distribution of chemicals are studied in this section. The case studies were designed to demonstrate the capability of proposed PD model in solving the transport processes in heterogeneous media which has permeable or impermeable inclusions. For this purpose, a 2D domain with randomly distributed circular inclusions is considered, as illustrated in Fig. 13.
In order to assign impermeable inclusions, the diffusion coefficient and hydraulic conductivity are assumed to be zero for the inclusions. Fig. 14 shows distributions of degree of saturations and concentration of chemicals in the domain at t = 40 d. The grey circles indicate the location of inclusions. The chemical transport and water flow are restricted only to the gaps between the inclusions. The PD model has successfully analysed the discontinuity in the domain. The chemical concentration and saturation of the impermeable circles maintain C = 0 and S r = 0, which are equal to the initial values. This shows that the implementations of PD bonds for heterogeneous system are able to solve the transport phenomena with strong impermeable inclusions. It is shown that migration of chemical and water is impeded by the existence of impermeable inclusions. For example, there is no chemical transport observed at the back of the circle and two ellipses on the left-hand side. The furthest transport distance of water and chemical is observed at y = 0.05 m. It can be found that y = 0.05 m is the furthest away from the impermeable inclusions so the shielding effects is the weakest at this line and a faster transport is observed. The results indicated that the   A similar domain is analysed but permeable inclusions are prescribed which have smaller diffusion and hydraulic conductivity coefficients than the rest of the domain (matrix). To prescribe the diffusion coefficient and hydraulic conductivity we use a random number (from a uniform distribution) ranging between (0, 1) for each inclusion. The black dashed circles indicate the location of permeable inclusions.    15 shows the concentration distribution and degree of saturation in the domain at t = 40 d. The results show that the randomly assigned diffusion coefficients and hydraulic conductivity lead to different levels of penetrations in different inclusions. For example, the back diffusion is found in the circle on the left-hand side. The diffusivity of chemical in porous media is larger than that in the inclusion, which leads to a faster transport of chemical in porous media. However, the concentration difference between porous media and permeable inclusions results in the back-diffusion phenomena observed in Fig. 15. It is noteworthy that the water saturation is almost 1 in one circle and two ellipses on the left-hand side of the model, while concentration is almost zero in these structures. This observation can be explained by the following two points. One is that water transport is faster compared to chemical migration. It is governed by the input diffusion coefficients and hydraulic conductivity. The same trends are observed in Figs. 8, 10-12 and 14. Therefore, water will be saturated in porous media in a shorter time than that of chemical and permeate into the inclusions. Another point is that the randomly assigned hydraulic conductivity in the inclusions is relatively large. It leads to a faster penetration in the inclusions. The penetrations distance of water and chemical is larger in Fig. 15 than that in Fig. 14. For example, the penetration distance of chemical is around 6.5 cm and 7 cm at y = 5 cm in Figs. 14 and 15m, respectively. This shows that the impacts of inclusions are governed by the hydraulic conductivity and diffusivity of the inclusions. A smaller hydraulic conductivity and diffusivity of the inclusions will show a more significant blocking effect. The results from Figs. 14 and 15 demonstrate how the inclusion in overall control of transport processes.

Discussion
Heterogeneities such as permeable/impermeable inclusions and cracks in a porous medium result in discontinuities in the concentration and flow field (see . The classical local model faces difficulties simulating such problems since spatial derivatives fail at the discontinuities. In classical approach such as FEM and FVM, special interface/boundary conditions must be introduced to satisfy conservation principles (Fig. 11). However, in non-local peridynamics model, transport equations are valid everywhere, regardless of the presence of material discontinuities/heterogeneities. This makes nonlocal models naturally applicable to sharp interface problems without any additional interface conditions.
The heterogeneities/discontinuities of porous media are defined by using different bonds across them without any additional modification to the transport formulations or any extra computational expense. The results from Fig. 10 show a jump of saturation and concentration around the discontinuous cracks region, which demonstrate that the discontinuous features in the domain are captured by the peridynamic model. For a dynamic crack growth porous media (Fig. 12), the PD shows a strong ability to capture the features induced by the propagating cracks without any refinements and modifications. It further demonstrates the advantages of PD for solving heterogeneities/discontinuities in porous media.
Important insights toward transport in heterogeneities porous media are obtained in Section 5.3. The simulations in Fig. 15 shows the backflow and back-diffusion phenomena which are caused by the existence of permeable inclusions, which promote the appearance of fast paths as well as regions of hold up. The extra modifications, such as stochastic approaches and continuum-scale dual porosity, are required to capture these phenomena by using the local method; in contrast, the non-local PD model are naturally applicable to such problems and capable to capture these abnormal phenomena. Development of multi-chemical reactive transport in discontinuous porous media, which is induced by stress or geochemical reaction, is of interest for future work.

Conclusions
Bond-based non-local formulations of water flow and chemical transport in porous media are proposed. Water flow is partially coupled with chemical transport by an implicit numerical scheme to allow for solving coupled problems of practical importance.
The formulations are implemented within a computational framework (Pyramid) and verified against a set of benchmarks. The verifications demonstrate excellent agreements between the proposed approach and the classical/local formulations for homogeneous media.
Extended capability of the model is demonstrated by analysis of coupled chemical transport and water flow problem in 2D domains with stationary and propagating cracks and various randomly distributed inclusions. The simulations illustrate how the model can capture naturally the effects of static and dynamics discontinuities and heterogeneities, which is problematic or computationally demanding for numerical schemes based on local formulations.
Overall, it is demonstrated that the proposed nonlocal formulationimplementation framework is efficient and effective in analysis of coupled water flow and chemical transport problems. The model is applicable to multi-physics problems in porous systems where coupled hydro-chemical degradations can create discontinuity or heterogeneity. The addition of mechanical deformation, fracture nucleation and erosion is a subject of ongoing work.

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.