Chemo-Mechanical Couplings at Granular Contact: The Effect of Mineral Dissolution and Precipitation across the Scales

: Strong interactions between mechanical deformation and chemical reactions may play a critical role in the response of geomaterials or geological systems to evolving environmental circumstances that may occur in both natural and engineered processes. The present study focuses on mineral dissolution and precipitation at the intergranular contact whose consequences are often manifested at the macro-scale where the mechanical and transport properties of the geomaterial may be altered. Discrete element modeling is employed to explore two applications involving such mineral transformations. The ﬁrst example is primarily focused on the chemo-mechanical coupling mechanisms of intergranular contact in the natural process of pressure solution and secondary compression. The effect of the mineral dissolution on the mechanical response at the grain contact is incorporated into the contact model. Discrete element simulations are performed to examine the overall mechanical response of particle assembles subject to mineral dissolution and the results demonstrate the important role of the kinetic rate characteristics of the dissolution process. The second part of the present study revolves around the effect of mineral precipitation in an engineered process known as microbially induced calcite precipitation for potential soil improvement. The kinetics of involved bio-chemical process is incorporated into on the contact model and the simulation results indicate considerable strengthening effect. Overall, the present study demonstrates the feasibility of discrete element approach as a numerical tool to model coupled chemo-mechanical phenomena across the scales.


Introduction
Micro-scale phenomena of mineral dissolution and precipitation in geomaterials may play a crucial role in the macroscopic processes of mechanical deformation and pore fluid transport. There have been growing interests in the behavior of geological materials or systems that are susceptible to the evolving environmental circumstances which often induce a strong interaction of chemical and mechanical processes. Examples include granular contact dissolution in sediment subsidence [1,2], bonding and cementation in rock weathering [3,4], mineral transformation triggering landslides [5,6], sulfate-induced subgrade expansion [7,8], and karstic rock dissolution in sinkholes [9], to name just a few. There is often a two-way coupling between chemically mediated deformation and mechanically enhanced reaction that may progress in a self-organized manner [2]. Such interactions often lead to the enhanced alteration of porosity, permeability, strength, stiffness, and other material properties.
It is worth noting that such chemo-mechanical interactions are not only ubiquitous in natural processes, but also may be considerably intensified or accelerated in anthropogenic processes enhanced by increased human activities, and in some cases even deliberately taken advantage of in engineering applications to accomplish certain desirable outcomes. Interests in technologies such as clay buffer in nuclear waste disposal [10], geological storage of carbon dioxide [11], and hydraulic fracturing in unconventional oil or gas recovery [12] continue to grow worldwide to seek specific applications of coupled multiphysics and multi-scale interactive processes. A prime example that has emerged in the late 2000s and undergone rapid development in recent years is the microbially induced carbonate precipitation (MICP) technology for potential soil improvement [13]. In this application, microorganisms are introduced as a catalyst to promote the precipitation of calcium carbonate for strengthening and stiffening sandy soils.
Mineral dissolution and precipitation was originally attributed to the process of pressure solution which is considered as a critical mechanism for sedimentary rock formation and evolution by geological society. Its theory, which dates back to the 1950s, focuses mainly on the geochemical aspect and does not address the specific mechanisms of various mechanical processes at the granular contact under typically very high contact stress [14][15][16][17]. However, such phenomena recently attracted strong interest from geotechnical engineering communities who recognized the role of chemo-mechanical granular contact in a number of important geotechnical contexts [18][19][20]. The coupled chemo-mechanical processes at the intergranular contact are investigated for the natural process of secondary consolidation of sediments in Hu and Hueckel [19] and Hueckel and Hu [2], concentrating on a series of multi-physics phenomena through multiple scales based on the extended geomechanics theory of chemo-plasticity [4,21,22]. These studies are primarily focused on the couplings between the chemical reactions and mechanical strength evolution and the involved acrossscale correlations are generally built on the simplified structure of particulate assembly. One of the major challenges encountered in the modeling of multi-physics multi-scale phenomena is the characterization of critical processes at each scale and the homogenization across scales. In the present study, two distinctive phenomena are explored by employing the discrete element modeling to investigate the behavior of particulate assemblies where chemically affected contact mechanisms play an important role. In the first example, examined in Section 2, a natural process of pressure solution and enhanced deformation is studied that may originate from the mineral dissolution occurring at grain contact which eventually affects the overall response of a granular assembly. The second part of the present study is dedicated to the engineered process of microbially induced carbonate precipitation (MICP) presented in Section 3, mainly from a computational perspective to investigate the effect of carbonate precipitation between sand grains. The principal goal of the present study is to examine the feasibility of discrete elements to model sequentially or simultaneously coupled chemo-mechanical phenomena.

