Estimating vertical and lateral pressures in periodically structured montmorillonite clay particles.

Given a montmorillonitic clay soil at high porosity and saturated by monovalent counterions, we investigate the particle level responses of the clay to different external loadings. As analytical solutions are not possible for complex arrangements of particles, we employ computational micromechanical models (based on the solution of the Poisson-Nernst-Planck equations) using the finite element method, to estimate counterion and electrical potential distributions for particles at various angles and distances from one another. We then calculate the disjoining pressures using the Van't Hoff relation and Maxwell stress tensor. As the distance between the clay particles decreases and double-layers overlap, the concentration of counterions in the micropores among clay particles increases. This increase lowers the chemical potential of the pore fluid and creates a chemical potential gradient in the solvent that generates the socalled 'disjoining' or 'osmotic' pressure. Because of this disjoining pressure, particles do not need to contact one another in order to carry an 'effective stress'. This work may lead towards theoretical predictions of the macroscopic load deformation response of montmorillonitic soils based on micromechanical modelling of particles.


INTRODUCTION
The fundamental mechanisms by which soils change volume largely depend on particle size. On one hand, volume change in coarse-grained soils is governed first by deformation at particle contacts, then by rotational displacement and packing of the grains, or by bending and rearrangement of flaky particles (e.g., mica) (Terzaghi 1943). On the other hand, fine-grained soils, particularly clays, experience volume change that is also dependant on the composition of pore fluid (Bolt 1956, Mitchell andSoga 2005). Further, there is no need to direct contact among clay particles under low confining pressures, as required for coarse-grained soils (i.e. in clays, an effective stress does not necessarily imply particle contact, as it does in coarse grained soils).
These observations have been captured in a variety of mathematical models. In this work, we restrict our analyses to diluted montmorillonite clay particles (i.e., solid phase) -water (i.e., solvent) -and a monovalent salt (i.e., electrolytes) system, and investigate the particle-scale responses of this system to externally applied loadings for a variety of particle configurations.
Shrinkage and swelling of clays is a matter of study in a wide range of engineering disciplines, and though generally not in the consciousness of the general public, they have a tremendous impact on society (Murad and Moyne 2002): in civil engineering, with structural damage to buildings due to differential soil volume change (particularly important where montmorillonite is present); in petroleum engineering, with regards to borehole stability in expansive shales (usually the cap of oil and gas reservoirs) particularly when using water-based drilling muds, which change the pore fluid composition near the borehole walls; and in geoenvironmental engineering, when using clay liners as containment for contaminants (i.e., ions and other substances) or sealants of radioactive materials.

GOVERNING EQUATIONS
This section summarizes the governing equations for ion concentration and electric potential distributions around clay platelets. Based on these equations we then compute attractive and repulsive forces for particular particle configurations. Formulation of force balance allows then the determination of the swelling pressure in clays.

