Modelling the piping-assisted erosion of clay barriers

Mass loss of bentonite-based barriers can occur due to erosion caused by hydro-chemical interactions (Baik et al., 2007; Reid et al., 2015; Alonso et al., 2018; Bouby et al., 2020; Komine, 2020), thus affecting the barriers’ hydraulic performance and reducing their sealing capacity. For example, in the context of a bentonite buffer in the geological disposal of high-level nuclear waste (HLW) in crystalline rock, the clay buffer can be eroded at the interface between the buffer and the host rock, resulting in gradual loss of sealing capacity and co-transport of radionuclides into the biosphere (Sane et al., 2013; Schatz et al., 2013; Missana et al., 2018; Bouby et al., 2020). Erosion induced by piping has also been reported to occur in the early phase of saturation of compacted bentonite when the groundwater flows into a deposition hole from fractures in the bedrock, causing sub-vertical mass transport into the deposition tunnel backfill (Borgesson & Sandén, 2006; Sane et al., 2013; Suzuki et al., 2013; Navarro et al., 2016; Sandén et al., 2017). Bentonite is also a potential material for plugging and sealing boreholes in the geological disposal of radioactive waste (Nagra, 2002; Posiva, 2006; SKB, 2010; Zeng et al., 2019; Kale and Ravi, 2021). Bentonite plugs have been suggested as an alternative material to cement in the abandonment of hydrocarbon wells (Towler et al., 2020; Aslani et al., 2022). The strong ability to swell and self-heal and the low hydraulic conductivity of bentonite plugs are the key characteristics that have attracted attention to their use in sealing deep boreholes and wells.

The present work focuses on the potential erosion of bentonite plugs for sealing boreholes or hydrocarbon wells induced by piping flow that occurs in the annular gap between the bentonite blocks and rock. Examples of such erosion may occur during the emplacement of the plug (lowering the plug to the required depth) (Sandén et al., 2017) and longer-term erosion that occurs when the plug is in its place and interacts with percolating water in the fractures. The authors present the first mechanistic modelling of the erosion of bentonite plugs assisted by piping, which accounts for the complex physical processes of deformation, damage and particle detachments.

Mathematical descriptions of erosion include a combination of formulations for the continuum deformation (swelling due to re-saturation of bentonite) and the discontinuous damage and erosion processes. Most research studies on the computational modelling of erosion are related to non-swelling materials, and their direct applications for the analysis of swelling clay erosion (including bentonite) have critical limitations (Sane et al., 2013). For example, studies on mechanically induced erosion have focused on the piping failure of earth dams or similar structures where the pressure differences and flow rates are substantial – for example, larger than 1 m³/s (Borrelli et al., 2011). In contrast, the groundwater flow rate in rock fractures is lower – for example, 10⁻¹² to 10⁻¹⁸ m³/s, and the hydraulic conductivity is in the range from 10⁻⁵ to 10⁻⁸ m/s (Suzuki et al., 2013). The erosion behaviour of swelling clays is complex as it involves clay–water interactions, large deformation, phase change (solid to gel or dispersed colloidal particles) and evolution of the piping channel (widening or narrowing due to combinations of erosion and swelling) (Sane et al., 2013; Yan et al., 2021).

Moreover, the main erosion mechanisms for swelling clays (e.g. bentonite) are fundamentally different from those formulated for non-swelling materials and coarse-grained soils. The initially unsaturated bentonite, once hydrated, expands to fill the gaps or voids, leading to a reduction in permeability. However, water pressure acting on bentonite may increase if the inflow is localised in fractures leading to the continuous wetting of the bentonite buffer (Borgesson & Sandén, 2006; Sandén et al., 2008). High hydraulic pressure causes a series of hydraulic phenomena, such as erosion and piping (Suzuki et al., 2013). Piping damage and erosion have been observed when the water pressure is larger than the hydraulic resistance of the bentonite buffer (Chen et al., 2016).

