Abstracts

PHYSICAL SCIENCES

Estimating vertical and lateral pressures in periodically structured montmorillonite clay particles

Guillermo A. Narsilio^I; David W. Smith^II; Peter Pivonka^III

^IDepartment of Civil and Environmental Engineering, The University of Melbourne Civil Engineering Block, D 317, Parkville, VIC 3010, Australia

^IIFaculty of Engineering, Computing and Mathematics, The University of Western Australia 35 Stirling Highway, Crawley, WA 6009, Australia

^IIISchool of Computer Science and Software Engineering, The University of Western Australia 35 Stirling Highway, Crawley, WA 6009, Australia

^{Correspondence to} Correspondence to: David W. Smith E-mail: davidw.smith@uwa.edu.au

ABSTRACT

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.

Key words: clay, Maxwell tensor, Poisson-Nernst-Planck equation, pore-scale, swelling, van't Hoff law.

RESUMO

Dada uma argila montmorilonítica de alta porosidade e saturada por counteríons monovalentes, investigamos as respostas da argila ao nível de partículas para diferentes cargas externas. Como soluções analíticas não são possíveis para arranjos complexos de partículas, empregamos modelos computacionais micro-mecânicos (baseados na solução das equações de Poisson-Nernst-Planck), utilizando o método de elementos finitos, para estimar counteríons e distribuições de potencial elétrico para partículas em diversos ângulos e distâncias uma da outra. Nós então calculamos as pressões de separação usando a relação de Van't Hoff e a tensão de cisalhamento de Maxwell. À medida que a distância entre as partículas de argila diminui e as duplas camadas se sobrepõem, a concentração de counteríons nos microporos entre as partículas de argila aumenta. Este aumento reduz o potencial químico do fluido nos poros e cria um gradiente de potencial químico no solvente, que gera a chamado pressão 'osmótica' ou de 'separação'. Devido a esta pressão de separação, as partículas não precisam de contato entre si, a fim de exercer uma 'tensão efetiva'. Este trabalho pode conduzir a previsões teóricas da resposta macroscópica a carga de deformação em solos montmoriloníticos baseado na modelação micromecânica das partículas.

Palavras-chave: argila, modelo de Maxwell, equação de Poisson-Nernst-Planck, escala de poro, inchamento, relação de van't Hoff.

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 and Soga 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 self-diffusion 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, ε₀ = 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., σ = 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 Solutions in which the activity coefficients are equal to one. , van't Hoff's Law relates the osmotic pressure Π to the electrolyteconcentrations c_i:

where x₁ and x₂ 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 ψ₀ 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):

where the double layer thickness is the background electrolyte concentration (which is the same for both ionic species in 1:1 electrolytes fulfilling electroneutrality), and ψ₀ 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 (Anandarajah and Chen 1994, Verwey et al. 1948):

where for simplicity, a dimensionless potential is defined as >= zF ψ/(RT), thus ₀ = zFψ₀ /(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₁ = x → - ∞, then as E_xx(x₁) = 0; which is also the case if the right edge of the box is situated at x₂ = 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 limitationsof 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.

Figure 5(a) shows the computed swelling pressure as a function of half-distance (e.g. the half of the vertical separation or 'distance' between platelets) for different edge-to-edge separations (e.g. horizontal separations; Refer to figure inset).

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 indicativeof 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₀, as at larger face-to-face separations the van der Waals forces become much smaller. In montmorillonitic clays, K₀ 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₀ 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₀ for this geometrical configuration does not fall within realistic K₀ 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₀) 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.

Manuscript received on April 30, 2008; accepted for publication on March 26, 2009

Selected paper presented at the IUTAM Symposium on Swelling and Shrinking of Porous Materials: From Colloid Science to Poromechanics - August 06-10 2007, LNCC/MCT.

Publication Dates

History