Geomechanical Background
The role of mineral dissolution has been long recognized by the geological society as a key mechanism in the process of pressure solution, which is defined as a deformation process resulting from transfer of solid material from highly stressed areas through diffusion in pore solution toward stress-free areas, and considered as one of the main mechanisms for compaction, cementation, and deformation of sedimentary rocks during diagenesis or medium grade metamorphism [15,16,23]. Figure 1a shows a schematic representation of pressure enhanced dissolution and diffusion processes between the granular contact. There are multiple proposed theories about the formation of different micro-structures where the dissolution occurs, such as water film diffusion [14,16], free-face solution [17], and island channel structure [24,25]; however, in general, the intricate mechanical characteristics under the contact are not extensively addressed in the relevant theories. The major interest in pressure solution studies is focused almost exclusively on the reduction in porosity [17].
Gradually, the significance of mineral dissolution is recognized in the study of various chemo-hydro-mechanical coupled processes in geotechnical engineering. Hueckel et al. [1] proposed the mechanism of mineral dissolution in the interpretation of the experimental observations from a series of soil secondary consolidation tests. Pietruszczak et al. [18] investigated the pressure solution of chalks. Cha and Santamarina [20] studied the localized mineral dissolution and formation of open conduits in cohesive media such as carbonate rocks. It is also worth noting that the enhanced deformation based on pressure solution theory is only appreciable after prolonged periods of time, usually in geological scales. However, laboratory soil tests suggest that the time scale of the underlying physico-chemical processes can be much shorter, e.g., in a matter of days or weeks [1]. Therefore, of particular relevance to these geotechnical investigations is the extension of the chemo-plasticity theory to model the coupled processes of mineral dissolution and deformation across multiple scales proposed by Hu and Hueckel [19], based on experimental evidences reported by Denisov and Reltov [26] and Baxter and Mitchell [27]. As shown in Figure 1b, at the micro-scale, the near-contact damage and irreversible deformation at the intergranular contact accelerates the mineral dissolution by the creation of free specific dissolution surface area which results from dilatant damage (microcracking) induced by the plastic yielding and, subsequently, the dissolved mineral diffuses within the grain away from the reaction sites. The consequences of such microscopic processes are manifested in the upper scales; at the meso-scale, intergrain diffusion of the dissolved minerals occurs in the pore water and subsequently the dissolved minerals precipitate on the free surface of adjacent or remote grains. It is worth noting that these phenomena can be complicated by a number of physical and chemical processes at different scales: the difference between the lithostatic and the hydrostatic pressure as an important part of the driving force behind the pressure dissolution can considerably affect the solution rate at the granular contact; accumulation of precipitated mineral around the free surface of the grains could potentially block or influence the diffusion path; evolution of the mineral concentration in the pore fluid may also strongly affect the diffusion process [28].
These effects eventually lead to not only enhanced compaction, i.e., reduction in pore space, but also alteration of other mechanical and hydraulic properties including stiffness and strength increase and permeability decrease at the macro-scale.
The multi-scale model discussed above typically demands that quantities at the macroscale need to be derived through cross-scale transfer functions (e.g., volume averaging). Indeed, the macroscopic representative elementary volume (REV) considered in the above mentioned work [2,19] is of prototype topology which adopts a periodic structure of grain assembly of uniform particle size to render relevant semi-analytical formulations for altered soil properties. This naturally leads to the consideration of employing an effective computational tool such as discrete element modeling (DEM) [29] to simulate a realistic assembly of a wide range of particle sizes.  [30,31]. (b) Schematic representation of an enhanced dissolution mechanism that considers the near-contact localized deformation and damage [19]. (c) DEM configuration of particle contact with geometry parameters indicated for subsequent modeling considerations.
In the present study, a degradation based contact model is proposed to address the effect of mineral dissolution on material weakening which, in turn, affects the local mechanical response in terms of force or stress, and displacement or strain ( Figure 1c). Obviously, the basic discrete element contact models are not capable of dealing with coupling effects and hence need to be modified in the numerical implementations. Still, employing a DEM scheme imposes certain restrictions on the flexibility of the possible coupled formulations, as the number of geometric parameters and the details of granular contact are rather limited in typical DEM contact models which, however, should not be perceived as a mere weakness of such an approach. The simplicity of its contact models at the particle level is also its main strength, making this numerical tool computationally affordable. Compromises have to be made to allow adequate representation of characterization of critical mechanisms of micro-mechanical responses and at the same time meet computational demands arising from analysis of large quantities of particles.
The proposed contact formulation needs to be simple but must be sophisticated enough to adequately represent the critical interactive chemo-mechanical mechanism. In DEM analysis, particles are often assumed rigid but can overlap with each other. The overlap is analogous to the deformation that occurs between real particles. It can be described by a contact-force formulation. Typical contact models, e.g., a linear model or Hertz contact model, can be used to describe the linear or nonlinear displacement response to the contact force at the contacts. In the present study, a linear contact model is used. The relationship describing the contact interaction is characterized by the normal contact stiffness and the shear contact stiffness. In the present study, a degradation effect is introduced to the linear contact model to account for the chemical effects on the granular contact responses. For the sake of simplicity, both normal and shear modulus (k) are modified in the same formulation as follows. (1) is formulated as a function of the solid mass removal parameter from the system, ξ, by adopting an exponential function form [32] to imposes a rational range for mechanical degradation, α, β 1 , and β 2 are constants; as such, regardless of the values of the three constants, Φ is always constrained in the following range: 0 ≤ Φ ≤ 1, which also represents the progression of mechanical degradation confined in a sensible range. Degradation functions described by Equation (2) are depicted in Figure 2 under different sets of parameters to show the progression of the degradation as a result of mineral dissolution. The condition α = 0 or α = 1 deactivates the second or third term of Equation (2), and reduces the needed parameters to two. The evolution of mass removal ξ is dependent on the progression of mineral dissolution typically described by the kinetic rate laws as discussed in what follows in Section 2.2.