Although still limited in numbers, experimental studies at the laboratory scale and underground research laboratories provide evidence for potential erosion and piping of the clay buffer during the re-saturation process. The experimental observations from the SKB and Posiva's BACLO projects (Sweden and Finland), the LOT tests at Äspö URL (Sweden), EPSRC's SAFE barriers (UK), the BELBaR project (EU) and the in situ tests at Horonobe URL (Japan) describe complex interactions governing the erosion and piping of clay barriers. Experimental investigations have primarily focused on understanding: (i) the conditions for piping formation, (ii) the evolution of piping channel, (iii) the effects of inflow rate, and (iv) the effect of piping on buffer properties and eroded mass (Sandén et al., 2008; Suzuki et al., 2013; Jo et al., 2019). In situ tests indicated that the eroded mass of the bentonite buffer induced by piping could reach several kilogrammes (Sandén et al., 2008). In addition, water flowing along the bentonite plug was reported to cause a large amount of mass loss in a short time (43·5% mass loss in 1 h) and to compromise the sealing performance of bentonite plug boreholes (Sandén et al., 2017).

A limited number of theoretical studies of piping-assisted erosion of swelling clay have been undertaken recently. Navarro et al. (2016) presented a coupled swelling and mechanical erosion model for compacted MX-80 bentonite. The mass loss calculation was evaluated using a simplified experimentally calibrated erosion model. The model was extended by Asensio et al. (2018) to account for the effects of water salinity. In previous works, the authors of this paper have introduced a new mechanistic modelling approach for analysing the erosion of swelling clays (Sedighi et al., 2021; Yan et al., 2021) based on the theory of peridynamics (PD) and by looking into forces acting between clay particles. These earlier works have been focused on quantifying the erosion of a clay buffer in artificial fracture systems where the clay buffer can swell and penetrate several metres into the fracture (e.g. maximum 10 m for 1 m/year water velocity) (Moreno et al., 2011; Huber et al., 2021).

This paper presents the mathematical approach explicitly developed to investigate the erosion of swelling clays assisted by piping flow. This could be especially important in HLW backfill and plugging of boreholes/wells where the expansion of swelling clays is limited to piping channels, which generally have sizes on the order of a few millimetres/centimetres (Sandén et al., 2008; Sane et al., 2013; Suzuki et al., 2013; Navarro et al., 2016). In addition, the piping channels/gaps may be widened if the erosion potential induced by piping flow is beyond the self-healing capacity of swelling clays; contrarily, the piping channels/gaps can be closed/filled due to the swelling clays’ strong self-healing characteristics (where hydro-chemical interactions of clay with groundwater do not hinder swelling). The work presented herein is a step toward understanding the role of piping in the erosion behaviour of swelling clays. However, the initiation and propagation of piping channels are not included in the current model. The model incorporating the fundamental physical interactions controlling the erosion, with PD formulations of swelling, detachment of diluted particles and detached particles transport processes, is presented together with model validation and analysis of erosion behaviour in piping channels. The model focuses on the characterisation of the erosion process, without integrating in the formulation the hydromechanical behaviour of the material that is eroded.

The key processes that occur during clay erosion caused by piping flow include: (i) clay swelling into the existing piping void; (ii) detachment of clay particles at the clay/water interface; and (iii) co-transport of the detached particles by the piping flow. These processes labelled (i) to (iii) are shown in Fig. 1.

The continuum and PD formulations for clay erosion in fractured rock (Sedighi et al., 2021; Yan et al., 2021) are extended here to the case of clay erosion assisted by piping. The main difference between the clay erosion in fractured rock and the clay erosion assisted by piping is in the water velocity; typical flow rates in a piping channel range between 10⁻³ and 1 m/s, much larger than typical flow rates in a fractured system/rock, which range between 10⁻⁸ and 10⁻⁵ m/s. Therefore, the mass loss by erosion in a fractured system combines chemical erosion (extruded mass) and mechanical erosion (particle detachment at the interface). In contrast, bentonite erosion by piping flow is dominated by mechanical erosion in the piping channels, which evolves during the erosion process; thus, the channels are either widening or narrowing depending on the balance between erosion and swelling. The overall behaviour is governed by the rate of swelling and the detachment rate, which vary based on the material properties (hydraulic properties, including the hydration rate) and flow conditions (flow rate, water chemistry).

