Compressibility Behavior of Bentonites by Stern Theory based on Constant Surface Charge Conditions

Abstract

The compressibility behavior of clays is governed by the electrical double layer formed around the clay particles. The Gouy-Chapman diffuse double layer theory is often utilized to predict the compressibility behavior of clay minerals. The theory does not consider the effect of the size of the cations, however, and thus predicts unrealistically small void ratios for compacted bentonites under large mechanical pressures expected in high-level nuclear waste-repository applications. In this study, the Stern layer was introduced to incorporate the cation size effect in the prediction of the compressibility behavior of bentonites. The overall diffuse double-layer thickness at large pressures was much smaller than the initially assumed Stern layer thickness based on the exchangeable cation size for all the bentonites studied. A compressible Stern layer was, therefore, considered for the first time in the prediction of the compressibility behavior of bentonites. The compression behavior of the Stern layer under the applied loading is influenced by the ratio of the mid-plane to the Stern potential, which is dependent on the type and composition of the exchangeable cations on the clay surface. Stern layer compression was initiated when the potential ratio was in the range 0.65–0.75 for bentonites with various surface cation characteristics. The incorporation of cation size and a compressible Stern layer provided significant improvements over the existing models in predicting the compressibility behavior of bentonites over a wide pressure range. The compressibility data predicted by the proposed model showed very good agreement with the data measured for five bentonites from the literature in the pressure range 0.1–42 MPa.

Introduction

Bentonites consist predominantly of the expansive smectite group of minerals and exhibit attractive features such as significant ion-adsorption capacity, large swelling capacity, and very poor hydraulic conductivity (Benson et al., 1994; Chen et al., 2016; Glatstein & Francisca, 2015; Kaufhold et al., 2015). Bentonites are used widely in various geotechnical and geoenvironmental engineering-field applications. Compacted bentonites have been considered as buffer and backfill materials for underground, high-level nuclear-waste repository systems in many countries (Bharat et al., 2009; Bharat et al., 2013; Butcher & Müller-Vonmoos, 1989; ENRESA, 2000; Ishikawa et al., 1990; Japan Nuclear Cycle Development Institute, 1999; Pusch, 2015; Tripathy et al., 2004; Zheng et al., 2017). These facilities are being planned at a depth of ~500 to 1000 m below ground level in various countries (Atomic Energy of Canada Limited (AECL), 2002; Enviros, 2003). The geostatic stress at such a depth is expected to be in the range 9–40 MPa (Tripathy & Schanz, 2007). The magnitude of stress in landfill liners and tailing impoundments is usually expected to be in the range 0.36–6 MPa (Peirce et al., 1986; Timmons et al., 2012). Several studies considered the stress range 3–42 MPa for studying the compressibility behavior of bentonites in these applications (Baille et al., 2010; Bharat et al., 2013, 2020; Marcial et al., 2002; Pusch et al., 2011; Tripathy & Schanz, 2007; Ye et al., 2014). Laboratory estimation of compressibility behaviors at high pressures, however, is very time-consuming and expensive as it requires specialized heavy equipment and loading mechanisms (Ng et al., 2006; Tripathy & Schanz, 2007).

Empirical models have been proposed in the past to predict the compressibility behavior of natural soils (Bharat & Sridharan, 2015a; Burland, 1990; Nagaraj & Murthy, 1986). The applicability of these models, however, is limited to a certain range of soil plasticity and lower ranges of applied pressures. The Gouy-Chapman model for interacting parallel clay-water-ion systems has been used to predict the compressibility behavior of clays (Bharat & Sridharan, 2015b; Bolt, 1956; Sridharan & Jayadeva, 1982; Tripathy and Schanz, 2007). Discrepancies have been observed between theoretical predictions and the measured compressibility data primarily due to the assumption of the parallel arrangement of the clay platelets and the treatment of the cations as point charges in the theory (Bolt, 1956; Warkentin et al., 1957). The effect of the size of the cations was incorporated (Stern, 1924) by introducing a thin and compact layer of cations next to the clay-platelet surfaces using the original Gouy-Chapman model (DDL). The Stern model has been utilized to study the electrical potential distribution of non-interacting clay platelet systems (Shang et al., 1994; Sridharan & Satyamurthy 1996; van Olphen, 1977; Verwey et al., 1948) and the compressibility behavior of bentonites based on the constant surface potential (CSP) condition (Tripathy et al., 2014). The Stern theory at constant surface charge condition (CSC) is, however, favored for the clays as their basal surfaces possess constant/permanent charges. Prediction of clay compressibility behavior using the interacting Stern theory under the CSC condition is not yet possible because the mathematical formulation of electrostatic potential distribution is yet to be established.