Kinetic Rate Laws for Geochemical Reactions
In the present study, the chemical variable ξ is proposed to describe the measure of mass dissolution of the mineral that weakens the material; it is taken as the relative mass removal, defined as the ratio of the dissolved mass (or mol) of mineral, M diss , to the original mass (or mol) of mineral, M 0 . Hence, it is also confined in the range of [0, 1], when ξ = 1 represents the complete depletion of the mineral. Its rate can be described as The time rate of the dissolution of solid mineral is typically described by a specific geochemical reaction rate law. In the presented study, two types of dissolution rate formulations are investigated, the first being based on the established theories of pressure solution processes and the other based on a damage-enhanced kinetic rate law. It is also noted that the presented example in this section is limited to the analysis of mineral dissolution only and does not consider the precipitation process, for the sake of simplicity, although it is entirely possible to include the mineral precipitation in the relevant formulations. The potential effect of precipitation is explored in the example discussed in Section 3.
The first type of mathematical formulations explored is based on theories of pressure solution where various forms of mathematical formulations [17,30,33] have been proposed with the central idea of directly linking the dissolution rate to the granular contact pressure or stress; in the present study, the following equation, as proposed by Yasuhara et al. [30], is used, where R is the ideal gas constant and T is the absolute temperature; ν m is the molar volume of the solid and k + is the dissolution rate constant; ρ g is the solid grain density, and d c is the diameter of the grain-to-grain contacts, from which contact area, A con , as indicated in Figure 1c, can be calculated. The effective stress σ eff can be readily computed as the average contact stress σ con = F n /A con in the numerical implementation, where the contact force F n at each contact is computed at each time step. Hence, the rate of ξ can be expressed aṡ where the average diameter of the two grains in contact,d = 0.5(d 1 + d 2 ), is used for the grain volume calculation; it is noted that the original mineral mass of the solid grain is calculated as M 0 = πd 3 /6 in the above equation. As such, the rate of dissolution is computed, and the corresponding dissolution and degradation are updated at each time step in the DEM simulation, depending on the contact normal force. In this formulation, the dissolution rate is strongly affected by the contact stress and hence referred to as the stress affected formulation hereafter. The second type of mathematical formulations explored in the present study takes into account the near-contact damage as the major factor. For a specific mineral-water reaction, e.g., silica-water reaction, a rate equation was proposed in Hu and Hueckel [19], which incorporated the irreversible plastic strain into the specific surface area in the kinetic rate law established by Rimstidt and Barnes [34] for silica-water interaction. Such more elaborate formulation is simplified in the DEM scheme where only simple geometric parameters are used. The overlap distance between the two grains, d over (Figure 1c), is used to describe the progress of strain development. In addition, the rate law proposed in Hu and Hueckel [19] is also simplified to yield a relative dissolution removal quantity with respect to the mole of the original mineral where a represents the activity of each reaction species, and d 0 describes the initial distance of the two grains just in contact, i.e., d 0 = d 1 + d 2 . The functionĀ represents the relationship between the dissolution reaction rate and the accumulated local deformation, characterized via a parameter, Γ, which entails the consideration of the opening of the micro-cracks around the near contact for specific surface area where dissolution occurs; however, in the present DEM scheme such local deformation is simply represented by the overlapping distance d over (Figure 1c) in the above equation. It should be noted that the overlapping distance is calculated at each time step, as the displacement response at each granular contact evolves with time when the contact is affected by the degradation as a function of mass dissolution. In this formulation, the dissolution rate is coupled with the contact deformation and hence referred to as the deformation affected formulation hereafter.

