Three-dimensional fully coupled hydro-mechanical-chemical model for solute transport under mechanical and osmotic loading conditions

Mechanical deformation and chemico-osmotic consolidation of clay liners can change its intrinsic transport properties in all direction and can alter fluid and solute transport processes in the entire model domain. These phenomena are described inadequately by lower-dimensional models. Based on the Biot’s consolidation theory, fluid and solute mass conservation equations, a three-dimensional (3D) fully-coupled hydro-mechanical-chemical (HMC) model has been proposed in this study. The impacts of mechanical consolidation and chemico-osmotic consolidation on permeability, hydrodynamic dispersion, solute sorption, membrane efficiency, and chemical osmosis are considered in the model. The model is applied to evaluate performances of a single compacted clay liner (CCL) and a damaged geomembrane-compacted clay composite liner (GMB/CCL) to contain a generic landfill contaminant. Effect of model dimensionality on solute spread for CCL is found to be marginal, but for GMB/CCL the effect is significantly large. After 50-year simulation period, solute concentration at the half-length of the GMB/CCL liner is predicted to be 40% of the source concentration during 1D simulation, which is only 6% during the 3D simulation. The results revealed approximately 74% over-estimation of liner settlement in 1D simulation than that of the 3D for GMB/CL system. Solute spread accelerates (over-estimates) vertically than horizontally since overburden load and consequent mechanical loading-induced solute convection occurs in the same direction. However, in homogeneous and isotropic soils, horizontal spread retards the overall migration of contaminants, and it highlights the importance of 3D models to study solute transports under mechanical and chemico-osmotic loading conditions in semi-permeable clays, especially, for damaged geomembrane-clay liners. The results show the utility of geomembranes to reduce soil settlement, undulation, and restriction of solute migration. Furthermore, application of geomembrane can inhibit development of elevated negative excess pore water pressure at deeper portion of a clay liner.


