Peridynamic modelling of clay erosion

Experimental observations of the erosion of clays indicate a high degree of coupling between the hydro-chemical processes in the clay-water system and the mechanics of erosion. In contrast with the...


INTRODUCTION
Robust prediction of clay erosion is important for a range of geotechnical engineering problems: from clay erosion in embankments (Fujisawa et al., 2009), through erosion of clay in the geological disposal of nuclear waste (Schatz et al., 2013), to erosion of geosynthetic clay liners (Ashe et al., 2014). Of particular importance and recent interest is the erosion of clay buffer and backfill in the context of geological disposal of high-level nuclear waste (HLW). The clay barriers are key components of the engineered barrier system (EBS) in geological disposal concepts and provide critical functions in the safe isolation of HLW. During the long-term operation of disposal facilities, the clay buffer/backfill can be subject to erosion induced by the hydro-chemical interactions at the interface between the clay and the fractured host rock. Localised interactions of the clay with percolating water at fractures/fissures/pre-cracks in the host rock, varying in size from micro-fractures (below 1 mm) to relatively large fractures (several millimetres and more), can result in gradual loss of the barrier and a co-transport of radionuclides into the biosphere. The aim of this paper is to address the need for a realistic and robust modelling of clay erosion. This is accomplished by presenting a set of new formulations, based on the non-local continuum mechanics of peridynamic (PD), which is utilised, for the first time, to couple the hydro-chemical processes with the erosion of swelling clays.
Experimental studies of the erosion of smectite-rich clays (highly swelling clays) have shown that the overall behaviour is strongly dependent on the fluid geochemical composition and exchangeable ions, and on the fluid velocity (Schatz et al., 2013;Smith et al., 2017). The erosion of swelling clays involves three physical processes: (a) free swelling, representing a continuous transformation of solid-type phase to gel-type material; (b) clay particle detachment; and (c) migration of detached particles with the water flow in the fracture. The semi-empirical or phenomenological models proposed to date do not fully address all phase changes accompanying the underlying processes. The existing models address mostly the individual processes, for example: generation of clay gels (Pusch, 1999); mechanical erosion caused by water flow (Grindrod et al., 1999); erosion of particles and the dispersion-flocculation behaviour of clay particles (Kurosawa et al., 1999); and expansion/swelling of the clay (Navarro et al., 2017;Asensio et al., 2018). Pioneering works by Neretnieks et al. (2009) and Liu et al. (2009aLiu et al. ( , 2009b introduced the dynamic force balance model in simulating clay expansion and extrusion. A version extending the model by flow and colloid diffusion was presented by Moreno et al. (2011), while a recent development (Neretnieks et al., 2017) addressed the strong hydro-chemical effects on erosion of smectite-rich clays observed in experiments (e.g. Schatz et al., 2013).
A critical missing element in the existing models is related to the detachment of particles from the gel phase at the solid/fluid interface and its strong coupling with hydrochemical effects. The current authors' modelling approach will introduce a new criterion for detachment based on the particle interaction forcesthat is, an intrinsic criterion, rather than one prescribed phenomenologically. The particles' detachment is inherently a discontinuous process taking place at the moving boundary of the continuously swelling phase (gel). The non-local formulations proposed in this work, based on PD theory, are suitable to tackle the complexity of this problem. The PD theoretical framework replaces the partial differential equations of the classical local theory by a set of integral-differential equations (Silling, 2000). This results in a mathematically consistent formulation, even where strong discontinuities appear due to breakage and fragmentation of the material. The PD theory has successfully been applied in different areas, such as hydraulic fracturing (Ouchi et al., 2015), thermomechanics of fracture development (Oterkus et al., 2014), membranes and fibres (Silling & Bobaru, 2005), as well as the dynamic fracture of nanofibre (Bobaru, 2007).
The aim of the work is achieved by developing and integrating bond-based PD formulations for free swelling of clay as a diffusion process and for transport of detached particles, with the important incorporation of the new detachment model. The PD formulations have been shown to provide the computational capability required to address diffusion-type problems for evolving discontinuities or phase changes (e.g. heat transfer in evolving cracking system, pitting corrosion damage and mechano-chemical damage) (Bobaru & Duangpanya, 2012;Chen & Bobaru, 2015;Jafarzadeh et al., 2019aJafarzadeh et al., , 2019b. Furthermore, the PD modelling approach allows for autonomous evolution of the moving interface without any extra modelling efforts (e.g. interface conditions), which is important for modelling the moving boundary of the continuously swelling phase. The PD formulation can offer a computationally robust solution for describing the free swelling where the parameters of the diffusion process vary by several orders of magnitude during the soil expansion. Approximation of such a level of non-linearity by classical numerical methods, such as the finite-element method (FEM) and the finite-volume method (FVM), may yield gross computational errors (Jafarzadeh et al., 2019c). The PD erosion model is validated against the results of a set of experimental benchmarks, first to establish the realism of the detachment model, and then to predict the behaviour of compacted smectite-rich clay under eroding environments.

LOCAL FORMULATIONS OF EROSION PROCESSES
From a macroscopic (continuum) perspective, the erosion of clay by water can be described as a series of phase changes: an initially solid phase is hydrated, leading to a continuous transformation into a swelling paste, which in turn is continuously transformed into a gel, and eventually discontinuously transformed into dispersed colloidal particles (sol) in the water. The material is initially partially saturated and the free swelling model approach describes the evolution of the solid volume fraction with time from the initial saturation stage until a fully hydrated system, with an assumption that the process follows a diffusion-type behaviour. In this model, the evolution of swelling in a relatively thick fracture (.1 mm) is considered, where the friction between the material and boundary does not induce significant effects on the overall process (Neretnieks et al., 2009). Figure 1 provides a schematic diagram of the erosion of swelling clay and its corresponding processes. The free swelling step is characterised by continuous reduction of the clay density. The distinction between the gel and sol phases can be made on the base of previous studies on smectite-rich clays.
Forces acting between the particles are shown in Fig. 1; these include the diffusion forces (F T ), attractive van der Waals forces (F A ) and repulsive electrical double-layer forces (F R ). These forces can be described by Liu et al. (2009a): where ϕ s is the clay (solid) volume fraction; k B is the Boltzmann constant; T is absolute temperature; A H is the Hamaker constant; S p is the surface area of the particles; h is the separation distance between the flat particles; δ p is the particle thickness; c is the ionic concentration in the pore Clay core Fig. 1. Conceptual description of clay erosion by swelling, particle detachment and transport of detached particles SEDIGHI, YAN AND JIVKOV water system; and R is the universal gas constant. The function y m represents the scaled potential in the midpoint between parallel plates of, for example, a fracture, as described by Liu et al. (2009a): where y m 1 ¼ 4 tanh À1 2 tanh The dimensionless surface potential of the isolated plate (y ∞ 0 ) and the Debye length (κ) in equations (4)-(6) are given by where F is the Faraday's constant; σ 0 is the surface charge density; z is the valence of ions in the pore system; and ε 0 ε R is the dielectric constant. It is noted that the calculation of F R is based on the 'compression' approach, which is applicable to a wide range of parameter conditions in comparison with other available approximate expressions (Laxton & Berg, 2006;Liu et al., 2009a).

Free swelling of clay
The free swelling is treated as a solid diffusion process following the approach proposed by Liu et al. (2009a). It is noted that the free swelling model proposed by Liu et al. (2009a) is developed for the colloidal system of clay. However, applying a diffusion-type description of free swelling for compacted swelling clay is not far from experimental observations (Moreno et al., 2011;Neretnieks et al., 2017). In the absence of gravitational forces and buoyant effects, the mass balance equation for clay solid diffusion can be simplified to where t is time, χ is the sum of the particles' energies, and f r is the friction coefficient between particles and water. The friction coefficient is given by Moreno et al. (2011) where r eq is the equivalent radius of the non-spherical particles; V p is the volume of the particles; k 0 is the pore shape factor; τ is the tortuosity of the flow channel in the clay gel; a p is the specific surface area per unit volume of particles; and η w is the dynamic viscosity of water. The particles' energy (χ) due to the interaction forces F T , F A and F R , is given by Liu et al. (2009a) The separation distance between the flat particles (h) is described by Laxton & Berg (2006) h where ϕ s max is the maximum clay solid volume fraction.

Particle detachment
The process of the detachment of particles from the solid/fluid interface is central to clay erosion. The detachment is caused by shear stresses on the particle-particle bonds of the swelled clay (gel) induced by the fluid flow along the interface. This shear stress can be related to the flow velocity by Stoke's law (Kurosawa & Ueta, 2001) where τ f is the shear stress; D p is the particle diameter; and v l is the fluid velocity. A detachment model is required to link dynamically the solid swelling and particles' release processes at the solid/liquid interface. A general approach is to compare the shear stress at the solid/liquid interface induced by water flow with a material parameter limiting the solid behaviour of the gel, typically called the yield stress (Laxton & Berg, 2006;Sane et al., 2013). At shear stresses lower than the yield stress, the gel behaves like a solid. At shear stresses above the yield stress, the clay particles exhibit sol behaviour and can be transported away by seeping water. This concept is in agreement with experimental observations (Eriksson & Schatz, 2015). However, its application in a continuum modelling framework requires knowledge of both the yield stress and the solid/liquid interface, which remains a challenge. Based on experimental observations, Neretnieks et al. (2017) developed a rim-region model, which imposes an external boundary as a detachment interface with a value of ϕ R (solid volume fraction at rim). Although their model was able to capture a number of experimental observations, it did not provide physical explanations of the mechanisms involved.
The detachment model proposed in this work is significantly different in two ways: (a) the yield stress is not prescribed but emerges from the balance between the shear stress on a particle and a cohesive stress arising from its interactions with the particles (in its horizon in the PD formulation) and (b) the solid/liquid interface is evolving naturally with the erosion process, rather than being assumed. With this formulation, the yield stress is a variable dependent on the local conditions at the solid/liquid interface, which equals the balancing shear and cohesive stresses.
The diffusion forces F T , given by equation (1), result from changes in the chemical potential and cause a corresponding stress expressed by Goodwin & Hughes (1992) Using the forces F A and F R given by equations (2) and (3), the cohesive stress is approximated by The detachment of clay (gel) particles is then postulated to take place if the shear stress induced by the water flow on a particular boundary particle exceeds the cohesive stress.

Migration of sol particles
The detached particles (clay sol) are transported and dispersed in the fluid by mechanisms of advection and dispersion where D d and D m are the diffusion and the mechanical dispersion coefficients of detached particles in the liquid, respectively. The fluid velocity can be expressed by Darcy's law where p is the fluid pore pressure; T w is the fracture transmissivity for water; and η s is the viscosity of clay sol/gel. The latter is given by Moreno et al. (2011) where the co-volume fraction (ϕ s cov ) is defined by Neretnieks The diffusivity of detached particles is calculated based on the Stokes-Einstein equation, which is valid for large spherical particles (or molecules). The dispersion coefficient is considered to be a function of the flow velocity. The sum of the diffusion and dispersion coefficients is given by Bird (2002) where α L is the dispersivity parameter. The clay erosion is a discontinuous phase change process with highly non-linear coefficients. A non-local formulation development will offer significant benefit by introducing a detachment functions, which enables the discontinuous detachment to evolve autonomously without the need of tracking surfaces.

BOND-BASED PD FORMULATION OF EROSION PROCESSES
In PD, a body occupying a region (R) is considered to be a collection of infinite particles, as illustrated in Fig. 2(a). In the bond-based PD, the term 'bond' refers to the interactions between the material points at spatial positions x and x′. The bond between two interacting material points is called a 'mass bond' ( f h ), see Fig. 2. The peridynamic mass flux per unit volume along an 'f-bond' depends on the distance between the points x and x′. The material point x interacts with (and is connected to) all materials points x′ within a certain finite region H x , called the horizon of material point x. The radius of the horizon is denoted by δ, and ξ represents the distance vector between the two material points x and x′that is, ξ ¼ (x′ À x).
Corresponding to the three distinctive processes involved in erosion (i.e. free swelling, particle detachment and detached particle transport), three different PD bonds are defined: solid-solid bonds, interfacial bonds and liquidliquid bonds (illustrated in Fig. 2(b)). The solid-solid bonds control the free swelling process and liquid-liquid bonds govern the transport of detached particles. The interfacial bonds represent the solid-liquid interaction at the interface and control the detachment process (Jafarzadeh et al., 2019a(Jafarzadeh et al., , 2019b. Once the shear stress exceeds the cohesive stress, the solid particles detach and become liquid particles. The erosion process will switch from free swelling to detached particle transport as the phase change proceeds. Accordingly, the interface evolution (detachment) is computationally processed in the PD model of erosion. It is noted that the free swelling and detached particle transport follow two different governing equations. Therefore, the detachment of particles at the moving interface is considered as a discontinuous process.

PD formulation for clay expansion and detachment
The mass conservation equation for the bond between two material points x and x′, describes the temporal variation of solid content by Zhao et al. (2018) whereφ s ðx; x′; tÞ is the average clay volume fraction along the bond between points x and x′; J s (x, x′, t) is the solid mass diffusion flux; and ξ= ξ k k is the unit vector along the bond. The gradient of clay solid content (clay concentration gradient) between two particles induces clay expansion along the bond connecting these particles. The rate of clay expansion along the bondthat is, the solid flux J s (x, x′, t)is, according to Bobaru & Duangpanya (2010) where ϕ s (x, t) and ϕ s (x′, t) are the clay volume fractions at points x and x′, respectively, andD s ðx; x′; tÞ is the average solid diffusivity between points x and x′that is, D s ðx; x′; tÞ ¼ 0Á5 D s ðx; tÞ þ D s ðx′; tÞ ½ . The conservation of mass at point x involves the fluxes in all the bonds adjacent to x (bonds to particles within the horizon, H x ) and is obtained by integrating equation (21) over the horizon, H where V x′ is the horizon volume (area in two dimensions, length in one dimension) of particle x′.
The relationship between the concentration at point x and time t and the average concentration in all the bonds adjacent to x is given by Bobaru & Duangpanya (2010) ð where V H x is the horizon volume of particle x.
Combining equations (22)-(24) leads to the following equation for the evolution of clay solid content at particle x The detachment of clay particles at the solid-liquid interface is related in this work to force bond damage, as previously modelled in PD (Silling & Bobaru, 2005;Bobaru & Duangpanya, 2012). A similar concept of using concentration-dependent damage to simulate a corrosion process by way of a peridynamic non-linear diffusion equation is presented by Jafarzadeh et al. (2019aJafarzadeh et al. ( , 2019b. A phase change process autonomously evolves when concentration in the solid drops below the concentration in the electrolyte, and the phase change evolving is governed by SEDIGHI, YAN AND JIVKOV solid-liquid bonds. Similarly, the solid-liquid bonds used in this study control the detachment processes. However, different definitions for solid-solid and liquid-liquid bonds are used in this study. Solid-solid bonds and liquid-liquid bonds govern the swelling process and the transport of detached particles, respectively. The detachment of clay particles is postulated to occur when the shear stress (τ f ) induced by the flow at the gel-fluid boundary exceeds the cohesive stress (τ c ) given by equation (15). To establish the relation, a detachment function is introduced here The detachment function implies that when the cohesive stress calculated from the mass bonds of solid-solid points is smaller than the shear stress induced by the seeping water, the solid material points change to liquid points and follow the governing equations of transport processes. Incorporating the detachment function μ(x, x′, t) in equation (25) leads to the peridynamic equation for coupled clay expansion and detachment of particles The PD microscopic diffusivity is defined as which transforms equation (27) into PD microscopic diffusivity is determined from measurable material diffusivity. Using the approach proposed by Bobaru & Duangpanya (2010), the relationships between PD microscopic diffusivity and macroscopic diffusivity are derived for one-dimensional (1D) and two-dimensional (2D) deformation as Details of the derivations of equations (30) and (31) are provided in the online supplementary material (S1).

PD formulation for transport of detached particles in water
The mass conservation equation, equation (16), for the bond along the two material points x and x′ can be presented within the PD formulation as (Zhao et al., 2018) @φ s x; x′; t ð Þ @t  where J l d (x, x′, t) is the mass flux caused by dispersion and J l a (x, x′, t) is the mass flux caused by advection. In a similar way to equation (22), the PD dispersion flux J l d (x, x′, t) is written as follows: where D l (x, x′, t) is the diffusivity of the 'm-bond' between the material points of x and x′. The PD advection flux J l a (x, x′, t) is given by the equations (Zhao et al., 2018) where V l (x, x′, t) is the advection of the 'm-bond' between the material points of x and x′.
Integrating over the horizon region of point x, H x , on both sides of equation (32) gives Substituting equations (24), (33) and (34) into equation (35) gives the PD description of transport of detached clay particles The PD microscopic diffusivity and microscopic advection of detached particles are defined as which transforms equation (35) into: The PD microscopic diffusivity for 1D and 2D cases is given in equations (30) and (31). By using the same approach for the proposed of microscopic diffusivity, the relationships between PD microscopic advection and macroscopic advection for 1D and 2D cases are derived as v l x; x′; t ð Þ¼ 4V s πδ 2 Details of the derivations for equations (40) and (41) are provided in the online supplementary material (S2).
Notably, the newly derived PD formulations for clay solid swelling and detached clay particles transport, equations (29) and (39), do not involve any spatial derivatives. This indicates that the proposed PD formulations can be applied everywhere in the domain, including in cases of discontinuities and heterogeneities of the parameters (i.e. d s , d l an v l ) and discontinuities and heterogeneities of the domain.

Numerical solution of the PD formulations
Numerical solutions of the equations describing clay expansion, equation (25), and transport of detached clay particles with the fluid, equation (35), are developed for an infinite number of particles. The domain is discretised using a uniform grid with spacing Δ, and particles are placed within the grid cells. A finite number of nodes is obtained as the result of the spatial discretisation. Each particle/cell has an associated volume V ip (length in one dimension and area in two dimensions). In this study, the horizon size of δ = 3Δ has been selected. Fig. 3 illustrates a discretisation in the vicinity of a node x i in a 2D domain. ξ ip represents the distance vector between the two material points x i and x p .
The spatial discretisations of equations (25) and (35), using the mid-point rule, are The forward Euler method is used for the time integration Online supplementary material S3 provides the details of the numerical solution applied to the 1D PD. The online supplementary material S4 also provides the flowchart of the erosion calculation. Fig. 3. Illustration of discretisation around node x in a 2D domain SEDIGHI, YAN AND JIVKOV

RESULTS AND DISCUSSION
Validation of the cohesive stress proposal First, the current authors present a validation of the proposed cohesive stress formulation by comparison with experimental studies reported by Eriksson & Schatz (2015). The cohesive stress induced by the different particle interaction forces is calculated by equations (2), (3), (14) and (15). The overall cohesive stress emerges as a combination of the forces. The experimental study reported by Eriksson & Schatz (2015) involved samples of clay with solid contents ranging from 0·1 to 15% in aqueous solutions containing 1 and 17 mM sodium chloride (NaCl). The materials' properties and parameters used in this section are taken from several literature sources and are shown in Table 1. Figure 4 shows the calculated cohesive stress as a function of particle distance for the case of the system containing 1 mM and 17 mM sodium chloride (as an example to demonstrate the component of cohesive stress model). The repulsive forces (F R ) are larger than attractive and diffusive forces (F A , F T ) at small distances, and become very weak as the distance between the particles increases. The overall cohesive stress (τ c ) emerges as a combination of the forces.    It is noted that τ c increases with the increase of the ion concentration as the repulsive forces (F R ) are governed by the ion concentration. It has been reported that the thickness of a single particle (δ p ) for sodium smectite (WyNa) ranges from 0·6 to 2·4 nm (Cadene et al., 2005). Therefore, the simulations were performed for a range of particle thicknesses including 0·6, 0·8, 1, 1·2 and 2 nm. Fig. 5 shows the cohesive stress calculated by the proposed formulation, equation (15), with solid lines, and the yield stress observed by experimental data with symbols. The experimental observations indicated that for solid contents below 4% (grey area in Fig. 5), the clay exhibits the properties of a sol, which is most susceptible to erosion by water flow (Abend & Lagaly, 2000;Michot et al., 2004;Neretnieks et al., 2009). The modelling results are in close agreement with the data, which provides confidence in the realism of the cohesive stress formulationa key component of the detachment model.

Critical solid volume fraction at the solid/liquid interface
The critical solid volume fraction is defined as the volume fraction at which the shear stress equals the cohesive stress. This clearly depends on the hydro-chemical conditions, so it is not a material-only dependent constant. Equating equations (13) and (15) provides a critical distance between the particles, h cr , from the equation The critical solid volume fraction, ϕ s cr , is obtained from the critical distance and equation (12) ϕ cr The critical soil volume fraction at the solid/liquid interface is calculated for sodium-rich smectite clay with parameters and properties given in Table 1 and for a number of hydro-chemical conditionswater velocity and ionic concentration. The results are presented in Fig. 6 and compared with those calculated by Neretnieks et al. (2017) for the same conditions and material properties. It can be seen that the calculated values are in very close agreement with the experiments. Small differences can be observed only for intermediate ionic concentrations, but these are in the order of previously reported ranges of critical fractions. In can be concluded that the proposed detachment model captures successfully the effects of ionic concentration and fluid velocity on the critical solid volume fraction.

Behaviour of compacted bentonite under eroding environment
A series of experimental work on erosion of compacted bentonite by Schatz et al. (2013) is used to test the accuracy of the prediction by the PD model. Erosion of 2 cm dia. discs of compacted sodium-rich bentonite is analysed when exposed to either stagnant or flowing water through a 1·0 mm thick artificial fracture (as shown in Fig. 7). Three experimental tests reported were selected for the simulations. A very low value for the ionic concentration (0·05 mM) is considered for this case (Neretnieks et al., 2017). (b) Test 2 is the expansion and subsequent erosion of the clay disc by flowing deionised water (average flow rate of 0·38 ml/min). (c) Test 3 represents the expansion and erosion of clay due to exposure to groundwater solution (0·04 g/l sodium chloride and an average flow rate of 0·09 ml/min).
The domain is 24 cm Â 24 cm square. The cell size, Δ, and the horizon, δ, are 0·12 cm and 0·36 cm, respectively. The material parameters used are provided in Table 1. The initial bentonite fraction in the domain outside the sample is set to zero. The solid fraction and height of the initial experimental sample is 0·6 and 2 cm, respectively. The solid fraction of the initial sample is assumed to be constant as the height of the sample is 20 times larger than the fracture thickness.
The left boundary (x = À12 cm) was kept at zero solid fraction and the particles would exit through the right boundary (x = 12 cm). Impervious boundary conditions were prescribed at the domain boundaries y = À12 cm and y = 12 cm. Water inlet and outlet boundary conditions include five input and output points in the domain. Constant heads are imposed at the left and right borders, creating a steady water flow from left to right in the domain. Figure 8 presents the solid volume fraction profile along the centre line at x = 0 for the three tests analysed by the PD model. The results are compared with experimental results reported by Schatz et al. (2013) for test 1 and test 3. The results of the model correlate with the experiments with regard to the extrusion distance. There are experimental data for test 2, but the calculations shown in Fig. 8(b) show a similar trend to those in Fig. 8(c). As shown in Fig. 8, (a) the solid volume fraction decreases with distance from the solid/liquid interface (as expected) and (b) the decrease is slower for stagnant water (test 1) compared to flowing water (test 2 and test 3) (which agrees with the experiments).
Figures 9 and 10 illustrate the evolution of swelling and erosion observed in the experiment and predicted by the model. It is noted that the images of the experimental tests used for comparison do not provide the evolution in the full domain. Although the simulations included the entire domain (i.e. outside detachment zone), the visual comparison is only possible for the areas where images were available.
The profiles of the solid volume fraction over time obtained by the experiments and those obtained by PD are in close agreement. Fig. 9 describes the evolution of the swelling front under the hydration process, while Fig. 10 illustrates the combined effects when the eroding solid has been co-transferred with the water. These observations are also in agreement with the results shown in Fig. 8(c). The discontinuous detachment interface is marked by a dashed line. It is shown that the modelling results are much smoother compared to the sharp interfaces that are clearly observed by the experiment (Figs 8(b), 8(c) and 10(c), 10(d)). The difference may be caused by the sedimentation of aggregate, the formation of flocs and experimental artefacts (e.g. transport of initially loosely bound material) (Neretnieks & Moreno, 2018), which are not considered in the current model. It is noted that in test 2 and test 3 there is a difference between the outflow boundary condition of the experiment (outlet boundary is effectively a narrow point) and the model (the entire outer boundary is open to flow). This difference is believed to have limited impact due to the length of the outlet from the solid/water interface.
Comparisons between the average extrusion distances into the fracture obtained by PD and experiment for test 1, test 2 and test 3 are shown in Fig. 11  the model's accuracy to account for the hydro-chemical effects on bentonite clay erosion. Theoretically, test 2 with deionised water should show larger expansion than test 3 with groundwater, because of the lower ionic concentration, which affects the repulsive double-layer forces in equation (3) and the cohesive stress in equation (15) and Fig. 1. However, the flow rate in test 2 is 4·2 times larger than that in test 3. This leads to a larger shear stress, which reduces the penetration distance, equation (13). This demonstrates how the proposed model captures the complex interplay between hydrodynamics and water chemistry to provide evolution in close agreement with the experimental data.

CONCLUSIONS
A bond-based PD formulation for clay erosion was presented. It includes descriptions of (a) the free swelling of clay as a diffusion process and (b) the transport of detached particles as an advection-dispersion process, with an important addition of a new detachment model for separation of particles from the gel.
The physical realism of the new detachment model was tested against a set of experimental data. It was shown that the model captured the impacts of ionic strength and flow rate on the detachment of a smectite-rich clay with good accuracy. The results highlighted that increasing the ionic concentrations and velocities of the eroding fluid led to increasing the critical clay volume fraction at the solid/liquid interface.
For validation, the full erosion model was applied to cases with available experimental data for erosion of compacted bentonite under three hydro-chemical conditions. First, it was shown that the diffusion model was appropriate for representing the free swelling, as the simulation results were in close agreement with the experimental data of test 1. Second, by comparisons with test 2 and test 3 it was shown that the full erosion model accounted realistically for the hydro-chemical conditions of the eroding fluid.  In summary, it has been demonstrated that the proposed non-local PD-based framework is efficient and effective for analysis of the hydro-chemical effect on clay erosion. Specifically, the model is able to capture the complexity of the problem arising from highly non-linear physics (input parameters vary over several orders of magnitude) and geometry (moving gel-fluid interface), as well as the discontinuous process of particle detachment. Noting the PD applications to a range of problems with discontinuity, large deformation, non-linearity and heterogeneity in solid mechanics, the existing applications to problems in soil mechanics are very limited and mainly related to granular materials. The model presented here is the first PD model related to clay behaviour, for which the effects of coupled phenomena on the overall behaviour are inherently stronger than on the behaviour of coarse-grained soils.