Numerical Results
The presented coupled chemo-mechanical model is implemented in the discrete element analysis package PFC2D [35] using its programming language FISH to simulate the enhanced deformation affected by grain contact dissolution of two different sample minerals, mimicking "quartz" and "calcite", the former characterized with a slow kinetic dissolution rate and the latter with a relatively faster dissolution rate. Each particulate assembly is prepared with an initial void ratio of 0.54 and a height-to-width ratio of 2:1, consisting of approximately 2500 round particles; the range of the particle size is between 0.75 and 1.1 mm. The normal and shear contact stiffness are 5 × 10 8 N/m and the inter-particle friction is 0.23.
For the sake of simplicity and comparison, the same forms of kinetic rate formulation discussed earlier for quartz-water reaction are used to describe the calcite dissolution, with the rate constant acting as the dominant controlling parameter. Some key parameters used for the simulation are: α = 0.25, β 1 = 0.252, β 2 = 1.02; ν m = 22.7 cm 3 /mol and 31.2 cm 3 /mol for quartz and calcite, respectively; k + = 3.98 × 10 −14 mol/m 2 /s and 1.55 × 10 −6 mol/m 2 /s for quartz and calcite, respectively, at the considered temperature T = 298.15 K. In each simulation, a typical biaxial compressive loading is applied to the grain assembly; the assembly is first brought to an initial isotropic stress state of MPa under zero friction condition. Subsequently the additional axial loading is applied under the constant confining stress of 1 MPa with friction coefficient activated while each grain is subject to enhanced mass dissolution at the contact. Figure 3 shows the stress-strain responses of the DEM samples. Geochemical reaction occurs simultaneously as the biaxial test is run. The strength decreases as a result of the ratedependent mineral dissolution process. Both formulations as described by Equations (4) and (5) are examined in Section 2.2. Both formulations show a similar general trend regarding the responses of samples subject to mineral-water reaction. Figure 4 shows the corresponding volumetric deformation. It suggests that the greater dissolution and degradation in calcite assembly seems to lead to the less dilative behavior. However, although much faster calcite dissolution appears to lead to a lower strength response, the differences in the presented figures appear to be rather small. This observation is due to the fact that the linear degradation model incorporated into the contact model appears to be too modest. As discussed in Utili and Nova [36], the contact stiffness needs to be altered by orders of magnitude (i.e., logarithmic scale) to render significant differences in DEM simulations, even though physical evidence or explanation for such drastic degradations remains to be explored. In addition, dissolution rate is very small and its effect is barely appreciable unless for a significantly prolonged period of time. It is of interest to examine the evolution of the degradation parameter averaged over the assembly as shown in Figure 5. The main difference between the two stimulated minerals is clearly due to discrepancy in the kinetic rates.   It is worth noting that, just as with most DEM simulations, the presented analysis should be considered as a numerical exercise rather than precise predictions of actual material response. A natural sand may consist of a wide range of particle sizes; the effect of particle size distribution may affect not only the overall simulated responses, but also the mechanical contact and chemical reactions between particles. In addition, as PFC2D enforces neither plane stress nor plane strain condition, such simulations should not be aimed to replicate responses in triaxial tests. However, this numerical study demonstrates that it may be feasible to adopt a DEM approach to link coupled processes from grain contact to the macroscopic material response. It should be also noted that early investigations on time-dependent DEM simulations were primarily focused on the soil creep behavior and largely achieved using explicit time-dependent constitutive functions without addressing any specific mechanism [37][38][39], the present study shows that it is possible to incorporate a intrinsic rate-dependent process or mechanism such as chemical dissolution into DEM simulations.

