An evaluation of contact models for particle-scale simulation of clay

Geotechnical engineers are well aware that the particle surface chemistry and the pore fluid composition can significantly influence the mechanical behaviour of clay. Reference is often made to the Derjaguin-Landau-Vervey-Overbeek (DLVO) theory, which enables the electrochemical interactions between charged particles to be estimated. Hitherto, the absence of an effective framework for particle-scale simulation of clay has inhibited a direct link between these electrochemical interactions and clay behaviour (e.g. load:deformation response) or fabric (i.e. the development of a disperse or flocculated fabric). Ebrahimi [1] demonstrated the viability of using molecular dynamics simulations where the clay grains are simulated as ellipsoidal particles whose interactions are described by an analytical expression called the Gay-Berne (GB) potential. While promising when compared to other approaches documented in the literature, Ebrahimi’s work considered only a single clay mineralogy and did not explicitly account for the pore fluid composition. This paper considers the use of the Gay-Berne potential in particle-scale modelling of clay from a more general perspective. Calibration of the GB model parameters to predict kaolinite particle interactions reveals a lack of generality in Ebrahimi’s approach. The Gay-Berne potential cannot simulate situations in which attractive and repulsive interactions co-exist, which lead to the classical “cardhouse” fabric, as is the case of kaolinite particles interacting via an acidic pore fluid.


Introduction
Inter-particle forces between clay particles can be categorised as contact and non-contact forces (e.g. [2]). The non-contact forces are electrochemical forces that combine the contribution associated with the van der Waals attraction ( ) and the electrostatic ( ) interaction. These non-contact forces are understood to dominate the interaction between clay particles ( [3], [4]). It is generally accepted in soil mechanics that the electrochemical forces can be described by the DLVO (Derjaguin-Landau-Vervey-Overbeek) model ( [2], [5]). This model was developed in the 1950s to explain the equilibrium behaviour of colloids in solution and describes the variation of the potential energy of the particle-particle interaction with separation distance.
Ebrahimi et al. [6] and Carrier [7] highlighted limitations of the ability of DLVO theory to accurately describe the interactions between clay particles. They argue that DLVO theory cannot accurately capture the stability of clay-water systems at short distances, where aggregates form. However, given its general acceptance in the geotechnical literature, it is important to have an effective means to incorporate this theory in particle-scale models of clay. Existing approaches to model clay at the particle scale that have considered DLVO theory include the contribution by Liu [8], who simulated the interaction between particles using the Discrete Element Method (DEM). Sjoblom [9] used molecular dynamics (MD) and employed a contact model containing both DLVO and contact forces that exist between charged particles at very short separation distances. Pagano et al. [10] proposed a 2D DEM model with a simple tri-linear contact formulation that considers the DLVO forces in a highly ideal manner and simulates also particle contact. In each of these contributions, a single kaolinite clay particle was modelled as an assembly of spherical sub-particles clustered together to form a rigid body. A large number of sub-spheres is required to create a model clay particle with a realistic aspect ratio. However, the use of a large number of sub-particles to describe a single particle is computationally costly and may inhibit achieving a representative element volume (REV) in MD or DEM simulations.
Ebrahimi et al. [6] developed a more computationally efficient approach using the MD code LAMMPS ( [11]). They modelled clay particles as ellipsoids and captured the interaction between them using the Gay-Berne (GB) potential. The GB potential describes the variation in the potential energy of the particle interactions as a function of separation for rigid ellipsoidal particles ( [12], [13]). By appropriate choice of semi-axes lengths, it is reasonable to approximate a clay particle as a flat ellipsoid. The approach proposed by Ebrahimi et al. [6] is significantly more computationally efficient than the clustered-sphere approach favoured by other authors, as discussed above. Ebrahimi et al. [6] chose their GB model parameters by fitting the energy-separation data obtained from atomistic MD simulations rather than using DLVO theory. They considered only one type of clay (montmorillonite) and did not explicitly account for the pore fluid composition.
This paper builds upon the contribution of Ebrahimi et al. [6] by assessing the general applicability of their approach to simulate clay particle interactions. The importance of calibrating the GB potential to simulate the particle interactions predicted by DLVO theory is firstly outlined. The selection of the DLVO model parameters used here is then discussed, prior to presenting a critical analysis of the ability of the GB potential to model kaolinite particles interacting through potassium chloride (KCl) solutions with a range of pH values.