Introduction
Clays and clayey soils are often used as buffer materials, such as landfill liners to protect the surrounding geology and groundwater around landfill sites. Detailed understanding of physical, chemical, and mechanical behavior of clay liners are crucial to ensure their effective and adequate performance (Chen 2014;Reddy et al. 2017). Numerous research efforts and investigations are available in literatures focusing on the material's behavior under various physical, chemical, and mechanical loading conditions Yan et al. 2021aYan et al. , b, c, 2022. However, the current work is focused on understanding the osmotic behavior of clays and its impact on solute transport and settlement behaviors of landfill clay liners. This is of significant importance where landfill leachate is in direct contact with the landfill liners (Zhang et al. 2018;. Studies of contaminant (or solute) transport in clay liners, usually assume that the liners are rigid porous media, and flow of solutes are driven by advection and/or diffusion mechanisms and influenced by sorption processes (Li et al. 2002;Islam and Singhal 2002;Bourg et al. 2003;Kooi et al. 2003;Dominijanni and Manassero 2005;Zhang et al. 2005;Reddy et al. 2017;Qiu et al. 2022). However, clay liners undergo mechanical consolidation due to the weight of landfill wastes altering its intrinsic transport properties (Huang et al. 2012;Zhang et al. 2018;Zheng et al. 2021). Influences of liner consolidation on solute transport have been extensively studied in Smith, (2000); Alshawabkeh and Rahbar (2006); Zhang et al., (2013;Xie et al., (2016); Yu et al., (2018); Pu et al., (2020); Yan et al., (2021aYan et al., ( , b, c, 2022; Qiu et al., (2022). Meanwhile, clayey soils that possess semi-permeable membrane properties promote osmosis, a process where movement of a solute occurs from higher to lower concentration region, while the solvent moves in the opposite direction in the presence of the solute concentration gradient across the semi-permeable membranes. The movement of pore fluid drives osmotic-consolidation which has been of particular interest over the past few decades. Initial theoretical developments of Greenberg et al. (1973), Barbour and Fredlund (1989), and Mitchell (1993) on chemical osmosis and osmotic-consolidation are extended and advanced in the works of Kaczmarek and Hueckel (1998), Smith (2000), Smith (2002, 2004), Malusis and Shackelford (2002), Lewis et al. (2009), Huang et al. (2012), Musso et al. (2013), Pu et al. (2020), Xie et al. (2016) and others. Zhang et al. (2018) conducted a comprehensive analysis of those models and revealed that majority of those models did not consider multiphase flow and deformation and their simultaneous interactions on osmotic processes. In other words, the models are mostly un-coupled or partially coupled in nature, which is a significant limitation, since clay liners in landfills are exposed to multiphase flows of fluids, chemicals, heat as well as mechanical deformation conditions. They are also largely limited to one-dimensional or two-dimensional problems and do not represent realistic or in situ design layouts of landfill liners. For example, mechanical deformation of clay liners can change its intrinsic transport properties in all direction, which eventually can alter fluid and solute transport processes in the entire model domain. These phenomena are described inadequately by lower-dimensional models. Within the scope of their work, Zhang et al. (2018) developed a fully coupled HMC model by addressing the limitations of the previous models, e.g., Kaczmarek and Hueckel (1998), Smith (2002, 2004), Malusis and Shackelford (2002), Musso et al. (2013), Pu et al. (2020). Nevertheless, the model of Zhang et al. (2018) was also one-dimensional. Huang et al. (2012) and Zhang and Fang (2016), Zhang et al. (2017) developed 3D coupled models describing solute transport in deforming soil taking account of the effect of consolidation on solute transport processes. However, to the authors' knowledge, a fully coupled, three-dimensional HMC model that investigates both mechanical consolidation and osmotic-consolidation is rarely available. Accurate assessment of clay liner performance is crucial to design landfill barriers for optimum containment of chemical solutes (Chen 2014). Variation of spatial dimension can influence the migration law of solutes (Li et al. 2011). Therefore, it is essential to establish a three-dimensional model to investigate solute transport processes under complex engineering and geoenvironmental conditions.
In this work, the previously published model of Zhang et al. (2018) is further extended to develop a fully coupled 3D HMC model with the aim to overcome the above limitations and accurately represent landfill clay liners. Both mechanical consolidation, associated with landfill waste, and chemicoosmotic consolidation, driven by solute concentration gradient, are included in the model formulations. Impacts of dynamic hydraulic permeability, hydrodynamic dispersion, solute sorption, membrane efficiency, and chemical osmosis on solute transport behavior are considered in the proposed model. The finite element software COMSOL Multiphysics is used to develop the 3D numerical model. The model is applied to evaluate performances of a single compacted clay liner and a geomembrane-compacted clay liner (GMB/CCL) to contain a generic landfill contaminant (or solute).

Coupling model and governing equation
In this section, the governing equations of the fully coupled 3D HMC model is presented. The following assumptions are made: soil is homogeneous, isotropic, and under isothermal condition; soil particles are incompressible and deformation is linear and elastic.