Mechanisms and Mathematical Formulations
The second example examined in the present study is an engineered process that takes advantage of the chemical precipitation for potential soil improvement. Microbially induced calcite precipitation (MICP) is a promising soil improvement technique that exploits the bio-chemical processes to induce chemical precipitation of calcite in granular soils, thus reducing their porosity and permeability and improving their strength and stiffness [13,[40][41][42][43][44][45][46]. Despite extensive research and successful implementations achieved in the laboratories, there are still a number of open questions remaining regarding the predictability and controllability of engineered soil properties, which are critical for potentially successful application of this technology in the field scale. As the involved hydro-biochemo-mechanical processes are inherently inter-connected, and these interactions often occur across multiple scales, an in-depth study of MICP would likely entail a multi-scale character. The spatial distribution and temporal evolution of bio-mediated precipitation and the altered soil properties must be understood through the key mechanisms and their interactions at and across each scale. In the present study, numerical simulations with DEM are investigated to assess the macroscopic strength behavior of treated soils by implementing the grain and pore scale phenomena into the meso or macro-scale grain assembly simulation.
The fundamental idea of MICP is the engineered promotion of calcium carbonate precipitation in the presence of micro-organisms such as ureolytic bacteria. The cementation solution of calcium chloride and urea is supplied into the soil in the presence of ureolytic bacteria (e.g., Sporosarcina pasteurii), which is typically introduced into the soil prior to the cementation solution. The hydrolysis of urea produces ammonium and carbonate ions, and the latter react with calcium ions to produce precipitated calcium carbonate; these two processes can be described by the two following reactions, The key to the effectiveness of this technology is the successful precipitation of calcite and the resulting pore reduction and particle bonding. However, it is worth noting that, at the grain or pore-scale level, some calcite mineral might precipitate and reside around a sand particle, forming a "uniform" coating which may not be effective in improving the soil property [13]; whereas if the precipitated mineral stays between particles, forming a bridging contact, this would be more effective in improving the soil strength and stiffness. It is uncertain whether certain control or influence can be exerted to induce a preferable outcome of the mineral precipitation at the grain scale, and much study is needed to explore such possibilities.
The present study is aimed to employ numerical techniques to assess or simulate the upscaling of strengthening effect of calcite precipitation. The discrete element method provides a useful tool to simulate the response of a grain assembly affected by calcite precipitation [47,48]. However, in these earlier efforts the time-dependent progress of calcite precipitation itself, which is obviously affected by the rate of chemical reactions (i.e., Equations (6) and (7)), was not considered. In the present study the progression of calcite precipitation is modeled with a Monod kinetic rate formulation; the reaction rate R reac for calcite precipitation is assumed to take the following form as suggested by van Wijngaarden et al. [49] where U m is the maximum reaction rate, K m is the half saturation constant, and C urea represents the concentration of urea. The term in the parenthesis represents the decay of the micro-organisms, describing the reduction in or loss of activity of bacteria; τ m is the maximum life for urease to be active. It is noted that other forms of decay functions such as an exponential function are also possible [50,51]. Based on the above calcite precipitation rate formulation, the following equation for the mass of calcite can be derived [49], C CaCO 3 (0) is the initial mass of calcite per unit volume, θ 0 is the initial porosity, ρ CaCO 3 and m CaCO 3 are the density and molar mass of calcite, respectively. It is noted that this derivation assumes a constant concentration of urea, i.e., spatial variation or transport of urea is not considered, hence the adopted formulation (Equation (9)) reflects a closedsystem approximation of calcite evolution. Calcite may precipitate either between granular contact or around the free surface; in the present simulation, only the former is considered; thus the contribution of precipitated calcite is reflected by the growth of the bonds. The bond stiffness and radius evolve as a function of calcite precipitation which itself progresses with time as described by Equation (9). The adopted formulation is based on experimental observations reported in Al Qabany et al. [52] which relates the stiffness increase and the shear wave velocity increase to the precipitated calcite. For simplicity, the bond radius (R b ), normal stiffness (Kb n ), and shear stiffness (Kb s ) are assumed to evolve as described in the same manner: The variable with subscript "0" indicates the original value of R b , or Kb n , or Kb s ; each parameter is represented by a generic symbol B. λ represents a parameter indicative of the effect of the calcite precipitation and is taken to be 0.0638 m 3 /kg based on the calibration of the data reported in Al Qabany et al. [52]. As Equation (9) is adopted for all contacts, calcite precipitation is also assumed to be spatially uniform in the assembly when, at each contact, the bond starts to grow at the same rate despite the fact that the initial bond radii are varied.

Numerical Simulations
The above formulations are implemented in PFC2D [35] with its programming language FISH. The same sample preparation and DEM parameters for the first example reported in the preceding section are used here for a similar biaxial loading test. As noted in Section 3.1, the assembly additionally features a parallel bond with both the initial bulk (Kb n ) and shear modulus (Kb s ) being set to be 100 GPa. The following constants are used for the calcite mineral, n 0 = 0.35, C CaCO 3 (0) = 0, ρ CaCO 3 = 2710 kg/m 3 , m CaCO 3 = 100 g/mol, K m = 1 × 10 −5 mol/L, C urea = 1 × 10 −3 mol/L, and τ m = 6.12 × 10 5 s. Figure 6 shows the simulation of a typical response of bio-treated sands in a deviatoric stress vs. axial strain curve; the trends of loose and dense sands such as behavior are strongly evident and agree well with the laboratory experimental observations reported in the literature [13,40]; it also shows the significant strengthening effect of the bio-treatment of mineral precipitation. A few comments about the implementation of the mineral precipitations are due. It should be noted that a typical DEM biaxial loading consists of two stages, an isotropic loading (under a certain confining pressure) followed by additional axial loading (shearing), as discussed in Section 2. In the first stage, calcite precipitation is kept to be progressing for a certain period of time before additional axial loading is applied. While precipitation occurs throughout the entire process, the strengthening effect is largely caused by the accumulated calcite produced during the first stage. This effect of accumulation of calcite can be also explored for longer periods of time as shown in Figure 7. The evolution of strength behavior is clearly affected by the growth of calcite accumulated with time. The presented simulation shows the potential of numerical tools for upscaling the evolution of calcite precipitation at the granular contact scale to its effects at the macroscopic scale.

Concluding Remarks
This present study explores a discrete element approach to modeling coupled phenomena of mechanical deformation and intrinsic chemical reaction processes. The first example investigated is a typical natural process of multi-scale interactions of mineral dissolution and mechanical deformation. It examines a modified linear contact model that describes the degradation effect caused by mineral dissolution, which evolves with time with its rate being affected by the local grain contact force or contact displacement. Both stress and deformation controlled formulations are used to describe the kinetic rate laws for enhanced dissolution. Simulations of quartz and calcite minerals, characterized by different rate constants and physical properties, indicate different softened stress-strain responses.
A second example explored is focused on the engineered process of bio-mediated mineral precipitation for potential soil improvement. It is different from the first example as it involves only a one-way coupling between chemical and mechanical processes in the absence of the mechanical effect on the chemical reaction. Hence, its formulation of kinetics of the involved bio-chemical process is straightforward and its effect on the mechanical properties of granular assembly can be conveniently assessed via DEM simulations. Simulations of this type are confined to closed systems and future investigations may seek innovative ways to incorporate the consideration of diffusion transport processes. The results show the considerable strengthening effect of mineral precipitation as observed in laboratory experimental studies. However, to better understand the intricacies of coupled processes across multiple scales, collective efforts in experimental and numerical investigations are still very much needed. Overall, the present study demonstrates the feasibility of DEM approach as a numerical tool to study the coupled processes at the grain contact.