Modelling particle interactions
DLVO theory can be used to derive expressions for the variation of interaction energy with separation distance either for infinitely long planar surfaces or for spheres. In either case the total energy of interaction is given by: The total energy is defined relative to the energy at infinite separation; positive energy values indicate a repulsive interaction, while a negative energy value indicates attraction. The negative slope of the energy-displacement curve gives the interaction force. Considering two infinite parallel plates, the van der Waals energy per unit area [ / 2 ] is defined as (e.g. [5], [14]): where [ ] is the Hamaker constant (which depends on the mineralogy of the interacting faces and the interacting medium), ℎ [ ] is the separation distance between the surfaces, and 1 [ ] and 2 [ ] are the thicknesses of the two interacting plates. This energy is negative, and therefore attractive, at all separation distances, as can be observed in Fig 1 (a)   The strengths of both and are influenced by the pore fluid chemistry (i.e. pH and salt concentration ), which affects the values of in Eqn. 2 and and in Eqn.3, respectively. Table 1 lists the DLVO parameters used to evaluate the energies presented in Fig.1.
The DLVO expressions given in Eqns 1-3 describe the interactions between two infinite parallel plates and so they cannot be directly applied to describe the interaction of clay platelets with a finite size that can interact at arbitrary relative orientations in 3D. Consequently, an analytical framework that can both capture the energyseparation relationship predicted by DLVO and simulate particle rotations is required.
The Lennard-Jones (LJ) potential ( [16]), which is widely used in colloidal science, can describe the interaction of spherical particles and capture the general shape predicted by DLVO theory. Integration of the LJ potential for ellipsoidal particles gives the GB potential, analytically defined as ( [13]): where [ ] is the energy scale, [ ] is the atomic interaction radius, ℎ 12 [ ] is the closest distance between particles and is related to their size and orientation and is the shift of the potential minimum.
The dimensionless quantities and are the shape and the energy anisotropies, respectively. The GB model parameters do not have a clear physical meaning and, therefore, cannot be directly linked with the clay surface chemistry and the pore fluid chemistry. Rather, the model parameters should be calibrated against data describing the energy-displacement relationship for a pair of particles. It is assumed that the intermediate orientations are correctly predicted by the GB model once the face-to-face and edge-to-edge interactions are calibrated correctly ( [6], [12]). The GB potential does not capture the short range attraction (attraction at short distances, Fig. 2) and so it is assumed that the energy threshold noted on Fig. 2 is not crossed. Here the GB model parameters were calibrated using the energy-displacement relationships predicted from DLVO theory at separation distances greater than the separation at which the energy threshold is observed.
Ebrahimi et al. [6] calibrated their GB potential parameters using atomistic MD simulations.

Kaolinite particle interactions
Kaolinite was chosen as the clay mineral to consider here as it is very common, data are available considering the influence of pore-fluid chemistry on kaolinite behaviour (e.g. [17], [18]) and recent experimental studies have directly measured key surface properties ( [8], [15], [19]). The sphere-cluster based DEM/MD studies noted above all considered kaolinite. Kaolinite particles are pseudo-hexagonal thin plates. The basal unit comprises an aluminium octahedron overlaying a silicon tetrahedron ( [2], [20]).
The different charges on the kaolinite particle surfaces are responsible for the interaction between particles in an electrolyte suspension ( [21]). Discussions in the literature on kaolinite particle interactions have considered the particle faces (silica and alumina) and edges. For many years, the face surface (silica and alumina) charges, measured using macroscopic techniques (e.g. potentiometric titration or electrophoresis), were believed to be negative regardless of the pH of the electrolyte solution and were attributed only to isomorphous substitutions ( [19], [22]). More recently, atomic force microscopy (AFM) has been used to measure the charge on each face. This technique has been applied to study the variation of the surface charge with the pH of the pore fluid ( [15], [19]) and with the salt concentration ( [8]). These data show that the silica faces are always negatively charged, while the alumina faces are positively charged for low pH values (< 6) and negatively charged otherwise ( [5], [8], [19]).
The surface charge on the particle edges has always been considered pH-dependent and mainly due to the protonation/deprotonation processes occurring in the hydroxyl groups ( [5], [22]). Uncertainties remain regarding the charge of the clay edges, as their surface areas are much smaller than those of faces, thus complicating the use of AFM to measure the charge ( [5]).

Table 2. Type of interaction (from [15]). Attraction is indicated with A, while R is used for Repulsion
One of six types of interaction can exist between two kaolinite particles: alumina face-alumina face (Af-Af), silica face-silica face (Sf-Sf), alumina face-silica face (Af-Sf), alumina face-edge (Af-E), silica face-edge (Sf-E) and edge-edge (E-E). Table 2 shows how the sign of the interaction (e.g. attractive/repulsive) changes with pH for each scenario. At high pH values (e.g. 8 and 10), the interaction is always repulsive, whereas at low pH values the interaction can be attractive or repulsive depending on the relative orientation between particles.
The variations in particle interaction noted in Table 2 have an influence on the load:deformation response ( [17], [18]). Using scanning electron microscope (SEM) images, Pedrotti and Tarantino [18] showed a sensitivity of the clay fabric to the pore-fluid pH, which can also be linked to the variation in particle interactions noted in Table 2. A particle-based model that aims to consider fabric effects on clay behaviour must therefore be capable of capturing interactions that are always repulsive and also interactions that can evolve from repulsive to attractive as particles rotate relative to each other.