The purpose of the present study was to present an improved predictive model for the compressibility behavior of bentonites by considering the effect of the sizes of cations. The Stern DDL theory under the CSC condition was utilized based on the postulation that clay platelets were in a parallel arrangement under large applied pressures. The Stern layer thickness, however, is still not well defined for clay minerals. Although most of the available studies consider the Stern layer to be incompressible, results from the present study showed that, based on the measured void ratios, the DDL thickness at large applied pressures is much smaller than the Stern layer thickness. The proposed model thus incorporates the compressibility of the Stern layer thickness which depends on the ratio of mid-plane to Stern potential as identified in the present study. The measured compressibility data of five different bentonites representing a wide range of surface area and surface cations in the pressure range 0.1–42 MPa were considered from the literature to validate the proposed model. The Gouy-Chapman model and the Stern model for CSP conditions were also considered for the comparative assessment.

Diffuse Double Layer Theory

The interaction of clays with water and electrolytes is important in understanding the engineering aspects of clays, such as the volumechange behavior, chemical sorption, and flow-related problems. Clay–water interaction involves physicochemical forces because of the electrochemical activity of the clay surface. The Van der Waals (VdW) attractions, capillary interactions, Coulombic attraction and repulsion, and long-range diffuse double-layer repulsive forces are the important surface forces that are known to exist in clays (Bishop, 1959; Bolt, 1956; Israelachvili, 2011; Lambe, 1960; Lambe & Whitman, 2008; Lu & Likos, 2006; Mitchell & Soga, 2005; Skempton, 1960; Sridharan & Rao, 1973; van Olphen, 1977, Verwey et al., 1948). While the Coulombic forces are negligible in expansive clays dominated by montmorillonite minerals, the capillary forces are absent at full saturation (Lu & Likos, 2006; Schubert, 1975). The VdW forces are significant at smaller interparticle separation distances (Israelachvili, 2011); their influence on the compressibility behavior is not well understood thus far, however. The VdW forces may be considered to be of passive type during compression as they are compressive in nature and will remain inactive during the compression loading. The magnitude of the VdW forces increases under compressive stress, however, due to the reduced separation distance or enhanced particle–particle interaction. The increased VdW forces are activated during stress removal and control to a significant extent the rebound or swelling response of the soil. The long-ranged diffuse double layer (DDL) repulsive forces, on the other hand, are predominant in saturated montmorillonite clays and primarily control the compressibility behavior (Bolt, 1956; Mitchell, 1960; Olson & Mesri, 1970; Sridharan & Rao, 1970, 1972). The Vdw forces are, thus, often ignored and the compressibility behavior of clays is obtained based on the equilibrium between the applied mechanical stress and the repulsive forces. The Gouy-Chapman DDL theory is commonly used to understand the clay–water–electrolyte interaction, which relates the repulsive forces to the electrostatic potential distribution in the clay–water system (Bharat & Sridharan, 2015a, 2015b; Bharat et al., 2013; Bolt, 1956; Honig & Mul, 1971; Komine & Ogata, 2003). A brief description of the diffuse double layer theory is presented below followed by the theoretical formulation of the compressibility behavior.

Net negative charges are available on the basal surfaces of clay platelets (montmorillonite minerals) due to the isomorphous substitution of Al³⁺ by Mg²⁺ in the crystal structure of the octahedral alumina sheet (Grim, 1968; Mitchell & Soga, 2005). Exchangeable cations are naturally present on the clay surface to compensate for the negative charges. In the presence of a water/electrolyte medium, the cations on the clay surface experience an additional diffusive type of force that tries to drive the cations away from the charged clay surface. The diffusive forces are developed due to the existing concentration gradient of the cationic species between the clay surface and the bulk electrolyte solution. An electric diffuse double layer (cation cloud) is formed around the clay platelets as a consequence of the competition between the strong electrostatic attraction between the cations and negatively charged clay surface, and the diffusive forces (Sparks, 1999; Van Olphen, 1977; Verwey et al., 1948). The DDL around the clay platelet primarily controls the interaction among the clay platelets and often leads to the parallel plate orientation in the saturated clays. For such a parallel plate configuration, the DDL around the clay platelets defines the separation distance among the clay platelets which is related to the macroscopic void ratio. Thus, the compressibility behavior of clays is predicted by the theoretical estimation of the thickness of the DDL around clay platelets (Bolt, 1956; Sridharan & Choudhury, 2002; Sridharan & Jayadeva, 1982; Tripathy and Schanz, 2007).