Governing equation of soil deformation
Given that soil deformation is marginal, the stress equilibrium equations follow: where ij is stress (i, j = x, y, z); is unit weight of pore water. According to the principle of effective stress (Terzaghi 1943), the stress equilibrium can be written as: where ′ ij is the effective stress. For isotropic linear poroelastic medium, the constitutive relations can be written in terms of effective stress σ ij ′, strain ε ij , and chemical (1 (Li et al. 2011;Zhang and Fang 2016) follows: where α c = m c (1 − )/( 1 + ), and m c is the coefficient of volume change due to chemical concentration change; G is shear modulus of soil medium; is Poisson's ratio; E is elastic modulus of soil medium; c s is the concentration of contaminant in the external environment; c is the concentration of contaminant in soil; v is volume strain; ij is the strain. In addition, the strain-displacement relations can be written as follows: where w i is the displacement in the i-direction (i = x, y, z). Substituting Eqs.

Governing equation of pore fluid flow
The flow continuity equation can be expressed as follows: where ρ f is the density of the pore fluid; n is porosity and q f is the absolute velocity of the pore fluid. By considering that the pore fluid is incompressible, Eq.(6) yields: Osmotic potential driven by solute concentration gradient across semi-permeable clay liner may lead to soil deformation and change in porosity. Following Kaczmarek and Hueckel (1998), Musso et al. (2013), and Zhang et al. (2018), the change of porosity can be expressed as follows: where m v is the coefficient of volume change due to the effective stress, m and is the coefficient of volume change due to osmotic potential. Substituting Eq. (8) into Eq. (7) yields: The absolute velocity of the pore fluid q f can be written as follows: where v s is the velocity of soil skeleton and can be written as follows: where w is the displacement tensor and can be written as w = (w x , w y , w z ) T . q r is the pore fluid velocity relative to the soil skeleton and it can be expressed as follows: where q is the Darcy velocity of the pore fluid and can be written as q = (q x , q y , q z ) T . Zhang et al. (2018) considered the influence mechanical loading and osmotic potential on q as follows: where ω is osmotic efficiency; R is universal gas constant; T is the absolute temperature, k x , k y , and k z are the permeability coefficients in x, y, and z directions. When the permeability of soil is isotropic, average permeability k = k x = k y = k z . Hart and John (1986) expressed in situ hydraulic conductivity in relation to the initial hydraulic conductivity and initial porosity as follows: where k 0 is the initial hydraulic conductivity and n 0 is initial porosity.

Governing equation of solute transport
Following the principle of mass conservation, solute transport in pore fluid can be written as follows: where Y is the sink/source; J f represents the solute flux tensor in the pore fluid. The term J f can be expressed as follows: in which J a is the solute convective flux tensor; J D is the solute hydrodynamic dispersion flux tensor; J m is the solute chemical diffusion flux tensor; J d is the solute mechanical dispersion flux tensor.
The solute convective flux tensor is given by the following: The solute chemical diffusion flux tensor can be expressed as follows: where D m represents the effective diffusion coefficient tensor of contaminant in porous medium and can be defined as (Malusis et al. 2014) follows: where D 0 is the diffusion coefficient for the contaminant in a free solution; τ is the tortuosity factor for soil and can be determined by the empirical formula, τ = (τ x , τ y , τ z ) T = (n m , n m , n m ) T ; m is the empirical parameter (Liu et al. 2004).
The solute mechanical dispersion flux can be written as follows: where D d is the mechanical dispersion coefficient tensor and can be expressed as follows: where α L is the dispersivity tensor and can be written as α L = (α T , α T , α L ) T , α T is transverse dispersion, and α L is longitudinal dispersion. In practice, the combination of molecular diffusion and mechanical dispersion is called hydrodynamic dispersion, and the hydrodynamic dispersion coefficient tensor can be written as follows: where D is the hydrodynamic dispersion coefficient tensor. The solute hydrodynamic dispersion flux tensor can be expressed as follows: Substituting Eqs. (18), (19), and (25) into Eq. (17) yields: Similar to solute transport in pore fluid, the mass conservation equation of a solute in the solid phase can be written as follows: where J s is the solute flux tensor in soil particles and S is the mass of contaminant adsorbed within the soil particles per unit mass of solid. Considering the adsorption isotherm is liner: For convenient calculation, the linear isothermal adsorption has been occupied as follows: The flux J s is described by the following: Substituting Eqs. (28)-(29) into Eq. (27) yields: Assume that the sink/source in the sorption process is the same as that in the desorption process. Equations (26) and (30) can be combined to achieve the governing equation of solute transport as follows: Equations (5), (16), and (31) constitute the threedimensional, fully coupled HMC model equations for solute transport including osmotic and mechanical loading in saturated, semi-permeable clays. The main variables in the fully coupled equation include the displacement vector w, the excess pore water pressure u, and the chemical concentration c. In addition, the use of the coupling parameters n, q, and D affected by mechanical load and solute concentration enables the governing equation to realize the full coupling process.

Model verification/evaluation
Laboratory experimental data to validate the 3D coupled model is not currently available in literatures. Therefore, in this section, the proposed model is evaluated against the analytical solution of Yan et al. (2021a).

Model application
Landfill liner system often used to restrict/contain leachate and solutes from seeping into the surrounding geology or groundwater. Single compacted clay liner and GMB/CCL are always used as the barrier system, and their service life are generally longer than the operation period of landfill site and the stabilization time of solid waste (which could be up to 50 years). In this section, the proposed model is used to predict evolution of excess pore water pressure, soil deformation, and solute migration behavior under the combined action of mechanical loading, and chemical loading in single compacted clay liner (CCL) and geomembrane-compacted clay liner (GMB/CCL) systems. When compacted clay is used as the contaminant barrier, CCL is directly exposed to leachate. Although GMB/CCLs provide higher containment than CCL, they often experience unwanted damage during installation in the landfill sites. By analyzing field data, Rowe and Brachman (2004) concluded that 70% of the sites reported torn or ripped geomembranes, even though the construction quality was strictly controlled (Bouazza, 2002;Rowe and Brachman 2004;Brachman and Gudina 2008). Therefore, in this study, two possible scenarios are investigated where a clay liner is exposed directly to the leachate (Case 1) or through a damaged portion of a GMB/CCL (Case 2). A conceptual physical system is illustrated in Fig. 2a, and the idealized model is presented in Fig. 2b and d for CCL (Case 1) and damaged or leaked GMB/CCL (Case 2), respectively.
The model domain is 4 m × 4 m × 1 m (Ω 1 ) which contains a circular leakage (Ω 2 ) of 0.1 m radius at the center of the top   Fig. 1 Comparison of the model predicted solute concentration results with that of the analytical solution of Yan et al. (2021a) surface in Case 2 (Fig. 2d). The model parameters are listed in Table 3. In this study, a low stress or applied waste load of 50 kPa is considered following Woodman et al. (2014) and Das and Bharat (2021) as well as to be consistent with the 1D simulations presented in Zhang et al. (2018) for allowing comparisons between 1 and 3D model predictions. Following Kaczmarek and   Hueckel (1998), in this study, it is assumed that the application of mechanical load is instantaneous meaning soil has no time to deform and the applied load is accommodated by the pore fluid. This is reflected at the initial pore fluid pressure considered in the simulations. The 3D model domain is discretized using tetrahedra mesh elements. For the Case 1 simulation, Fig. 2c, uniform mesh elements of 0.14 m is used. Whereas, in Fig. 2e, near the circular leakage, mesh element varies between 0.0008 and 0.08 m and in the rest of the domain between 0.0006 and 0.14 m. The simulation runtime is 50 years, and the discretized time-step between 0 and 1 year is 3.65 days and, between 1 and 50 years, it is 36.5 days.
The clay liner properties, of the presented simulations, are obtained from Kaczmarek and Hueckel (1998) and Peters and Smith (2004). The model application example of Kaczmarek and Hueckel was extended by Peters and Smith to study NaCl leachate through a typical landfill liner. However, additional parameters were required to support the model simulations presented in this study. For example, the shear modulus data was obtained from Bengtsson and Sällfors (1983), who measured the value for soft, plastic Swedish clays representative of a typical landfill clay liner. The longitudinal and transverse dispersion coefficients, for a natural clay liner, were obtained from Zhang et al. (2015), and the empirical parameter, m, to calculate tortuosity (Eq.21), was obtained from the effective diffusion coefficient and porosity data.

Barrier system including single compacted clay liner (Case 1)
In this case, the CCL is directly exposed to leachate, and the solute concentration at any point on the upper boundary (z = 0 m) of the liner is same as that in the leachate. The initial and boundary conditions of the simulations are presented in Table 4.
Evolutions of excess pore fluid pressure in the Case 1 model domain are shown in Fig. 3. The cloud diagrams of u correspond to the results after 0.1, 1.0, 10, 20, and 50 years of simulation.
The cloud diagrams of clay liner settlement at different times are shown in Fig. 4(a)-(d). The settlement profiles along the length of the liner depth with time are presented in Fig. 4(e). Two stages, e.g., (I) the increasing period and (II) the decreasing period, can be clearly seen from Fig. 4. The reason why the settlement evolution presents two stages is because mechanical consolidation and chemico-osmosis consolidation are key factors controlling soil deformation which occur successively with simulation time. As mentioned previously that, at the initial stages, mechanical consolidation is the dominant factor affecting soil settlement, and with the dissipation of positive excess pore water pressure, soil deformation continues. However, with time, mechanical consolidation is completed, and chemico-osmosis consolidation caused by chemical loading dominates. Meanwhile, as the solute transport continues into the clay liner, solute concentration difference across the semi-permeable clay gradually decreases, which reduces the degree of chemico-osmosis consolidation, leading to the rebound of soil deformation predicted in stage (II).
The cloud diagrams of solute concentration distribution at 0.1, 1.0, 10, 20, and 50 years are shown in Fig. 5. The results predict that the solute transport continues downward with time and the penetration distance of the solute at 10 years is 0.62 m, and at 50 years, the solute penetrates the entire clay liner. The penetration distances of solute are one of the most important parameters for estimating the liner performance, which is defined as the distance where the in situ solute concentration reaches 10% of the source concentration, i.e., c/c 0 = 0.1 (Fig. 5e). Solute concentration at the half-length of the liner depth reaches to 45% of that of the source after 50-year simulation period.

Barrier system including geomembrane-compacted clay composite liner (Case 2)
In this case, the clay liner is exposed to leachate via a damaged GMB/CCL, which is illustrated in Fig. 1b. The initial and boundary conditions of the simulations are presented in Table 5. The cloud diagrams of excess pore fluid pressure evolution under local leakage condition are presented in Fig. 6. The results show that during the early simulation period (0.1 year), excess pore fluid pressure is positive everywhere except around the leakage area. The positive pore fluid pressure is associated with the mechanical consolidation. However, the excess positive pressure disappears rapidly, and negative pressure develops due to chemical-osmosis. Figure 7 shows the settlement along the depth of the clay liner at different simulation time under local leakage condition. The settlement of clay liner directly below the leakage area, which is controlled by mechanical consolidation and chemical-osmotic consolidation, is greater than that of surrounding soil. Furthermore, with the completion of consolidation, the soil settlement increases gradually, reaches the maximum in about 1 year, and seems to change marginally afterwards. Figure 8 shows the spread of the solute in the model domain for various Case 2 simulation periods. After 50 years, horizontal spread of the solute, at the upper surface of the clay liner, reaches to 0.29 m, and vertical penetration to around 0.382 m. The results show significant improvement of liner performance due to consideration of the geomembrane. The breakthrough concentration (c/c 0 = 0.1) reaches only up to 0.4 m of the depth of the liner for the duration of the simulation period that is significantly less than the predicted concentration in Fig. 5(e). Solute breakthrough, across the depth of the clay liner, did not occur during the 50-year simulation period.

Influence of geomembrane on hydraulic-mechanical-chemical response of clay liner
In order to investigate the influence of geomembrane on hydraulic-mechanical-chemical response of clay liner, evolution of excess pore water pressure, soil deformation, and solute transport at the center line of clay liner for the Case 1: compacted clay liner system and Case 2: geomembrane-compacted clay composite liner system, have been compared here.
The comparison of excess pore fluid pressure for the two cases are presented in Fig. 9. It shows that the negative pressure corresponding to Case 2 is less than that of Case 1. This is because the chemical loading at the upper boundary of clay liner in Case 2 is much smaller than that in Case 1. Furthermore, it is worth noting that the dissipation rate of excess pore pressure in Case 2 is slower, and the peak negative pressure occurs at depth around 0.1 m during the simulation period while, in Case 1, it moves progressively from 0.1 to 0.45 m along the depth of the liner. This variation reflects the importance of multi-dimensional hydrodynamic dispersion in solute transport, which occurs in both vertical and horizontal directions in Case 2 and weakens the spread or reduction   The results demonstrate the usage of geomembrane on reducing the soil settlement, rebound, and undulation. Figure 11 shows the graphs of solute concentration with depth in Cases 1 and 2. Expectedly, the solute transport is faster in Case 1 than Case 2. The reason is that in Case 1, contaminant loading is uniform at the xy-plane which results into higher contaminant flux than Case 2 where the load is applied as a point load. Also, in Case 2, both longitudinal and transverse spreads result into slower movement of the contaminant which experiences only one directional (longitudinal) spread, i.e., along the z-axis in Case 1 due to uniform loading condition.

Significance of model dimensionality on predicted results
To evaluate the influence of dimensionality on the coupled HMC response of the clay liner, excess pore water pressure, solute transport, and/or soil deformation results under onedimensional and three-dimensional conditions are compared.

Comparison of 1D vs 3D coupled HMC model for compacted clay liner system
The comparison of excess pore fluid pressures, soil deformation, and solute spread results for the 1D and

Comparison of the simulation results of 1D and 3D HMC coupled models in geomembrane-compacted clay composite liner system
The comparison of excess pore fluid pressure developments for 1D and 3D HMC coupled models in the geomembranecompacted clay composite liner system is presented in Fig. 13a. It shows that, at the initial stages (0.1 year), the negative pressure development in the 3D coupled model is expectedly greater than that of the 1D model due to higher dissipation of pore fluid through the more available fluid boundaries (Fig. 13a). Although with simulation time, the variations between the two pore fluid pressure profiles (1D vs 3D) are reduced, the deviations, especially at the vicinity of the contaminant source, remain largely noticeable for the damaged GMB/CCL system. In comparison, for the CCL system, the deviations between the two profiles became negligible (Fig. 12a).
The predicted solute concentration values corresponding to the 1D HMC coupled model are greater than those of the 3D model (Fig. 13b). The reason is that contaminants migration in the transverse direction is not considered in 1D models which eventually accelerate its spread only in the longitudinal direction (in this case, the vertical z-direction). This suggest that the 1D coupled models significantly overestimate solute transport comparing to the 3D models.
Results of the time-dependent settlement behavior of the GMB/CCL liner are presented in Fig. 13(c). It shows that the peak settlement of the liner, in the 3D HMC coupled model, is 6 mm which is roughly 26% of the peak settlement (23 mm) observed during the 1D simulations. That means, analogous to the contaminant spread, liner settlement is also overestimated during 1D simulations for damaged GMB/CCL systems.

Conclusions
In this study, a fully coupled 3D hydro-mechanical-chemical (HMC) model is presented to investigate long-term solute transport, soil settlement, and excess pore water pressure evolution in landfill clay liners. The multi-physics finite element software COMSOL Multiphysics was used to develop the model and conduct model simulations. The application of the model focused on solute transport and liner deformation behaviors in a single compacted clay liner and a geomembrane-compacted composite clay liner subjected to leachate leakage conditions.
The results suggest that under field conditions, during osmosis, fluid can flow in all directions across semi-permeable membranes (i.e., clays) and the phenomenon is accurately represented by a 3D model. Lower dimensional models are limited to specific directional flow and therefore lacks true representation of field conditions. This is further reflected on solute spread results predicted in the liner through a damaged geomembrane. Despite the same source concentration, spread of solute at a certain period of time in 3D simulations seems to be less than that of the lower dimensional model simulation. Effect of model dimensionality on solute spread for CCL is found to be marginal, but for GMB/CCL, the effect is significantly large. After 50-year simulation period, solute concentration at the half-length of the GMB/CCL liner is predicted to be 40% of the source concentration during 1D simulation, which is only 6% during the 3D simulation. Solute spread accelerates vertically than horizontally, since overburden load and consequent mechanical loading-induced solute convection occur in the same direction despite the soil being homogeneous and isotropic. Moreover, the results revealed approximately 74% over-estimation of liner settlement in 1D simulation than that of the 3D for GMB/CCL system. This highlights the importance of a three-dimensional hydro-mechanical-chemical model to study solute transport under mechanical and chemico-osmotic loading conditions in semi-permeable clays. The application of the model revealed the utility of geomembrane in compacted clay liners on reducing the soil settlement, rebound, and/undulation as well as restricting the spread of the solute in the liner.

Data availability
The data presented in this manuscript are available on request to the first author and the second author.

Declarations
Ethics approval Not applicable.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.