Piping-assisted erosion is described by coupling three formulations based on the bond-based PD theory: free swelling of clay (Fig. 1(a) and (ii)), particle detachment; and detached particle transport (Fig. 1(a) and (iii)). The theory of PD considers a material domain as a collection of PD points (Fig. 1(b)). The PD points are not geometrical points but possess physical properties that depend on the problem at hand. As shown in Fig. 1(b), an arbitrary point x interacts with all points, x′, in a finite spatial region called horizon, H_x. The horizon radius is denoted by δ. The distance vector between points x and x′ is ξ = (x′ − x), and a PD bond describes the interaction between them. The horizons and bonds are shown in Fig. 1(b). Three types of PD bonds describe the three distinct interactions involved in piping-assisted erosion. These are solid–solid bonds (S–S) which control the swelling process; solid–liquid interfacial bonds (S–L) that control the detachment of particles at the clay/water interface; and liquid–liquid bonds (L–L) that control the transport of detached particles with the piping flow.

Free swelling is considered as a solid diffusion process, following the theoretical developments based on the dynamic force balance method by Liu et al. (2009). However, the solid diffusivity may change by five orders of magnitude as the clay solid content (density) changes (void ratio increases during expansion). This makes the problem highly non-linear, for which the PD is a suitable solution.

The solid diffusion process in the PD framework is given by Yan et al. (2020, 2021) and Sedighi et al. (2021) as follows

where ϕ_s(x, t) is the clay solid content; t is the time; V_x′ is the horizon volume of particle x′; and d_s(x, x′, t) is the PD microscopic diffusivity.

The relationships between PD microscopic diffusivity and measurable macroscopic diffusivity, D_s, for one-dimensional (1D) and two-dimensional (2D) problems are (Yan et al., 2020):

where D_s is the macroscopic solid diffusivity, which is calculated by

where χ is the particle's energy and f_r is the friction coefficient.

Equation (4) is obtained by developing a dynamic force balance, including the diffusion forces (F_T), attractive van der Waals forces (F_A) and repulsive electrical double layer forces (F_R) between clay particles (Sedighi et al., 2021).

The friction coefficient and particle's energy are given by Liu et al. (2009) and Neretnieks et al. (2009) 

and

where r_eq is the equivalent radius of the non-spherical particles; V_p is the volume of the particles; k₀ 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; η_w is the dynamic viscosity of water; k_B is the Boltzmann constant; T is absolute temperature; δ_p is the particle thickness; and h is the separation distance between the flat particles.

The calculations of attractive van der Waals forces (F_A), repulsive electrical double layer forces (F_R) and separation distance of particles (h) are provided in the online supplementary material.

The external force causing particle detachment is the shear force (τ_f) induced by the water flowing in the pipe channel. Detachment occurs when τ_f becomes larger than the interparticle forces that maintain interface solid particles attached to their surrounding solid particles (Laxton & Berg, 2006; Sane et al., 2013). The interactions of an interface solid particle with its surrounding solid particles provide a cohesive strength to the particle. This is related to the interaction forces with other particles described by (Sedighi et al., 2021)

With these settings, a detachment function is given by

Equation (8) is a criterion that ensures that when the cohesive strength calculated from S–S bonds is smaller than the shear force, the solid particles transform into liquid particles and follow the governing equations for transport processes (Yan et al., 2021). Incorporating the detachment equation (8) into equation (1) leads to PD formulation for coupling the free swelling and particles detachment:

The transport of detached particles by piping flow is described by the advection and dispersion equation (ADE). The PD formulation of transport of detached clay particles is given by Yan et al. (2020, 2021) and Sedighi et al. (2021) as

where d_l(x, x′, t) and v_l(x, x′, t) is the PD microscopic diffusivity of the detached particles in water and the microscopic liquid velocity, respectively. Similarly to equations (2) and (3), PD microscopic diffusivities and advection for 1D and 2D cases are given by

and

where D_l is dispersion coefficient, which is the sum of diffusion (D_d) and mechanical dispersion (D_m) coefficients of detached particles in water as follows:

where α_L is the dispersivity parameter.

The liquid velocity (V_l) 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 soil viscosity.

A sequential approach is adopted to solve the equations for free swelling and particles detachment – equation (9) – and detached particles transport processes – equation (10). The domain of interest is discretised into subdomains using uniform linear and square sub-grids (length in one dimension and area in two dimensions), respectively. The Euler method is used for time integration. The discretisation of domain and integration of time are presented in Sedighi et al. (2021) and Yan et al. (2021). The swelling equation is solved first for each time step where the clay solid volume fraction profile is obtained to compute the cohesive strength of particles. The detachment interface is, therefore, automatically updated by equating the cohesive strengths and shear forces. The transport of detached particles is then calculated. The eroded mass is obtained by summating detached particles after each time step. Details of numerical implementation and computations can be found in Yan et al. (2021) and Sedighi et al. (2021).