An illustration of the interacting parallel-plate clay-water-electrolyte system under applied mechanical pressure by the Gouy-Chapman DDL model (Chapman, 1913; Gouy, 1910) is presented in Fig. 1. As the two clay platelets approach each other under the action of applied mechanical stress, repulsive pressure (P _R) develops between the clay platelets due to the interaction of the similarly charged DDL around the clay platelets. The separation distance between two neighboring clay platelets continues to decrease under the applied mechanical stress until it reaches an equilibrium state. Under this condition, the repulsive pressure is equal to the applied mechanical pressure (P) on the system (Bharat & Das, 2017; Bharat & Sridharan, 2015a). The electrostatic potential distribution within the interacting system is represented by the curve y(x) with a minimum potential ( $u$ ) at the mid-plane as shown in Fig. 1. The midplane potential is related to the repulsive pressure between the two platelets or the applied mechanical pressure at equilibrium as per the equation given by Langmuir (1938) presented below (Mitchell & Soga, 2005; van Olphen, 1977):

(1) $P_R=P=2nkT\left(\cosh ,\left(u\right),-,1\right)$

where n is the concentration of cations in the electrolyte solution in ions/m³, k is the Boltzmann constant ( $1.38\times 10^{-23}J/K$ ), and T is the temperature in K. The mid-plane potential is uniquely related to the separation distance between the clay platelets and represents the degree of interaction in the system. The relationship between the mid-plane potential and the separation distance is given by the Poisson-Boltzmann equation as (van Olphen, 1977):

(2) $t_{\mathrm{DDL}}=-\frac{1}{\kappa }{\int }_{y_0}^u(2\cosh \left(y\right)-2\cosh \left(u\right){)}^{-1/2}dy$

where t _DDL is the DDL thickness, which is equal to half of the separation distance (d) in the Gouy-Chapman DDL model, $\frac{1}{\kappa }$ ( $\kappa =\sqrt{\frac{8\pi q^2{v}^{2}n}{\epsilon D_{0}kT}}$ ) is the Debye length (Bharat & Sridharan, 2015b; Bharat & Sridharan, 2015b; van Olphen, 1977; Verwey et al., 1948), q is the charge on the electron ( $=1.6\times {10}^{-19}C$ ), v is the valence, ε is the dielectric constant, $D_0$ is the dielectric permittivity of vacuum ( $8.854\times 10^{-12}{C}^2/N-m^2$ ), and $y_0$ is the normalized electrostatic potential at the clay surface. The surface potential estimated for a given soil surface and electrolyte properties are given by:

(3) ${\left(\frac{\mathrm{dy}}{d\xi }\right)}_{x=0}=\sqrt{2\cosh \left(y_0\right)-2\cosh \left(u\right)}=\sigma \sqrt{\frac{1}{2\epsilon D_{0}nkT}}=0.96352\frac{C_e}{S_a}\sqrt{\frac{1}{2\epsilon D_{0}nkT}}$

where ${\left(\frac{\mathrm{dy}}{d\xi }\right)}_{x=0}$ is the slope of the potential distribution curve near the clay surface, σ is the total surface charge density on the clay surface, C _e is the cation exchange capacity (meq/100 g), and S _a is the specific surface area (m²/g). The macroscopic void ratio is related to the inter-platelet separation distance as per the following equation (Bharat & Das, 2017; Bolt, 1956):

(4) $e=G{\rho }_w{S}_a\frac{d}{2}$

where  $d/2$ is half of the separation distance between two parallel clay platelets and G is the specific gravity. The void ratio at any given pressure, thus, can be estimated using Eqs. 1–4 assuming a parallel plate orientation of the clay platelets. The integral involved in Eq. 2 is elliptic, which is to be solved numerically to establish the relationship between mid-plane potential and separation distance (Bharat et al., 2013).

Fig. 1 Illustration depicting the interaction of DDL of two approaching clay layers under applied mechanical pressure

Stern DDL Theory

The finite size of the cations at the particle surface limits the closest approachable distance to the charged clay surface (Stern, 1924). This results in a relatively compact and immobile layer of counter-ions close to the surface, which is followed by a diffused layer of the counter-ions. Thus, the electric double layer in a clay-water system is characterized by the Stern layer and the outer diffused layer consisting of the Gouy layer. The center of the spherical cations in the Stern layer is at a distance approximately equal to their hydrated radius away from the clay surface (Guven & Pollastro, 1992), which is defined as the outer surface of the Stern layer (Shang et al., 1994). The charge within the Stern layer is constant and the electrostatic potential varies linearly from a maximum value ( $y_0$ ) at the clay surface to $y_{\delta }$ at a distance equal to the Stern thickness (δ) at the Stern-Gouy interface where the potential is termed the Stern potential. The dielectric constant of water within the Stern layer is reduced significantly to 3–6 as the water molecules are bound tightly to the clay surface (Hunter, 1981;Shang et al., 1994; Sposito, 1984; Sridharan, 1968; Sridharan & Satyamurty, 1996; van Olphen, 1977; Verwey & Overbeek, 1955). A graphical illustration of the Stern model for the interacting clay–water system is presented in Fig. 2. The inter-platelet separation distance, d, in the Stern model is the summation of the Stern layer thickness, δ, and the thickness of the Gouy diffuse layer, t _DDL.