GB Model Calibration
The DLVO expressions given in Eqns 1-3 and the GB potential in Eqn 4 were implemented in a MATLAB model which was verified using data available in Gupta et al. [5] and Ebrahimi et al. [6]. Two ellipsoidal kaolinite particles of diameter = 600 and thicknesses 1 = 2 = 11.2 ( Fig. 3) were considered to compare in the initial analysis of the calibration process. These particle dimensions are taken from Gupta et al., [5] who analysed scanning electron microscope (SEM) images of kaolinite particles. The calibration considered interactions in solutions with a salt (KCl) concentration = 1 at pH=4 and pH=8. Referring to Table 2 above, these different pH values were selected to enable consideration of a scenario where both the face-face and edge-edge interactions are always repulsive as well as a scenario where attractive face-face interactions co-exist with repulsive edge-edge interactions. Six GB parameters, , , , 1 , 2 and 3 , were calibrated to capture the energy:separation data predicted by DLVO in a trial and error process informed by a parametric study. Two interaction directions were considered: faceto-face and edge-to-edge, shown in Fig. 4 (a) and (b), respectively. In both cases, the interactions energies were computed for 0 nm ≤ h ≤ 60 nm.

Case A: pH=8
Fig . 5 illustrates the interaction energy (per unit area):separation distance relationship predicted by DLVO theory at pH=8 for both alumina face-silica face (Af-Sf) and edge-edge (E-E) interactions. In both cases, the total energy of interaction is positive for all separation distances considered and so the interactions examined are always repulsive. The calibrated GB model parameters for pH=8 are shown in Table 3. Fig. 6 compares the DLVO energies with the GB potential data. As expected, the GB potential does not reproduce the attraction that characterises the DLVO model at small distances. However, relatively good agreement is attained at h values exceeding h at the energy barrier. The data on Fig. 6 therefore support the hypothesis of Ebrahimi et al. [6] that the GB potential can effectively model clay particle interactions.

Case B pH=4
The DLVO prediction of the interaction for the lower pH value of 4 is illustrated in Fig. 7. In this case, only the edge to edge interactions are repulsive; the face to face interactions are attractive as opposite charges are carried by alumina and silica surfaces. These contrasting interactions cannot be captured by the GB potential using a single set of model parameters. Table 4 lists the optimal GB model parameters for each interaction type. Referring to Fig. 8, when compared with the DLVO model, use of the GB potential with the Set 1 and Set 2 parameters gives a good fit to the face to face and edge to edge interactions respectively. However, there are clear differences in the parameter values. Referring to Fig. 9 (a), Set 1 parameters give a very poor approximation to the edge to edge interaction; similarly, Set 2 parameters give a very poor prediction of the face to face interaction ( Fig. 9 (b)).

Conclusions
Particle-scale simulations have the potential to advance the understanding of clay behaviour. Ebrahimi et al. [6] proposed modelling clay particles as ellipsoids using the Gay Berne potential. Recognising the computational benefits of this approach, this contribution has examined whether it can be applied to clays in general, i.e. it has extended consideration of the model beyond the montmorillonite particles considered previously and explicitly considers the pore-fluid chemistry. Kaolinite particles were considered with interactions predicted by DLVO theory. The following observations can be made: 1. The GB potential is reasonably successful at capturing the energy:separation relationship considering the higher pH (=8) condition examined here, where both the face to face and edge to edge interactions are repulsive. 2. For the low pH (=4) value considered, attractive face to face interactions co-exist with repulsive edge to edge interactions. In this case the GB potential can capture the face to face and edge to edge interactions only if different parameter sets are used for each interaction type. This is not viable; to use the GB potential within a MD code a single parameters set should be capable of modelling all interactions and capturing the variation in interactions as particles move and rotate relatively to each other. 3. SEM shows that the material fabric will vary with pH and so particle-scale DEM or MD simulations that consider the link between fabric and mechanical behaviour need to be able to capture the pH sensitivity of particle interactions observed in kaolinite.
The data presented here indicate that, while the GB potential has some advantages from the perspective of multi-scale modelling of clay, it cannot be used without modification. Future research should explore alternative modelling strategies. The success at modelling the purely repulsive cases and the fact that, when decoupled, the GB can capture both attractive and repulsive interactions supports the development of a modified GB potential to achieve effective particle scale simulations of clay. This research is funded by the Leverhulme Trust.
(a) (b) Fig. 9. Inability of GB potential to catch face-face and edge-edge interactions at low pH values with a single set of parameters