The model is tested against the results of a pinhole experiment reported by Sane et al. (2013) (shown in Fig. 2). This exercise allows for validating the coupled swelling and erosion model when an MX-80 Wy-bentonite sample is exposed to hydration and erosion by piping flow. Referring to the experimental set-up shown in Fig. 2, the liquid flow was through a pinhole at the centre of cylindrical samples (100 mm dia. and 100 mm long). The constant flow rate was 0·1 l/min, representing the potential erosion of compacted bentonite buffer at the possible inflow rate in a deposition hole for the case of a repository in Finland (Juvankoski et al., 2012). The liquid used had a low ionic concentration (e.g. 0·58 mM). The initial dry density of MX-80 Wy-bentonite was 1700 kg/m³. A steel rod with a smooth surface and an initial diameter d_h (6 or 12 mm) was inserted into the cell (Navarro et al., 2016). An automated effluent collector collected effluent samples. The initial and boundary conditions applied in the simulations of the pinhole experiment are also shown in Fig. 2. The shear stress induced by the flowing water on bentonite particles is obtained by assuming that a quasi-steady laminar flow has fully developed in the central hole of the sample (Navarro et al., 2016). The hydraulic shear stress, therefore, is calculated by (Navarro et al., 2016):

where η_w is the dynamic viscosity of the water; Q is the water flow rate (0·1 l/min); and r_h is the characteristic mean radius of the central hole, which is a function of time and can vary along the cylinder due to the surface heterogeneities of the sample.

In the absence of information about the surface heterogeneities of the sample, it is assumed that the central hole has a constant radius r_h (see Fig. 2(a); Navarro et al., 2016). However, during the simulations, the radius of the piping channel may change due to the combined action of swelling and erosion. The simulation results of the bentonite solid volume fraction with time across the width of the sample are provided in the online supplementary material (see Fig. S1). The shear stress calculated from equation (17) for r_h = 3 mm and r_h = 6 mm is 0·078 Pa and 0·00975 Pa, respectively. The experimental data for r_h = 3 mm are labelled as s1a, s1b and s1c, and for the r_h = 6 mm case they are labelled as s3a in the following discussions and comparisons with numerical results.

The material properties and parameters used for the simulations are adopted from the literature. They are given in Table 1. The cell and horizon sizes are denoted by Δ = 0·1 mm and δ = 0·3 mm, respectively.

Navarro et al. (2016) presented an erosion model to describe mass loss at the water/clay interface. The mass loss was obtained using an erosion rate (k_e) from experimental data. The erosion rate for r_h = 3 mm and r_h = 6 mm was reported to be 0·004 s/m. Fig. 3 compares the cumulative mass loss per unit length from experimental data, including the modelling results by Navarro et al. (2016) and the PD model presented in this paper. The bentonite mass loss rate (N_erosion) from the detachment at the clay/water interface in the PD simulations was calculated by (Moreno et al., 2011; Neretnieks et al., 2017):

where ϕ_s,R is the critical clay solid volume at the detachment interface, which can be obtained by equating the cohesive stress (equation (7)) and the shear stress (equation (17)); L is the length of the piping channel; r_R is the radius of the piping hole; and D_R is the diffusion coefficient for bentonite particles at the clay and liquid interface (e.g. r = d_t, see Fig. 2(b)). A value of D_R =2 × 10⁻¹⁰ m²/s was used for the response to ionic solutions at low concentrations (e.g. <1 mM) case.

