Stern DDL Theory

The finite size of the cations at the particle surface limits the closest approachable distance to the charged clay surface (Stern, Reference Stern1924). 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, Reference Guven and Pollastro1992), which is defined as the outer surface of the Stern layer (Shang et al., Reference Shang, Lo and Quigley1994). 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, Reference Hunter1981;Shang et al., Reference Shang, Lo and Quigley1994; Sposito, Reference Sposito1984; Sridharan, Reference Sridharan1968; Sridharan & Satyamurty, Reference Sridharan and Satyamurty1996; van Olphen, Reference van Olphen1977; Verwey & Overbeek, Reference Verwey and Overbeek1955). 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, Reference Grim1968), 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., Reference Verwey, Overbeek and Van Nes1948):

(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., Reference Verwey, Overbeek and Van Nes1948; van Olphen, Reference van Olphen1977):

(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, Reference van Olphen1977) 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 (Reference van Olphen1977) 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., Reference Tripathy, Bag and Thomas2014). Because montmorillonites have constant surface charge (Grim, Reference Grim1968), 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 (Reference van Olphen1977), 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., Reference Bharat, Sivapullaiah and Allam2013). 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., Reference Verwey, Overbeek and Van Nes1948; van Olphen, Reference van Olphen1977; Shang et al., Reference Shang, Lo and Quigley1994; Sridharan & Satyamurhty, Reference Sridharan and Satyamurty1996). 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, Reference Sposito2008), 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, Reference Delage and Lefebvre1984). 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, Reference Sposito2008) 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. (Reference Tripathy, Bag and Thomas2014), ^bMarcial et al. (2002), ^cDi Maio (Reference Di Maio, Santoli and Schiavone2004), ^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., Reference Jiang, Zhou, Zhu, He and Yu2001). 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, Reference Das and Tadikonda2021). 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.