ION CONCENTRATION AND VOLTAGE DISTRIBUTIONS
As no analytic solution to the disjoining pressure is possible for realistic arrangements of particles, we employ relatively sophisticated micromechanical analyses based on the numerical solutions of the microscale Poisson-Nernst-Planck equations to determine electrical potential ψ and ion concentration c distributions in the porespace, assuming complete saturation of the soil (Cussler 1997, Smith et al. 2004: where the sub-index i refers to the i th ion (i.e., Na + and Cl − ), j i is the molar flux density, D 0i is the selfdiffusion coefficient (i.e., the diffusion coefficient of the i th ion in pure water), c i is the ion concentration, z i is the ionic valence, F = 96, 485.3 C • mol −1 is the Faraday constant, R = 8.3144 J • mol −1 • K −1 is the universal gas constant, and T is the absolute temperature. In the Poisson equation, ε 0 = 8.854 × 10 −12 F • m −1 is the permittivity of free space, and ε f is the relative permittivity or dielectric constant of the pore fluid (e.g., ε f ∼ 78 for water). It is assumed that the charge on particle surfaces is constant, i.e., An Acad Bras Cienc (2010) 82 (1) σ = const. The electrostatic boundary conditions applied to a charged particle surface can be expressed as (Stratton 1941, p. 36): where D f and D s are the electric displacement in the pore fluid phase ( f ) and solid particle phase (s) respectively, and α is either s or f .

SWELLING PRESSURE ESTIMATION
The swelling pressure in this system is estimated from the balance between osmotic and Maxwell pressures of electrostatic nature. For long-range interactions, the electrostatic forces are the ones dominating shrinkage and swelling behaviour (Murad and Moyne 2002). In a constant volume cell, this behaviour will imply in decreasing and increasing swelling pressures, respectively. For short-range interactions, attractive van der Walls forces need also to be considered in the force balance equations in addition to osmotic and Maxwell pressures (Anandarajah 1997, Santamarina et al. 2001). In the following, we restrict ourselves to long-range interaction between clay particles and, hence, neglect van der Walls forces. The electrically charged clay particles (with surface charge σ ) induce higher counter-ion concentrations in the micropores between particles, particularly as double layers overlap. This increase in ion concentration lowers the chemical potential of the pore fluid and creates a chemical potential gradient that causes pore fluid to want to flow into the micropores, forcing apart the clay particles through the osmotic pressure . For ideal solutions 1 , van't Hoff's Law relates the osmotic pressure to the electrolyte concentrations c i : where x 1 and x 2 are position vectors employed for formulation of force balance over an arbitrary volume element also denoted as 'tensor box' (see Fig. 1). The Maxwell stress tensor is used to estimate pressure contributions of electrostatic nature: where E = −∇ψ is the electric field vector and n is the unit vector normal to the boundaries of the tensor box, and I is the identity matrix. Using these definitions we can formulate force equilibrium at the tensor box which allows determining the swelling (or disjoining) pressure as the difference between equations (3) and (4) (Grodzinsky 2000, McCormack et al. 1995, Murad and Moyne 2002, Phillips 1999.

NUMERICAL RESULTS -DISCUSSIONS
The governing equations previously summarized are solved numerically using the finite element methods in a customized script written in the Comsol Package (COMSOL-AB 2007). The numerical models basically compromise pore fluid domains (sodium chloride solutions) that surround (rigid) 2D clay platelets. The boundary conditions are set such that symmetries are accounted for, and reference electrical potential and background concentration are prescribed in the far field to make the problem well-posed. Mesh refinements are required in the areas surrounding the particles and close to symmetry regions to achieve consistent convergence. A parametric study on optimal number of elements was performed to achieve such goals (not shown here).

VALIDATION: INFINITE PARALLEL PLATELETS
To validate the numerical finite element results obtained for ionic concentration and electrical potential distributions given by equations (1), we compare these numerical results with closed-form solutions: 1) a non-interacting double layer (i.e., single double layer) and 2) interacting double layers (i.e., electric potential between two close flat clay platelets). The non-linear Poisson-Boltzmann equation correctly describes ion and voltage distribution around charged clay particles for surface potentials ψ 0 less than 100mV and monovalent salt concentrations less than 100 mol • m −3 (0.1M) (Phillips 1999, Sposito 1984. This equation is valid within those surface potential and concentration limits, even under the following assumptions: 1) ions are considered point charges, 2) uniform clay surface charge distribution, 3) one dimensional condition (i.e., the particle surface is large relative to the double layer thickness), and 4) uniform effective pore fluid permittivity. In the non-interaction double layer case, the analytical solution for the electrical potential ψ as a function of the distance x normalized by the double layer thickness ϑ is (modified from Verwey et al. 1948): An Acad Bras Cienc (2010) 82 (1) where the double layer thickness ϑ = RT ε 0 ε f 2F 2 z 2 c ∞ , c ∞ is the background electrolyte concentration (which is the same for both ionic species in 1:1 electrolytes fulfilling electroneutrality), and ψ 0 is the surface potential.
In the case of two interacting double layers, that is, when two flat parallel clay platelets are close enough (i.e., less than approximately 5 times the single double layer thicknesses for typical montmorillonite clay surface charges), the closed-form solution to the mid-plane potential ψ d can be found from the following elliptic integral of the first kind Chen 1994, Verwey et al. 1948): where for simplicity, a dimensionless potential φ is defined as φ = z Fψ (RT ), thus φ 0 = z Fψ 0 (RT ) is the dimensionless surface potential, φ m is the mid-plane dimensionless potential, and d is the half-distance between clay platelets. Equation (7) needs to be solved numerically or with the help of tables (e.g. Jahnke and Emde 1945). Figure 2 shows some representative examples of the excellent agreement between our numerical results and the ones from equations (6) and (7). When the repulsive and attractive forces between infinitely long parallel clay platelets balance, equilibrium is reached leading to an equilibrium separation distance 2d. Conversely, given a particular equilibrium separation distance 2d, it is possible to find the corresponding swelling pressure T sw of external (e.g. external applied load) or internal origin (e.g. van der Waals -although these are less relevant at the long-range), in a way that it balances the repulsive forces, both electrical T e and osmotic (Grodzinsky 2000). Consider the one-dimensional case (x-direction) of two infinitely long parallel platelets with interacting double layers, and place a stress tensor box surrounding the left platelet (recall Fig. 1). If the left edge of the box is located at x 1 = x → −∞, then T e x x (x 1 ) = 0 as E x x (x 1 ) = 0; which is also the case if the right edge of the box is situated at x 2 = d due to symmetry. However, due to different concentrations at x = ±∞ and x = d, the osmotic pressure at (d) also differs from the osmotic pressure in the bulk solution (±∞). Therefore, the surface stress balance in the x-direction can be written as: The osmotic pressure terms in equation (8) for binary monovalent electrolytes can be shown to be a function of the background electrolyte concentration c ∞ and the voltage at the midpoint between clay platelets ψ m (Sposito 1984), recall equation (3): Thus, the swelling pressure as defined in equation (8) is: which indicates an entirely 'osmotic swelling pressure'. If the left and right edges of the tensor box were located infinitely close to the clay platelet surfaces, an identical expression would be found but, in this case, the swelling pressure would appear as an entirely 'electrical repulsive pressure'. In general, if the tensor box edges were located at arbitrary positions x's, the general result would still be the same as equation (10), suggesting that the variations in osmotic pressure are balanced by variations in the electric forces acting in the same planes on the fluid between the clay particles (Bolt 1956). Equation (10) is in fact the closed form solution of the non-linear Poisson-Boltzmann equation (PB) for infinitely long parallel charged particles (Van Olphen 1977). We first use this solution to validate the numerical results for equations (3) and (4). This is done by comparing the solution provided by equation (10) with the numerical result of the following equation (10): Figure 3(a) exemplifies the agreement obtained between the swelling pressure computed using the analytical solution (eq. (10)) and the numerical solutions (eq. (11)), as a function of the half-distance d between two clay platelets (see Fig. 3(a)-inset). The background electrolyte concentration in this example is 10 mol • m −3 , and the clay surface charge is −0.108 C • m −2 . The solution is valid for separation between particles greater than a few nanometres. This limitation reflects: (1) the theoretical limitations of the approach: most of the evidence support that the nonlinear PB equation is accurate for surface potential less than 100 mV and salt concentrations less than 100 mol • m −3 (since the ions are considered as point charges, disregarding that their finite size is more important at very short distances from the clay surfaces) (Van Olphen 1977); (2) given the high non-linearity of concentration and potential profiles near the clay wall, the number of nodes necessary to accurately resolve them makes it impractical. Additionally, we check for consistency of the numerical solutions for different choice of tensor 'boxes' (see Fig. 3(b)-inset) and background electrolyte concentrations and surface charges (Fig. 3(b)). In fact, the divergence theorem shows that force (per unit area) can be calculated by integrating over any surface containing a particle, not just on the particle surface itself or at the midpoint between particles (Phillips 1999). We cover the typical range values for montmorillonitic clays (i.e. σ = −0.108 C • m −2 to −0.266 C • m −2 ).