Figures 3(a) and 3(b) show a comparison between the eroded mass obtained experimentally, the simulation results by Navarro et al. (2016) (dashed black line), and the simulation results from the PD model (solid black line). It can be seen that the eroded mass obtained with the present model is in good agreement with the experimentally measured value in two of the experiments with similar records (s1a and s1b). All pinhole experimental results exhibit a level of dispersion (Sane et al., 2013; Navarro et al., 2016). It is noted that the erosion behaviour of three MX-80 batches was different, as seen in Fig. 4(a). The erosion rates measured by Exp.s1a and Exp.s1c differ by a factor of 2. The experimental datasets Exp.s1a, Exp.s1b and Exp.s1c are, in principle, identical. However, the recorded data appear to vary substantially due to the limitations of the erosion test equipment and the heterogeneity of the clay samples. The predictions reported by Navarro et al. (2016) and those from the present study are in close agreement. The main difference is that the results by Navarro et al. (2016) were obtained after calibrating the erosion rate in their model against experimental data, while with the current model, the cumulative mass loss was calculated without any calibration – which is an advantage of the coupling of the swelling and detachment formulations. Specifically, the erosion process is controlled by the balance between the cohesive strength (which depends on the swelling) and shear force (which depends on water velocity and chemistry). For the case with r_h = 6 mm, the value predicted by the model mass loss is larger than the experimentally measured value between 10 h and 65 h. This over-prediction is likely to be due to the assumption of a constant radius of the pipe (r_h). In addition, Navarro et al. (2016) pointed out that local effects at the top and bottom sample boundaries may influence the erosion behaviour of bentonite for shorter samples in the pinhole test. The boundary effects on the erosion behaviour are not included in the current model. However, the model is seen to predict the trend correctly and potentially can give the mass loss at longer times. In all cases, even when the experimental results exhibit a non-negligible dispersion, the model correctly reproduces the trend of the mass erosion rate. It is noted that no artificially created channels are expected in clay buffer and backfill in the context of geological disposal applications. Therefore, understanding and modelling the initiation and propagation of pipes is important for the performance assessment of bentonite barriers. Navarro et al. (2016) reported that the naturally created piping channels might close, while artificially created channels did not close. Considering that the results of this work are obtained with an artificially created channel, it is expected that the behaviour of a bentonite buffer in the field may be different.

Bentonite in a compacted form has been proposed for plugging and sealing boreholes (Pusch & Ramqvist, 2004, 2008; Borgesson & Sandén, 2006; SKB, 2010; Arnold et al., 2011). Also, considerations are given to compacted bentonite as an alternative to cement for sealing off hydrocarbon wells (Towler et al., 2008, 2020; Aslani et al., 2022). During the installation of the bentonite plug in the borehole, water passes at the annular gap between the plug and the rock at relatively high velocities (e.g. 1 l/s, depending on the lowering rate of the bentonite plug during installation). Such a rate of fluid flow can potentially cause erosion. Understanding the potential swelling and erosion of the plug is, therefore, important to ensure that the damage during installation is minimal (the installation rate is not too large) and swelling is kept limited (the installation rate is not too small). Erosion could also occur when the plug is fully installed at the gap between the plug and the rock (damaged or undamaged plug). The current authors present a set of simulations using the proposed model and compare the results with experimental observations of a similar problem by Sandén et al. (2017). Sandén et al. (2017) reported a series of experiments on bentonite plugs that were eroded in the tests by flushing water at a 3 mm gap between cylindrical samples of compacted bentonite and the experimental cell. Fig. 4 shows the bentonite block prior to installation in the cell (Fig. 4(a)) and the eroded sample after the test (Fig. 4(b)).

Cylindrical MX-80 bentonite blocks with a diameter of 80 mm and length of 0·25 m were used in this study. The gap between the wall and the block was d_g = 3 mm. The bentonite block was flushed with tap water at a 1 l/s rate for 1 h. The shear stress at the water/clay interface induced by flowing water was approximately 0·5 Pa (τ_f ≈ ηV_l/d_g; Eriksson & Schatz, 2015). For bentonite plugs, a piping gap exists between the plug and the rock, hence separate consideration of initiation is not necessary. The parameters used for the simulation are provided in Table 2. A clear visible erosion profile can be observed in Fig. 4(b).