Fig. 2 Electric potential distribution in the Stern DDL model

(5) $d=2(\delta +t_{\mathrm{DDL}})$

The electrostatic potential distribution within the Stern layer is dependent on the surface charge density and the dielectric properties of the pore fluid as given in Eq. 6

(6) $\sigma =\frac{{\epsilon }^{{}^{\prime }}kT}{4\pi \delta vq}\left(y_0,-,y_{\delta }\right)$

Within the Gouy layer, the electrostatic potential distribution varies non-linearly between the Stern potential (y _δ) at x = δ and the mid-plane potential (y _d) at x = d/2, which is represented by the Poisson-Boltzmann equation:

(7) $\kappa t_{\mathrm{DDL}}=-{\int }_{t_{\mathrm{DDL}}-\delta }^{t_{\mathrm{DDL}}}d\xi =-{\int }_{y_{\delta }}^u(2\cosh \left(y\right)-2\cosh \left(u\right){)}^{-1/2}dy$

Evaluation of the Stern potential ((y _δ) is a prerequisite for the estimation of Gouy layer thickness (t _DDL). The relationship between the charge density and potential distribution is utilized to determine the Stern potential for a given clay-water-electrolyte system. The Stern layer charge (σ₁) and Gouy layer charge (σ₂) together balance the total negative surface charge (σ) on the clay platelets:

(8) $\sigma =-({\sigma }_1+{\sigma }_2)$

The clay platelets are treated as constant-charged plates (Grim, 1968), for which the total surface charge density can be expressed as:

$\sigma =0.96352\frac{C_e}{S_a}\text{C/}{\text{m}}^2$

where C _e is the cation exchange capacity of the clay mineral (meq/100 g) and S _a is the specific surface area (m²/g). The charge density in the Stern layer can be obtained as (Verwey et al., 1948):

(9) ${\sigma }_1=\frac{N_{1}vq}{1+\left(\frac{N_A}{\mathrm{Mn}}\right)\exp \left(-,\left(y_{\delta },+,\frac{\psi vq}{kT}\right)\right)}$

where N ₁ is the number of adsorption sites per 1 cm² area of the clay surface, v is the valence of ions, N _A = Avogadro’s number, M = molecular weight of the solvent (water), y _δ is the Stern potential at the plane separating the Stern and Gouy layers, ψ = specific adsorption potential on the counter-ions at the surface. The charge in the Gouy layer can be derived as (Verwey et al., 1948; van Olphen, 1977):

(10) ${\sigma }_2=\sqrt{2nkT\epsilon }\sqrt{2\cosh y_{\delta }-2\cosh u}$

Combining Eqs. 6 and 8–10, the Stern potential can be expressed as a function of the mid-plane potential and the pore fluid parameters, which in turn can be used along with Eqs. 1 and 6 to estimate the void ratio at a given pressure. Equation 9 is only valid for non-interacting systems (van Olphen, 1977) and requires modification if applied to the interacting clay-water system. The number of available sites in the bulk solution, N ₁, depends on the volume of the diffuse layer which is subjected to change with the change in the degree of interaction under the applied pressure. Van Olphen (1977) presented the following equation to replace Eq. 9 for an interacting system,

(11) $\frac{{\sigma }_1}{{\sigma }_2}=\left(\frac{\delta }{\left(\frac{d}{2},-,\delta \right)}\right)\exp \left(y_{\delta },+,\left(\frac{\psi }{kT}\right)\right)$

The inter-particle distance can be estimated from Eq. 12,

(12) $\kappa t_{\mathrm{DDL}}=-{\int }_{\delta }^{d/2}d\xi =-{\int }_{y_{\delta }}^u(2\cosh \left(y\right)-2\cosh \left(u\right){)}^{-1/2}dy$

Equation 11 is derived based on the assumption that statistical charge distribution between the Stern and Gouy layer is proportional to their respective volumes or their respective thickness. The equation leads to an erroneous estimation of electrostatic potential distribution, however, as it does not utilize the correct volume for the Gouy diffused layer.

The interacting Stern model is utilized for predicting the compressibility behavior of clays by assuming the clay surface potential to be constant (Tripathy et al., 2014). Because montmorillonites have constant surface charge (Grim, 1968), the assumption of constant surface potential may not be tenable. Therefore, estimation of Stern potential and DDL thickness requires modification of the Stern theory, which has been dealt with in the present study.

Stern-Gouy Model for Constant Surface Charge

According to van Olphen (1977), the charges in the Stern and Gouy layers can be proportional to their respective areas under the electrostatic potential distribution curve (Eq. 13).

(13) $\frac{{\sigma }_1}{{\sigma }_2}=\frac{\text{Area of the Stern layer}\left(A_s\right)}{\text{Area of Gouy layer}\left(A_{\text{Gouy}}\right)}$

Estimation of a Stern layer charge (σ₁) in an interacting system for CSC conditions is not possible due to the difficulties involved in estimating the two parameters, such as the number of available adsorption sites and the specific adsorption potential of the counter-ions at the clay surface (see Eq. 9). The influence of platelet interaction on the number of available adsorption sites is not understood yet. The specific adsorption potential, ψ, for a given clay-water-electrolyte system is difficult to estimate under varying DDL interaction. Here, the estimation of σ₁ was eliminated and the following equation was developed to determine the Stern potential at a given pressure by knowing the clay mineral surface and pore-fluid properties. (A detailed derivation is presented in the Appendix.)

(14) $\sigma A_{\mathrm{Gouy}}=\sqrt{2nkT\epsilon }\sqrt{2\cosh y_{\delta }-2\cosh u}\left(\left(y_{\delta },+,\left(\frac{vq\sigma 2\pi \delta }{{\epsilon }^{{}^{\prime }}kT}\right)\right),\delta ,+,A_{\mathrm{Gouy}}\right)$

The area under the hyperbolic potential distribution curve for the Gouy layer (A _Gouy) is a function of Stern and mid-plane potentials which can be calculated by following the method of slices (see Appendix). Therefore, Eq. (14) provides an implicit solution for the Stern potential. For a known value of mid-plane potential, Stern potential is obtained through optimization. The objective function to determine the Stern potential, based on Eq. (15), is given as:

(15) $f\left(y_{\delta }\right)=\sigma A_{\mathrm{Gouy}}-\left(\sqrt{2nkT\epsilon },\sqrt{2\cosh y_{\delta }-2\cosh u},\left(\left(y_{\delta },+,\left(\frac{vq\sigma 2\pi \delta }{{\epsilon }^{{}^{\prime }}kT}\right)\right),\delta ,+,A_{\mathrm{Gouy}}\right)\right)$

The ‘fminbnd’ function, which is based on the golden section search and parabolic interpolation method, was used to obtain the optimized value of the Stern potential from the objective function (Eq. 15) in Matlab. The mid-plane potential, u, was used as the lower boundary in the optimization, to which was added a small value (~10^–9) to avoid the singularity (Bharat et al., 2013). The upper boundary was fixed at 30 for the studied pressures and pore-fluid concentrations. Further, the nature of the objective function was studied at three different applied pressures and four different pore-fluid concentrations. The important parameters used in the objective function evaluation are presented in Table 1. For all the cases, the true minima are preceded by a minimum at the lower boundary (i.e. u). The true minima approach the first minimum (u) with the increase in the applied pressure (Fig. 3a) and pore-fluid concentration (Fig. 3b). Overall, the local minimum was observed in the range 0–10 for all the cases considered.

Table 1 Parameters used to establish the electrostatic potential distribution in the Gouy-Chapman model and the proposed Stern model (Figs. 3, 4)

Fig. 3 a Nature of the objective function at 0.0001 M pore-fluid concentration for three different pressures. b Nature of the objective function at 10 kPa pressure for four different pore fluid concentrations

Electrostatic Potential Distribution

The potential distribution for the entire DDL in the Gouy-Chapman model follows the Poisson distribution from a maximum value at the surface to a minimum at the mid-plane. In the case of the Stern model, the potential starts with a maximum value at the surface and decreases linearly to the Stern potential at the Stern–Gouy interface. Beyond this, the potential follows the Poisson distribution within the Gouy-layer to a minimum value at the mid-plane. The mid-plane potential, u, at a given pressure was determined using Eq. 1. The surface potential was obtained using Eq. 3 for the Gouy-Chapman model and Eq. 6 for the Stern model after determining the Stern potential. The Stern potential was estimated through optimization using Eq. 15 for known u.

The potential distributions for both models from the surface to mid-plane distance under two different applied pressures at equilibrium are shown in Fig. 4a. The parameters considered in the computation are presented in Table 1. At the lesser applied pressure (0.01 MPa), the influence of the size of the cations was only visible near the clay surface up to a distance of ~20 Å. At the greater applied pressure, the effect of cation size on the potential distribution was more pronounced due to the increased DDL interaction as the separation distance reduced significantly.

Fig. 4 a Computed electrostatic potential distribution in the clay-water system (n = 0.0001 M) based on the GC and Stern models under two different applied mechanical pressures. b Variation in the Stern and mid-plane potentials in the clay-water system (n = 0.0001 M) with changes in the degree of interlayer interaction (separation distance) under loading

The variations in Stern potential, mid-plane potential, and the DDL thickness (on the third axis) of an interacting clay–water system with the applied pressure by the proposed Stern model are presented in Fig. 4a. The values of all parameters considered in the simulation are presented in Table 1. The applied pressure varied in the range of 0.01 to 40 MPa. The DDL thickness decreased exponentially with the applied pressure and attained a minimum thickness of 7.9 Å (hydrated radius of Na⁺ cation), equivalent to the Stern thickness at ~5 MPa. This indicated the full compression of the diffused Gouy layer leading to compact layers of cations around the clay platelets. The DDL thickness thus remained constant with a further increase in the pressure. The mid-plane potential increased linearly with the increase in the applied pressure on the semi-log scale due to the increased DDL interaction. The Stern potential, on the other hand, showed a slower increase in magnitude as compared to the mid-plane potential. The two potential curves eventually converged above ~7–8 MPa pressure as the mid-plane coincided with the Stern boundary after the elimination of the Gouy-layer. The ratio between the Stern potential and the mid-plane potential, thus, is a useful parameter to understand the compressibility behavior of the DDL under the applied pressure.

Stern Layer Thickness at Large Pressure

Choosing an appropriate thickness of the Stern layer is crucial for predicting the pressure-void ratio relationship, especially in the higher-pressure range, where the Gouy diffused layer is compressed significantly. The minimum possible void ratio (i.e. the minimum separation distance between two interacting clay platelets) in clays is controlled by the thickness of the Stern layer. A well-defined value for the Stern layer thickness, however, is not available for the clay minerals (Verwey et al., 1948; van Olphen, 1977; Shang et al., 1994; Sridharan & Satyamurhty, 1996). The type of exchangeable cations, charge distribution, and size and shape of the hexagonal cavity on the surface of the montmorillonite influence the adsorption of cations (Sposito, 2008), which can significantly influence the Stern layer thickness. Most of the available studies consider the Stern layer to be incompressible and equivalent to the radius of the hydrated cations, but the interaction between the particles considered in such studies is weak or negligible. The behavior of the Stern layer under significant applied pressure has yet to be studied, however.

Experimentally determined compressibility data of seven different bentonites from the literature were considered in the present study to understand the minimum achievable separation distance between the clay platelets (DDL thickness) under the applied mechanical pressures. The relevant properties of the bentonites are presented in Table 2. The DDL thickness was derived from the experimental void ratio for these bentonites from the literature using Eq. 4 by considering the parallel plate assumption. The void ratio of the bentonites considered here was in the range ~0.4–1 under the studied pressure range. The parallel arrangement of the clay platelets has been reported for heavily consolidated clays at such small void ratios (Delage & Lefebvre, 1984). Strong DDL repulsion brings the clay platelets toward a parallel arrangement as the clay mineral is heavily compressed at large pressure. The DDL thickness was plotted against the applied pressure (Fig. 5a) and revealed that the thickness of the DDL was compressed to the smallest value of 2.3 Å in the pressure range of 10—40 MPa for different bentonites. The minimum possible separation distance between the two clay platelets surrounded by a rigid Stern layer is shown in Fig. 5b(i). When the DDL thickness or half of the separation distance is decreased beyond the value equivalent to the diameter of the exchangeable cation, the Stern layer thickness consequently becomes compressed. The Stern layer compression is facilitated by the penetration of the surface cations into the hexagonal cavities of the clay surface (Fig. 5bii) once the diffuse layer is eliminated from the system under a large applied pressure. Generally, the diameter of the hexagonal cavity is ~2.6 Å, which is ~1/3 of the hydrated size of Na⁺ (Sposito, 2008) and similar in size to that of a water molecule. The Stern layer, thus, becomes compressed under a large applied pressure to facilitate further volumetric compression of the clay once the Gouy-layer is compressed significantly.

Table 2 Relevant bentonite properties used in the theoretical prediction of compressibility behavior

^aTripathy et al. (2014), ^bMarcial et al. (2002), ^cDi Maio (2004), ^dLow (1980)

Fig. 5 a Variations in the DDL thickness with applied mechanical pressure derived from the measured compressibility data (using Eq. 4 with the parallel plate assumption) of bentonites having varying plasticity. b Illustration showing cations penetrating the hexagonal cavities on the clay surface, resulting in compression of the Stern layer under very high applied mechanical pressure

The total DDL thickness at large pressure was, therefore, corrected through the incorporation of the Stern layer compression. The compressibility of the Stern layer was, however, incorporated only into the void ratio computation in Eq. 4. The effect was not considered in the computation of the Stern potential and the thickness of the Gouy-layer because the theoretical formulation for such a complex interaction is not available.

The ratio of mid-plane to Stern potential (u/y _δ) during the compression of the DDL thickness under the applied pressure was studied for three different bentonites, namely Na-Kunigel, Ponza, and Na-Ca-MX80 (Fig. 6). The three bentonites represented a wide range of surface cation characteristics and surface charge densities (σ). The ratio between the two potentials indicated the degree of interaction between the two interacting clay platelets (Jiang et al., 2001). The relevant properties of the respective bentonites and other parameters related to the pore-fluid and Stern layer used in the estimation are presented in Tables 2 and 3, respectively. A cation concentration of 0.0001 N was used to represent water as a pore fluid (Das & Tadikonda, 2021). The Stern layer thickness was taken as the hydrated radius of Na⁺ for Na-dominated bentonite. For the divalent-dominated and mixed-valence bentonites, the larger cationic size (i.e., Ca²⁺) was considered as the Stern layer thickness.

Fig. 6 a Variation in the double-layer thickness with pressure at various degrees of interlayer interaction for Ponza bentonite. b Variation in the double-layer thickness with pressure at different degrees of interlayer interactions for Na-Kunigel bentonite. c Variation in the double-layer thickness with pressure at various degrees of interlayer interactions for Na-Ca-MX80 bentonite

Table 3 Parameters used in the prediction of compressibility behavior of the bentonites considered by the three models

^*weighted average valence (see Table 2), ^#1 for Na-dominated and 2 for divalent-dominated, ^$initial Stern layer thickness at zero pressure (see Table 2), ^&for Na⁺

The potential ratio increased with the applied pressure for all three bentonites as the DDL thickness was compressed, resulting in a greater degree of DDL interaction. The Stern layer compression began when the potential ratio was ~0.65 for the divalent-dominated Ponza bentonite as well as the mixed-valence Na-Kunigel bentonite (Fig. 6a). On the other hand, Stern layer compression was observed at ~0.75 for the Na-dominant Na-Ca-MX80 bentonite. The observed differences in the potential ratio at the beginning of the Stern layer compression for different bentonites was related to the variation in the surface cation characteristics and surface charge density. Stern layer compression started early for the bentonites containing greater surface-charge densities and divalent cations. Overall, the Stern layer compression began when the potential ratio was in the range 0.65–0.75 for various bentonites.

An S-curve relation between the Stern layer thickness (δ) and the ratio between the mid-plane and Stern potential was assumed to predict the void ratio, which follows the typical compressibility behavior of clays, as given by

(16) $\delta =r{\exp }^{-\left(a,\left(\frac{u}{y_{\delta }},-,R_y\right)\right)}\text{;}\frac{u}{y_{\delta }}-R_y>0$

where r is the hydrated radius of the cation (Å); R _y is the potential ratio at the onset of compression of the Stern layer and varies in the range 0.65–0.75, depending on the surface-charge characteristics of the clays. The parameter, a, defines the slope of the curve, which was determined based on the observed minimum achievable thickness of the Stern layer (~2.3 Å) when the potential ratio becomes unity, as given by:

(17) $a=(R_y-1)ln\left(\frac{2.3\stackrel{\circ }{A}}{r}\right)$

The above correction for Stern layer thickness was incorporated into the void ratio computation through Eq. 5 when the potential ratio reached a specified value under the applied pressure for a given clay mineral.

Validation of the Proposed Stern Model

The proposed compressible Stern model was validated using five different bentonites from the literature (Fig. 8). A flowchart for the computation of the pressure-void ratio relationship based on the proposed approach is presented in Fig. 7. The predicted compressibility data of these bentonites from the Gouy-Chapman model and the Stern model for CSP conditions (Tripathy et al., 2014) are also presented along with the proposed model to carry out a comparative analysis.

Fig. 7 Flow chart showing the computation of the void ratio at a given pressure using the proposed Stern model

The relevant properties of the bentonites used in the prediction of the compressibility data of the bentonites using the three models are presented in Table 2. The valence was taken as 1 for Na-dominated bentonites and as 2 for Ca or other divalent cation-dominated bentonites in the Stern interacting CSP model as considered by Tripathy et al. (2014), while weighted average valence (v _avg) was considered in the prediction by the Gouy-Chapman and the proposed model. The initial Stern layer thickness (δ₀) at zero pressure was considered to be equivalent to the hydrated radius of the Na⁺ cation for the Na-dominated Na-Ca-MX80 bentonite and Mexico montmorillonite in the proposed model. For the other bentonites, which have either mixed-valence or divalent-dominated surface cations, the hydrated radius of the Ca²⁺ cation was taken as the Stern layer thickness, being the largest among the available exchangeable cations. The Stern layer compression was applied when the potential ratio (u/y _δ) reached a specified value (R _y ) for the given bentonite. However, a fixed value of 5 Å was used as the Stern thickness for all the bentonites in the CSP model by Tripathy et al. (2014). The pore-fluid parameters and other relevant parameters used in the three models presented in Table 3.

The compressibility data predicted by the proposed model and the two existing models were compared with the measured data for five different bentonites (Fig. 8). The CSP model by Tripathy et al. (2014) was not in good agreement with the measured data in the pressure range for the bentonites studied here. The CSP model was in close agreement with the measured data for the Na-Kunigel bentonite briefly in the lower pressure range 0.1–0.5 MPa; the model deviated significantly at higher pressures (Fig. 8b), however. For the other bentonites studied, the CSP model was far away from the measured data as compared to the proposed model (Fig. 8a, 8c–8e). The observed discrepancies were attributed to the issue with the assumption of constant surface potential conditions in the model, as discussed earlier. Overall, the Gouy-Chapman model was relatively close to the measured data as compared to the CSP model; it severely underestimated the void ratios at large pressures, however. This was primarily due to not considering the effect of the size of the cation in the Gouy-Chapman theory. The proposed model, based on the Stern theory under CSC conditions, showed a better agreement with the measured data in the studied pressure range compared with the two existing models. The proposed model improved significantly the prediction at large pressures as the predicted void ratios were greater than the Gouy-Chapman model due to the incorporation of the size of the cations and were very close to the measured data. Treatment of the compressible Stern layer provided a realistic void ratio variation with pressure at very large pressures, unlike the earlier Stern model under CSP conditions. Overall, the proposed model showed a better prediction at pressures > 0.1 MPa for the studied bentonites, but overpredicted at the lower pressure range (0.01–0.1 MPa). The discrepancies observed at lower pressures were attributed to the dominant presence of the edge–face clay platelet orientation as the theory is based on the parallel arrangement of the clay platelets.

Fig. 8 Theoretically predicted and measured pressure-void ratio data for a Na-Ca MX80 bentonite (Marcial et al., 2002); b Na-Kunigel bentonite (Marcial et al., 2002); c MX-80 bentonite (Tripathy et al., 2014); d Ponza bentonite (Di Maio, 2011); e Mexico Montmorillonite (Low, 1980)

Conclusions

The effect of cation size was incorporated into the prediction of clay compressibility behavior using the Stern theory under the condition of constant surface charge (CSC) for the first time. A mathematical model was developed to establish the potential–distance relationship for the interacting Stern model at the CSC condition. The compressibility of the Stern layer was further incorporated into the theory for the first time to provide a more realistic prediction of the compressibility behavior of bentonites in the high-pressure range. Based on the detailed analysis of the compressibility data of various bentonites from the literature using the proposed Stern theory, the following conclusions were drawn:

Under the application of load, the diffused Gouy layer initially undergoes significant compression, while the Stern layer remains unaffected. The compression of the Stern layer starts in the pressure range 0.5–1 MPa for different bentonites once the thickness of the Gouy layer is significantly reduced. The thickness of the Stern layer reaches a minimum value equivalent to the water molecule size at a pressure of ~40 MPa. The hexagonal cavities on the surface of the clay platelets accommodate the cations at such high pressure to facilitate the Stern layer compression. The void ratio corresponding to the minimum Stern layer thickness at such high pressure is ~0.4.

The ratio of the mid-plane to Stern potential, which represents the degree of interlayer interaction in the clay-water system, influenced the compression behavior of the Stern layer under the load. The potential ratio of the clay-water electrolyte system at any given pressure is dependent on the type and composition of the exchangeable cations on the clay surface. The Stern layer compression starts when the potential ratio is in the range 0.65–0.75 for bentonites with different surface cation characteristics.