FINITE SIZE PARALLEL PLATELETS
In the previous section, we validated our numerical models for a range of values of interest, and that are within the limits of applicability of the non-linear Poisson-Boltzmann equation, a particular case of the microscale PNP equations. In this section, we use our numerical models to explore cases in which analytical solutions do not exist.
Consider the following cases: stacked finite size parallel platelets (Fig. 4-left) and 'shifted' finite size parallel clay platelets in a 'stretcher bond' pattern ( Fig. 4-right). All the platelets are 1 nm thick and 100 nm long, the size of typical montmorillonitic clays (Santamarina et al. 2001). The interplay among background electrolyte concentration, surface charge density and swelling pressure (through interacting double layers) will determine the equilibrium separation between particles.
We compute the balancing vertical swelling pressure T sw for given vertical and horizontal separations between clay particles. Different spacings can be interpreted as macroscopic swelling and shrinkage behaviours. The domains of the respective numerical models are represented in Figure 4 by the dotted lines and only compromise pore fluid (sodium chloride solutions) that surround (rigid) 2D clay platelets. The boundary conditions are: symmetry (Boundaries (1) in Fig. 4), prescribed background concentrations and voltages (Boundaries (2) in Fig. 4) located in the far field where no influence of double layer is noticeable and surface charge density on the clay platelets boundaries.   of half-distance between platelets for edge-to-edge distances of 120, 40, 6 and 0 nm (infinitely long parallel platelets); (b) as a function of edge-to-edge separation for 2.2, 2.3, 2.4, 2.5, 3, 3.5, 4, 5, and 10 through 50 nm half-distances.
The figure also shows the limiting case of nil edge-to-edge separation, that is, the previously discussed case of infinitely long parallel platelets. The latter case is indeed the upper bound swelling pressure of this stacked fabric configuration. This bound has been corroborated via Monte Carlo molecular dynamic simulations (Phillips 1999, Sposito 1984 and the exponential decay trends with platelet separation have been also observed experimentally (Bolt 1956, Low 1979, Sposito 1984. Figure 5(b) depicts more detailed information, with swelling pressure as a function of edge-to-edge separations. As the separation between particles increases, the computed swelling pressure decreases. This change is more sensible for face-to-face separations where the influence in the overlapping double layers is larger, that is, when the vertical 'half-distance' between platelets decrease (refer to figure inset).
Figure 6(a) shows the computed swelling pressures in the shifted fabric, for edge-to-edge separations of 40, 20, 10, 6 and → 0 nm (e.g. infinite parallel). A trend similar to Figure 5(a) is observed; however, the comparison of the stacked parallel fabric with the shifted parallel fabric shows that the swelling pressure is larger for the stacked configuration, and that the difference of the increment of swelling pressure increases as the edge-to-edge separation increases. That increment ranges from 15 to 30% at a small horizontal separation (for the cases that we show here), and up to 70% for larger edge to edge separations, again due to increasing overlapping double layers. (Note to readers: the figure-inserts are just indicative of the geometric configuration of particles being considered). To conclude the discussion of these results, we also compute the horizontal balancing force per unit area for both fabrics for the horizontal swelling pressure. In doing so, the tensor box must be placed surrounding the clay particle, making sure that it does not intersect the particle. We can now estimate the horizontal-to-vertical pressure ratios. This ratio represents the well-known 'coefficient of lateral pressure at rest' or K 0 , as at larger face-to-face separations the van der Waals forces become much smaller. In montmorillonitic clays, K 0 ranges typically between 0.45 and 1.1, therefore the region between dotted lines in Figure 6(b) represents the only realistic configurations for these fabrics under the analyzed conditions (c ∞ = 1 mol • m −3 σ = −0.266 C • m −2 , edge-to-edge separations of 6 and 10 nm), other states may be reached upon shear deformation. Estimations of K 0 provide a powerful and relatively simple tool to predict fabrics that are realistic at rest. For example, a face-to-face separation of 20 nm (e.g. 10 nm half-distance) and edge-to-edge separation of 10 nm would probably not be found in nature in a clayey soil deposit under 'at rest' conditions since the estimated coefficient of lateral pressure at rest K 0 for this geometrical configuration does not fall within realistic K 0 ranges of clay deposits.

FINITE SIZE PERPENDICULAR PLATELETS
Finally, we analyze the case of edge-to-face fabrics. Only the average pressure between two particles is explored. We consider an unrealistic low surface charge (for montmorillonite) of −0.010 C • m −2 and a more sensible value of −0.108 C • m −2 . Figure 7 summarizes the results for an edge-to-face repeated configuration. As far as we are aware, we show for the first time the average constant vertical swelling pressure under this flocculated particle arrangement. This figure suggests that the swelling pressure increases as the surface charge of clay particles increases as a result of stronger electrical interactions. The absolute value of the estimated pressures are much lower than the corresponding parallel case, just hinted in Figure 3(b), as a result of smaller effective interacting surfaces (double layers).
As with the previous simulations, caution must be taken in the accuracy of the electrical potential gradient estimation. Higher order shape functions and finer meshing overcomes this problem, but increases the computation demand.

CONCLUSIONS
The following conclusions are drawn from this work: • Particle-scale modelling using coupled microscale Poisson-Nernst-Planck system of equations to determine concentration and electrical potential distributions, together with the van't Hoff relation and Maxwell stress tensor, allows estimating swelling pressures, and could be also used to estimate potential energy for a variety of periodic fabrics, numerically. However, this analysis is limited to diluted montmorillonite (clay) -water (solvent) -monovalent sodium chloride (electrolytes) systems (though this could be extended to more complex salt solutions using molecular dynamics tools).
• High accuracy in the determination of electrical potential gradients is needed, which can be achieved with higher order shape functions in the finite element formulations and finer meshing. Nevertheless, computational demand increases substantially.
• Swelling pressure decreases as edge-to-edge and edge-to-face separations increases. The increment is larger as face-to-face separations increases, due to larger influence in overlapping diffusive double layers. The possibility of estimating vertical and horizontal swelling pressures and their ratio (K 0 ) has been explored and, for some particle spacings, the computational model did yield at-rest earth pressure coefficients that are observed in practice.
Future work will include a more detailed quantitative study of the influence of background electrolyte concentration and the application of this methodology to more realistic flocculated and disperse clay particle arrangements. Extension of this methodology to 3D pore-scale geometries is possible; however, computational demand and definition of complex 3D 'tensor boxes' represent a major challenge.

ACKNOWLEDGMENTS
The authors acknowledge the Australian Research Council's support for this work under the Discovery Project scheme. Partial support was also provided by the University of Melbourne, ECR grant scheme.