Figure 5 compares the eroded mass obtained by the experiments and the numerical simulations. The present model reproduced the trend observed in the experiments. The model with the assumption of homogeneity (see black solid line in Fig. 6) underestimated the eroded mass amount at the end of the test (t > 50 min). It can be observed from Fig. 4(b) that the erosion profile of the bentonite is non-uniformly distributed, especially at the bottom of the sample. The radius of the lower bentonite block was much smaller after erosion than that of the upper block, as the water was injected from the bottom of the sample. Previous studies (Liu et al., 2009; Neretnieks et al., 2009, 2017) demonstrated that the swelling and erosion of bentonite is sensitive to the variation of the clay plate thickness (δ_p), which has a range from 0·6 nm to 2·4 nm (Cadene et al., 2005; Liu et al., 2009; Yan et al., 2021). It makes a difference to the strength of repulsive forces (F_R) and the critical clay solid volume (ϕ_s,R). The clay heterogeneity was consequently accounted for by considering that the thickness of the bentonite plate (δ_p) followed a Weibull distribution with shape parameter 5 (Tang et al., 2016; Wang & Wang, 2022). Taking into account the variability of the bentonite plate, the model demonstrated a better performance in the prediction of the eroded mass at the end of the test (see dashed line in Fig. 5). The amount of bentonite mass loss at the end of the test was 10·3% of the original mass at installation (see Table 2). The large amount of bentonite mass loss induced by erosion can be detrimental to the sealing capacity of bentonite plugs in boreholes. The possible loss of material should be carefully investigated at the planning and design stages using the above analysis approach.

The erosion can considerably reduce the clay's initial density during the installation of the plug in the borehole. It is assumed that the clay's maximum allowable mass loss by erosion is 10% of the original mass at installations with lengths up to 1000 m (Pusch & Ramqvist, 2004). Studies on the performance and quality assessment of boreholes have emphasised the effects of flow velocity and diameter of boreholes on erosion during the emplacements of the plug. Bentonite seals developed by SKB are primarily intended to be used in boreholes ranging in diameter from 56 mm to 100 mm and in boreholes of up to 1000 m depth. However, a larger diameter (e.g. d_h = 160 mm) was used in boreholes because more flexibility is given for all operations in the framework of the plugging phase (Chaplow, 2011). The measured inflows in the field range between 10⁻⁴ l/s and 1 l/s (Winberg et al., 2000). To assess the erosion of bentonite plugs with different diameters and flow rates, a set of simulations are carried out based on the geometry and properties of a bentonite plug described by Sandén et al. (2017). MX-80 bentonite with an initial dry density of 1787 kg/m³ is used as a plug material for sealing deep boreholes (e.g. 716 m) following Sandén et al. (2017). The gap between the wall and the block is assumed to be d_g = 3 mm. The time required to finish the plug emplacements depends on the flow rate passing the annular gap between the rock and the bentonite. For example, a flow rate of 1 l/s during 1 h corresponds to lowering bentonite blocks into a borehole with a diameter of 80 mm to a depth of 716 m (Sandén et al., 2017). Decreasing the flow rate of 1 l/s by a factor of 10 can increase the time for erosion from 1 h to 10 h. The parameters used for the simulations, in this case, are presented in Table 2.

Figures 6 and 7 show the mass loss and the erosion rate of the plug considering different borehole diameters (e.g. d_b = 56 mm, 80 mm and 160 mm) and as a function of flow rate (Q). It can be observed that the erosion rate increases with the flow rate, as larger flow rates cause higher shear stress at the clay and water interface. This is compatible with the trends reported by Schatz et al. (2013) and Yan et al. (2021). It is found that the differences observed for the erosion rates of the boreholes with different diameters (d_b) are minor, while d_b shows a significant impact on the final eroded mass after the emplacements of the plug. For example, the eroded mass for d_b = 56 mm is 3·1 times lower than that for d_b = 160 mm at Q = 1 l/s. The observed increase in eroded mass for a larger diameter is due to the increase of erosion surface areas, see equation (18). In addition, Fig. 6 shows that the amount of eroded mass can be reduced by increasing the flow rate. This is because larger flow rates require a shorter time to finish the emplacements, although higher erosion rates exist (see Fig. 6).

A bentonite plug's sealing performance depends on its final dry density. Fig. 7 shows the mass loss fraction (from the original sample) at installation after the emplacements of the plug, which is calculated by

where, m_o and m_e are the original sample's initial and eroded mass.

Based on the results presented in Fig. 7, increasing the flow rate and the diameter can reduce the mass loss fraction. For example, the eroded material increased from 3·37% for d_b = 160 mm and Q = 1 l/s to 43·9% for d_b = 56 mm and Q = 0·01 l/s. This is because the initial mass of the sample for a larger diameter is bigger, which decreases the mass loss fraction, see equation (19). This outcome agrees with the findings of Pusch et al. (2016),  which showed that the final density and tightness of the clay plugs increased significantly with increased borehole diameter. If the maximum allowable erosion mass is set as 10% of the original mass, the minimum flow rate required is 0·8 l/s, 0·5 l/s and 0·045 l/s for d_b = 56 mm, 80 mm, and 160 mm, respectively. The results in Figs 6 and 7 show one potential application of the proposed model to reduce the uncertainties related to erosion in the performance of the clay plug in the boreholes.

Assessing clay performance as buffer/backfill or sealing material requires robust predictive models that include expected environmental conditions and multi-physics phenomena involved in the erosion process. This paper presented a non-local model for coupling the swelling and erosion of clay in piping flow. The model was validated by comparing its predictions to a series of experimental tests, including free swelling tests and erosion assisted by piping flow.

The results demonstrated that the coupled swelling and erosion model quantified the amount of mass loss of swelling clay under piping flow. The model correctly reproduced the trend observed in the experiments, including the pinhole test and erosion of plug material in the borehole. Compared to the model presented by Navarro et al. (2016) using an experimental fitted erosion rate, the cumulative mass loss obtained by the new model is calculated directly by establishing the relationship between the internal force among the bentonite particles and the shear force induced by flowing water. The amount of mass loss of plug material at the end of the test was 10·3% of the installed mass. The significant mass loss induced by erosion is expected to reduce the sealing capacity of bentonite plugs in boreholes. In addition, the diameter of plugs shows a substantial impact on the final eroded mass after the emplacements of the plug. For example, the eroded mass for d_b = 56 mm is 3·1 times lower than that for d_b = 160 mm at Q = 1 l/s. The observed increase in eroded mass for a larger diameter is attributed to the increased eroding surface areas caused by a bigger diameter of clay plugs. In addition, increasing the flow rate can reduce the amount of eroded mass. This is because larger flow rates require a shorter time for finishing the emplacements.

H. Yan acknowledges the financial support through a joint PhD scholarship by the China Scholarship Council (CSC no. 201808350074) and the Department of Mechanical, Aerospace and Civil Engineering at the University of Manchester. M. Sedighi acknowledges the support of the Royal Society, UK, by way of IEC\NSFC\181466. A. P. Jivkov acknowledges the support of the Engineering and Physical Sciences Research Council (EPSRC), UK, by way of grant EP/N026136/1.

A_H

Hamaker constant

a_p^

specific surface area

D_d

molecular diffusion

D_l

dispersion coefficient

D_m

mechanical dispersion

D_p

particle diameter

D_R

diffusion coefficient for bentonite particles

D_s

macroscopic solid diffusivity

d_b

diameter of boreholes

d_l(x, x′, t)

peridynamic (PD) microscopic diffusivity of the detached particles

d_s(x, x′, t)

PD microscopic diffusivity

F

Faraday constant

F_A

attractive van der Waals forces

F_R

repulsive electrical double layer forces

F_T

diffusion forces

f_h

mass bond

f_r

friction coefficient

H_x

horizon of material point x

h

separation distance

k₀

pore shape factor

k₀τ²

Kozeny constant

k_B

Boltzmann constant

k_e

erosion rate

L

length of the piping channel

m_o

eroded mass

m_e

initial mass

N_erosion

mass loss rate

p

fluid pore pressure

Q

water flow rate

R

gas constant

r_eq

equivalent radius of the non-spherical particles

r_h

characteristic mean radius

r_R

radius of the piping hole

S_p

particle surface area

T

temperature

T_w

fracture transmissivity

t

time

V_l

liquid velocity

V_p

volume of the particles

V_x′

horizon volume of particle x′

v_l(x, x′, t)

PD microscopic liquid velocity

α_L

dispersivity

δ

radius of the horizon

δ_p

particle thickness

ε₀

permittivity of vacuum

ε_e

mass loss fraction

ε_r

relative permittivity of water

η_s

soil viscosity

η_w

viscosity of water

μ(x, x′, t)

detachment function

ξ

distance vector between the two material points

σ⁰

surface charge of particles

τ

tortuosity

τ_c

cohesive force

τ_f

shear force

ϕ_s,R

critical clay solid volume at the detachment interface

ϕ_s(x, t)

clay volume fraction at points x

ϕ_s(x′, t)

clay volume fraction at points x′

$ϕ¯sx,x′,t$

average clay volume fraction

χ

particle's energy