Advanced Geomechanical Model to Predict the Impact of CO₂-Induced Microstructural Alterations on the Cohesive-Frictional Behavior of Mt. Simon Sandstone

Abstract

We investigated the influence of CO${}_2$-induced geochemical reactions on the cohesive-frictional properties of host rock within the context of CO${}_2$ storage in a saline aquifer and focused on the Mt. Simon sandstone. The research objective was to model geo-mechanical changes due to host rock exposure to CO${}_2$-saturated brine while accounting for heterogeneity, double-scale porosity, and granular structure. We formulated a three-level multi-scale model for host rocks. We conducted scanning electron microscopy analyses to probe the microstructure and grid nanoindentation to measure the mechanical response. We derived new nonlinear strength upscaling solutions to correlate the effective strength characteristics and the macroscopic yield surface to the micro-structure at the nano-, micro-, and meso-scales. Specifically, our theoretical model links CO${}_2$-induced microstructural alterations to a reduction in the size of the yield surface, and a drop in the value of the friction coefficient. In turn, regarding the Illinois Basin Decatur Project, the CO${}_2$-induced drop in friction coefficient is linked to an increase in the risk of fault slip and a higher probability of induced microseismicity during and after the end of CO${}_2$ underground injection operations. The theoretical model presented is essential for the geo-mechanical modeling of CO${}_2$ underground injection operations at multiple length-scales.

1. Introduction

Geological sequestration of carbon dioxide (CO${}_2$) into saline aquifers is a broad geophysics application that has fascinated the scientific community for nearly two decades. A specific example is the Illinois Basin Decatur project, which is the first publicly-funded large-scale CO${}_2$ underground storage project in the United States with a total budget of $200 M [1]. There was a surge in microseismic events during CO${}_2$ injection in the Mt. Simon sandstone formation. A significant number of residual microseismic events was also recorded after cessation of CO${}_2$ underground injection activities [2]. To explain these observations, the prevailing hypothesis advanced is that CO${}_2$–rock interactions induce a weakening of the mechanical response with the potential of activating faults [3]. For the Illinois Basin Decatur project, this prevailing hypothesis remains unproven, many years later.

The issue at hand is the prediction of the impact of chemo-mechanical interactions on the nonlinear mechanical response of host rock materials. Many studies have formulated multiscale mechanical models to describe the linear elastic regime. For instance, Bemer et al. [4] developed a composite sphere model to predict the poro-elastic properties of host rocks as a function of porosity and mineralogical content. They validated their theoretical model using tests on limestones, poorly cemented sandstones, and shales. Nguyen et al. [5] formulated a multiscale micromechanical model to represent the influence of double-scale porosity and granular structure on the bulk and shear moduli. Their model was able to reproduce the mechanical weakening observed experimentally following acid-based geochemical alteration. However, the distribution of the porosity, between micro- and macropores, remained unclear. Arson and Vanorio [6] also developed a micromechanical model for the elastic regime that accounted for a double-scale porosity, the pore geometry, and the granular structure of mudstone and packstone. However, their assumption of a monolithic solid phase becomes invalid for highly heterogeneous rocks. We recently formulated a multiscale model to describe changes in elastic constants in Mt. Simon sandstone, while accounting for double-scale porosity, heterogeneity, and granular structure [7]. Our model was validated using nanoindentation tests, helium pycnometry tests, and well wave velocity logs. Herein, we propose to expand our theoretical model to describe the cohesive-frictional response of the Mt. Simon sandstone, when exposed to CO${}_2$/brine mixtures.

The underlying scientific challenge is to model the influence of chemo-mechanical interactions on the strength response. Bemer et al. [4] used empirical relations to model the impact of porosity on the cohesion and friction coefficient of host rocks. Espinoza and coworkers [8,9,10,11,12] formulated a series of numerical models, based on the discrete element method, to represent the influence of granular structure, heterogeneity, and grain size, on the strength response of host rocks. They validated their theoretical models using triaxial tests, scratch tests, and indentation tests. However, their models do not account for the effect of porosity and grain aspect ratio on the mechanical behavior. Moreover, their focus has thus far been on modeling the evolution of strength: they did not investigate the impact of chemo-mechanical interactions on the yield surface and friction coefficient. Thus, there is a need for a physics-based multiscale model to represent the evolution of the yield surface and frictional behavior following rock–CO${}_2$/brine interactions.

The related mechanics-of-material fundamental challenge is to formulate closed-form nonlinear solutions to upscale the strength characteristics of porous granular materials. There is a need for advanced analytical solutions better suited to the complex physics of granular materials with double-scale porosity. Willis and Casta neda built a broad framework to upscale the nonlinear behavior of composite materials [13,14]. However, they considered simple cases for detailed solutions. Dormieux and coworkers [15] formulated a poro-elasto-plastic model for double-scale porous materials with a monolithic solid phase. Pichler and Hellmich [16] developed a multiscale strength upscaling model for granular materials without accounting for dual-scale porosity. To this end, Ulm and coworkers focused on materials with double-scale porosity such as cement [17,18] and shale [19] and formulated cohesive-frictional upscaling models. However, their models rest on a set of restrictive assumptions, including: coupling between the elastic and plastic properties at the smallest lengthscale, and idealized non-porous solid phases at larger scales. Recently, we tackled this difficult problem of nonlinear upscaling for granular composites with dual-scale porosity, focusing on geopolymers [20]. However, geopolymers significantly differ from host-rocks, physics-wise. Thus, we need new fundamental poro-elasto-plastic solutions to predict the mechanical response of host rocks.

We want to correlate induced microseismicity to rock-fluid interactions within the context of CO${}_2$ underground storage. To this end, we formulated a multi-scale poro-mechanics model that connects CO${}_2$-induced alterations to changes in poro-elastic constants, and yield surface characteristics. The outline of this paper is as follows. First, we describe the specimens and experimental protocols. Second, we present our multi-scale micro-mechanics model. Third, we apply our model to understand the impact of rock-fluid interactions on the microstructure, poro-elastic constants, and strength characteristics. Finally, we focus on the Illinois Basin Decatur Project, and we investigate the potential impact of rock exposure to CO${}_2$-saturated brine on the likelihood of induced microseismicity.

2. Materials and Methods

2.1. Materials

Mt. Simon sandstone specimens were extracted from the CO${}_2$-storage zone, via Verification Well 1, in the Illinois Basin—Decatur Project. Figure 1a displays the site location of the Illinois Basin-Decatur Project: It is located in East-Central Illinois and in northern part of the Illinois Basin. Two cores were extracted from Verification Well 1 at two depths: 6925 ft. (2110.7 m), and 6927.55 ft. (2111.5 m). The two cores were sub-sampled in quarter cylinders of a radius of 40 mm, and a height of 80 mm. Figure 1 displays a representative digital photographic image of a core (as extracted).

2.2. CO${}_2$ Alteration Experiments

The core sub-samples were organized into two sets: reference samples and samples for incubation with CO${}_2$/brine. The latter samples were jacketed in a Teflon liner and incubated in CO${}_2$-saturated brine at a temperature of 50 ${}^{\circ }$C and under a pressure of 17.2 MPa for one week. The experimental set-up for sample incubation is shown in Figure 2 [7]. A hollow stainless steel cylindrical tube was utilized for sample incubation. The teflon liner, surrounding the core, served to prevent any interaction between the CO${}_2$-saturated liquid brine and the stainless steel casing. After sealing the incubation cell using Teflon O-rings, four ports on the cell were utilized to inject the brine, to inject supercritical CO${}_2$ using a Teledyne Isco fluid pump, to allow for pressure release, and to monitor the fluid pressure, respectively. Finally, the temperature of the incubation cell was monitored using a set of thermocouples.

In the remainder of the manuscript, we use the notation 6925-U and 6927.55-U, to indicate unaltered Mt. Simon sandstone samples from depth 6925 ft. (2110.7 m) and depth 6927.55 ft. (2111.5 m), respectively. We use the notation 6925-AS1 and 6927.55-AS1 to indicate incubated Mt. Simon sandstone samples extracted from a depth of 6925 ft. (2110.7 m) and 6927.55 ft. (2111.5 m), respectively.

Table 1. Unaltered and altered Mt. Simon sandstone samples that were considered in this study.

Sample Notation	Burial Depth	Alteration Procedure
6925-U	6925 ft.	None
6925-AS1	6925 ft.	Incubation in CO${}_2$-saturated brine at 50 ${}^{\circ }$C and 17.2 MPa for one week
6927.55-U	6927.55 ft.	None
6927.55-AS1	6927.55 ft.	Incubation in CO${}_2$-saturated brine at 50 ${}^{\circ }$C and 17.2 MPa for one week

summarizes the samples used in this study.

Both unaltered and altered specimens were further characterized using optical microscopy imaging, scanning electron microscopy imaging, and grid nanoindentation testing.

2.3. Grinding and Polishing, Optical Microscopy, and Scanning Electron Microscopy Imaging

Before mechanical testing and microstructural imaging, unaltered Mt. Simon sandstone samples and altered Mt. Simon sandstone samples were embedded in epoxy resin and machined into 3-mm thin sections using a low-speed diamond saw. Each thin section was mounted on a 30-mm aluminum disc using cyanoacrylate glue and polished and ground using a semi-automated grinder and polisher. We ground the specimens using silicon carbide pads of different fineness, from an abrasive size of 60 $\mathsf{\mu }$m down to 18 $\mathsf{\mu }$m. Afterward, we polished the specimens using woven polishing pads along with polycrystalline diamond slurries of different grades, ranging from 15 $\mathsf{\mu }$m down to 0.25 $\mathsf{\mu }$m. In between each step, we cleaned the specimens using N-decane in an ultrasonic bath for 2 min.

We assessed the quality of the grinding and polishing procedure using optical microscopy and surface contact profilometry. Figure 3 displays representative optical microscopy images of unaltered Mt. Simon sandstone. Optical microscopy was conducted using a Nikon optical microscope with a 200× optical lens. For all Mt. Simon sandstone specimens, the microstructure is heterogeneous, porous, and granular with grains in the range of 50–500 $\mathsf{\mu }$m. After grinding and polishing, we stored the polished specimens under vacuum for further micro-analysis and nano-mechanical testing.

The micro-structure was further resolved using an environmental scanning electron microscopy FEI Quanta 650 equipped with a backscatter detector under low vacuum at an accelerating voltage of 20–25 kV and a walking distance of 10 mm. Figure 4 shows backscattered environmental scanning electron micrographs of an unaltered Mt. Simon sandstone sample. A multi-scale, heterogeneous, granular, and porous micro-structure is observed. Figure 4a shows grains with a size ranging from 50 to 500 $\mathsf{\mu }$m. Quartz grains are shown in dark gray in Figure 4a–c. Pre-fractured K-feldspar grains are shown in light gray in Figure 4a–c. Figure 4a,b shows microscopic pores with a pore throat size of 20–50 $\mathsf{\mu }$m. Figure 4b,c shows clay minerals surrounding both quartz and feldspar grains. Figure 4d shows micropores with a pore throat size in the range 500–700 nm.

2.4. Grid Nanoindentation Testing

Grid nanoindentation tests were conducted on both unaltered and altered Mt Simon sandstone specimens using an Anton Paar nano-hardness tester with a Berkovich probe. Each indentation grid consisted of an array of 21 × 21 indents with an inter-indent spacing of 50 $\mathsf{\mu }$m. Each indentation test was conducted at a maximum load of 50 mN, a loading/unloading rate of 100 mN/min, and a holding period of 5 s. Figure 5 shows representative load–depth curves for indents on siderite, quartz, and K-feldspar. For each indentation grid, a Gaussian mixture model [7] was applied to assess the distribution of microstructural phases, given the measured values of the indentation modulus and indentation hardness.

The specimens were air-dried before indentation testing. After cleaning in N-decane using an ultrasonic bath, samples were dried using compressed air. Furthermore, prior to testing, the specimens were stored in vacuum desiccators in a dry environment with a controlled relative environment RH = 40%. Afterward, the specimens were tested in a closed environmental chamber with RH = 40%. The resulting penetration depths were small, i.e. less than 600 nm, and probing a volume of characteristic size less than 1.2 $\mathsf{\mu }$m. Moreover, the specimens were stored and tested under ambient conditions. Thus, it is sensible to consider drained conditions to analyze the nanoindentation experiments.

3. Poromechanical Model

3.1. Conceptual Multiscale Model for Host Rock

Figure 6 illustrates our multi-scale conceptual model for the Mt. Simon sandstone from the molecular level to the core level, at the centimeter lengthscale. At the molecular level, or Level 0, our fundamental unit is a potassium feldspar crystal. At the nanometer lengthscale, feldspar is intermixed with micropores. The morphology is statistically disordered. At the microscopic scale, porous feldspar grains are mixed with macropores, quartz, and siderite grains. Finally, at the centimeter lengthscale, the core resembles a homogeneous block.

Our research objective is to connect the mechanical response at Levels 0 (Molecular level), I (Pore level), and II (Core level) to the nano- and microstructure.

Table 2. Mechanical constants and structural parameters at different structural levels.

Level	Mechanical Constants	Structural Parameters
0	$(m_s,{\nu }_s,c_s,{\alpha }_s)$
I	$({\mathbb{C}}^I,{\mathit{b}}^I,N^I,f^I)$	$\phi$
II	$({\mathbb{C}}^{II},{\mathit{b}}_I^{II},{\mathit{b}}_{II}^{II},N_{i,j\in \{1,2\}}^{II},f^{II})$	($\varphi$,${\varphi }_i$, $\delta$)

describes the relevant mechanical constants and structural parameters at each level. At Level 0, the mechanical behavior is characterized via the elasto-plastic constants of the solid skeleton: plane strain modulus $m_s$, Poisson’s ratio ${\nu }_s$, cohesion $c_s$, and internal friction coefficient ${\alpha }_s$. The solid skeleton is assumed to exhibit a linear elastic isotropic behavior and obey a Drucker–Prager yield criterion. At Level I, an important structural parameter is the microporosity $\phi$. The effective stiffness tensor at Level I is denoted ${\mathbb{C}}^I$, and the yield function is denoted $f^I$. At Level II, we consider the macroporosity $\varphi$. We also consider the volume fraction ${\varphi }_i$ of each phase i. The subscript i designates each of the remaining microconstituent, including: microporous feldspar, quartz, and siderite. Finally, we account for the granular aspect ratio, $\delta$. The effective stiffness tensor is ${\mathbb{C}}^{II}$, and the effective yield criterion function is $f^{II}$.

3.2. Upscaling of Poro-Elastic Properties

3.2.1. Molecular Level (0) → Pore Level (I)

At Level I, the poro-elastic stress-strain response is that of a porous material with porosity $\phi$. The poro-elastic constitutive response is given by:

$$\begin{array}{cc}\hfill \mathbf{\sigma }& ={\mathbb{C}}^I:\mathbf{\epsilon }-{\mathit{b}}^{I}p\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \phi -{\phi }^0& ={\mathit{b}}^I:\mathbf{\epsilon }+\frac{p_I}{N^I}\hfill \end{array}$$

(2)

where ${\mathbb{C}}^I$ is the macroscopic stiffness tensor, ${\mathit{b}}^I$ is the tensor of Biot coefficients, and $N^I$ is the skeleton Biot modulus at Level I, that are defined as [21]:

$$\begin{array}{cc}\hfill {\mathbb{C}}^I& =(1-{\phi }^0){\mathbb{C}}_s:{\langle {\mathbb{A}}^I\rangle }_{V_s}\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\mathit{b}}^I& ={\phi }^0\mathbf{1}:{\langle {\mathbb{A}}^I\rangle }_{V_{\phi }}\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \frac{1}{N^I}& =\mathbf{1}:{\mathbb{C}}_s^{-1}:\left(-{\phi }^0\mathbf{1}+{\mathit{b}}^I\right)\hfill \end{array}$$

(5)

where ${\mathbb{C}}_s=3{\kappa }_s\mathbb{J}+2{\mu }_s\mathbb{K}$ is the stiffness tensor of the solid skeleton. ${\kappa }_s=\frac{m_s(1-{\nu }_s^2)}{3(1-2{\nu }_s)}$ is the bulk modulus of the skeleton material. ${\mu }_s=\frac{m_s(1-{\nu }_s)}{2}$ is the shear modulus of the bulk skeleton material. $\mathbf{1}$ is the second-order unit tensor, whereas $\mathbb{J}=\frac{1}{3}\mathbf{1}\otimes \mathbf{1}$ is the spherical projection tensor. $\mathbb{K}=\mathbb{I}-\mathbb{J}$ is the deviatoric projection tensor and $\mathbb{I}$ is the fourth order identity tensor. ${\langle {\mathbb{A}}^I\rangle }_{V_s}$ is the strain concentration tensor of the solid phase and ${\langle {\mathbb{A}}^I\rangle }_{V_{\phi }}$ is the strain concentration tensor of the porous phase. The strain concentration factors can be determined by solving the Eshelby inclusion problem [22]. For instance, if we consider a granular morphology, then the homogenized solid is the reference phase, and the strain concentration tensor of the porous phase is defined by:

$$\begin{array}{cc}\hfill {\langle {\mathbb{A}}^I\rangle }_{V_{\varphi }}& ={\left(\mathbb{I}-\mathbb{S}\right)}^{-1}:{\left((1-{\phi }^0)\mathbb{I}+{\phi }^0{(\mathbb{I}-\mathbb{S})}^{-1}\right)}^{-1}\hfill \end{array}$$

(6)

Here, $\mathbb{S}$ is the Eshelby tensor. For spherical pores within an elastic isotropic matrix of bulk modulus ${\kappa }^0$ and of shear modulus ${\mu }^0$, $\mathbb{S}$ is given by: $\mathbb{S}=\frac{3{\kappa }^0}{3{\kappa }^0+4{\mu }^0}\mathbb{J}+\frac{6({\kappa }^0+2{\mu }^0)}{5(3{\kappa }^0+4{\mu }^0)}\mathbb{K}$ [23]. We solve Equations (3) and (6) using a self-consistent scheme due to the granular and statistically disordered morphology at Level I.

3.2.2. Pore Level (I) → Core Level (II)

At Level II, the material presents itself as a two-scale double-porosity material with a heterogeneous solid phase. Our starting point is Dormieux et al.’s solution for a material with a double-scale porosity and a monolithic solid phase [21]. We then derive new solutions to account for the heterogeneity of the solid phase at Level II. Three phases coexist: macropores with a volume fraction $\varphi$, microporous feldspar with a volume fraction ${\eta }_I$, and quartz with a volume fraction ${\eta }_2$. The fraction of micropores in the microporous feldspar phase is $\phi$. $p_{II}$ is the pressure prevailing in the macroporosity, whereas $p_I$ is the pressure prevailing in the microporosity. The state equations for the effective poro-elastic behavior can be written as [21]:

$$\begin{array}{cc}\hfill \mathsf{\Sigma }& ={\mathbb{C}}^{II}:\mathit{E}-{\mathit{b}}_I^{II}{p}_I-{\mathit{b}}_{II}^{II}{p}_{II}\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \varphi -{\varphi }^0& ={\mathit{b}}_{II}^{II}:\mathit{E}+\frac{p_{II}}{N_{11}^{II}}+\frac{p_I}{N_{12}^{II}}\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \phi -{\phi }^0& ={\mathit{b}}_I^{II}:\mathit{E}+\frac{p_{II}}{N_{21}^{II}}+\frac{p_I}{N_{22}^{II}}\hfill \end{array}$$

(9)

Thus, there are seven poro-elastic parameters: the drained stiffness tensor ${\mathbb{C}}^{II}$; Biot coefficients associated to the microporosity ${\mathit{b}}_I^{II}$ and macroporosity ${\mathit{b}}_{II}^{II}$, respectively; and the skeleton Biot moduli, $N_{11}^{II}$, $N_{12}^{II}$, $N_{21}^{II}$, and $N_{22}^{II}$. Moreover, the expressions for the Biot coefficients and the Biot skeleton moduli are given by:

$$\begin{array}{cc}\hfill {\mathit{b}}_I^{II}& ={\eta }^I{\mathit{b}}^I:\langle {\mathbb{A}\rangle }_I^{II}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\mathit{b}}_{II}^{II}& =\mathbf{1}:\left(\mathbb{I}-{\eta }^I\langle {\mathbb{A}\rangle }_I^{II}-{\eta }_2\langle {\mathbb{A}\rangle }_2^{II}\right)\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \frac{1}{N_{11}^{II}}& =\frac{1}{3{\kappa }^I}\mathbf{1}:\left({\mathit{b}}_{II}^{II}-{\varphi }^0\mathbf{1}\right)\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \frac{1}{N_{12}^{II}}& =\frac{1}{3{\kappa }^I}\mathbf{1}:\left({\mathit{b}}_I^{II}-{\eta }^I{\mathit{b}}^I\right)\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \frac{1}{N_{21}^{II}}& =\frac{1}{3{\kappa }^I}{\mathit{b}}^I:\left({\varphi }^0\mathbf{1}-{\mathit{b}}_{II}^{II}\right)\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \frac{1}{N_{22}^{II}}& =\frac{1}{3{\kappa }^I}{\mathit{b}}^I:\left({\eta }^I{\mathit{b}}^I-{\mathit{b}}_I^{II}\right)+\frac{{\eta }^I}{N^I}\hfill \end{array}$$

(15)

From Equations (10)–(15), the knowledge of the strain concentration tensors and of the macroscopic drained stiffness tensor, at Level II, is enough to derive the remaining macroscopic poro-elastic parameters. We apply Berryman’s model to compute the macroscopic drained bulk and shear moduli. Considering a three-phase model consisting of macropores, feldspar, and quartz, the drained bulk modulus ${\kappa }^{II}$ and drained shear modulus ${\mu }^{II}$ are solutions of the implicit equation [24]:

$$\begin{array}{c}\hfill \sum _i{\eta }_i({\kappa }_i-{\kappa }^{II})P^{*i}=0;\sum _i{\eta }_i({\mu }_i-{\mu }^{II})Q^{*i}=0\end{array}$$

(16)

where ${\kappa }_i$, ${\mu }_i$, and ${\eta }_i$ are the bulk modulus, shear modulus, and volume fraction of phase i, respectively. The functions $(P^{*i},Q^{*i})$ are defined by:

$$\begin{array}{cc}\hfill P^{*i}& =\frac{{\kappa }^{II}+\frac{4}{3}{\mu }_i}{{\kappa }_i+\frac{4}{3}{\mu }_i+\pi \delta {\lambda }^{II}}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill Q^{*i}& =\frac{1}{5}\left(1+\frac{8{\mu }^{II}}{4g_i+\pi \delta ({\mu }^{II}+2{\lambda }^{II})}+2\frac{{\kappa }_i+\frac{2}{3}{\mu }_i+\frac{2}{3}{\mu }^{II}}{{\kappa }_i+\frac{4}{3}{\mu }_i+\pi \delta {\lambda }^{II}}\right)\hfill \end{array}$$

(18)

where ${\lambda }^{II}={\mu }^{II}\frac{3{\kappa }^{II}+{\mu }^{II}}{3{\kappa }^{II}+4{\mu }^{II}}$. Equations (16)–(18) are solved using an iterative numerical scheme. $\delta$ is a parameter that accounts for the aspect ratio of the grains. Finally, we can write the expression of the local strain concentration tensor for each phase (microporous feldspar or quartz): ${\mathbb{A}}_i^{II}=\frac{1}{{\eta }_i}\frac{\partial {\kappa }^{II}}{\partial {\kappa }_i}\mathbb{J}+\frac{1}{{\eta }_i}\frac{\partial {\mu }^{II}}{\partial {\mu }_i}\mathbb{K}$. The implicit equation (Equation (16)) yields:

$$\frac{\partial {\kappa }^{II}}{\partial {\kappa }_i}=1+({\kappa }_i-{\kappa }^{II})\frac{1}{P^{*i}}\frac{\partial P^{*i}}{\partial {\kappa }^i}.$$

3.3. Casta neda and Willis Variational Framework for Strength Upscaling

We apply the linear comparison-composite approach to derive the yield criterion functions. The fundamental idea is to employ a variational approach to obtain the energy function of a nonlinear material [13,25,26]. Specifically, for the strength behavior, we seek to maximize the energy dissipation potential using a comparison-composite medium as a reference [14,27,28]. We use here the dual definition of the plastic domain: using (a) a yield function; and (b) through a plastic potential [29,30]. In what follows, we apply a variational approach to derive the yield criterion for Levels I and II using three steps: (1) derivation of the effective elastic strain rate energy; (2) determination of the local non-linearity functions and local plastic potential; and (3) optimization to yield the global energy dissipation potential.

3.3.1. Molecular Level (0) → Pore Level (I)

We derive new solutions for the homogenized yield criterion of a porous composite. Prior solutions, developed by Ulm and coworkers [18,19] using the linear comparison-composite method, are inapplicable here because the solutions presume a coupling between elasticity and plasticity at the skeleton level along with the presence of rigid inclusions at larger lengthscales. At Level I, we compute the homogenized dissipation potential of a composite consisting of a solid phase and air voids. The solid phase bulk modulus is ${\kappa }_s$, and its shear modulus is ${\mu }_s$. The solid phase obeys a Drucker–Prager criterion with a cohesion $c_s$ and an internal friction coefficient ${\alpha }_s$. The porous phase is characterized by a volume fraction $\phi$ and a packing density $\eta =1-\phi$. The homogenized plastic dissipation potential ${\mathsf{\Pi }}^I$ is the stationarity value of the functional:

$$\begin{array}{c}\hfill {\mathsf{\Pi }}^I=\mathtt{stat}{W}_0^I\left(\mathbf{D}\right)+\eta {\nu }_s\end{array}$$

(19)

$W_0^I$ is the homogenized strain rate energy at Level I. ${\nu }_s$ is the degree of nonlinearity function of the solid phase and is given by:

$$\begin{array}{c}\hfill {\nu }_s=\frac{1}{2{\mu }_s}\frac{{\left(\frac{c_s}{{\alpha }_s}-\tau \right)}^2}{\frac{1}{{\alpha }_0^2}+\frac{1}{{\alpha }_s^2}}\end{array}$$

(20)

with ${\alpha }_0^2=\frac{{\mu }_s}{{\kappa }_s}$. We further define the functionals:

$$\begin{array}{c}\hfill {\kappa }^I=\frac{{\kappa }^I}{{\kappa }_s};{\mathcal{M}}^I=\frac{{\mu }^I}{{\mu }_s}\end{array}$$

(21)

where ${\kappa }^I$ and ${\mu }^I$ are the bulk and shear moduli of the porous composite at Level I, respectively. The homogenized strain rate energy for the porous material reads:

$$\begin{array}{c}\hfill W_0^I\left(\mathbf{D}\right)=\left[\frac{1}{2}\frac{{\kappa }^I}{{\alpha }_0^2}{D}_v^2+{\mathcal{M}}^I{D}_d^2\right]{\mu }_s+{\kappa }^I\tau D_v-\frac{1}{2{\kappa }_s}\left(\eta -{\kappa }^I\right){\tau }^2\end{array}$$

(22)

where $D_v=\mathbf{1}:\mathit{D}$ is the volumetric strain rate and $D_d=\sqrt{\left(\mathit{D}-\frac{1}{3}{D}_v\mathbf{1}\right):\left(\mathit{D}-\frac{1}{3}{D}_v\mathbf{1}\right)}$ is the deviatoric strain rate. The homogenized dissipation potential at Level I is then defined as:

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}^I& ={\mathtt{stat}}_{{\mu }_s,\tau }\left(\left[\frac{1}{2}\frac{{\kappa }^I}{{\alpha }_0^2}{D}_v^2+{\mathcal{M}}^I{D}_d^2\right]{\mu }_s+{\kappa }^I\tau D_v-\frac{1}{2{\kappa }_s}\left(\eta -{\kappa }^I\right){\tau }^2+\frac{\eta }{2{\mu }_s}\frac{{\left(\frac{c_s}{{\alpha }_s}-\tau \right)}^2}{\frac{1}{{\alpha }_0^2}+\frac{1}{{\alpha }_s^2}}\right)\hfill \end{array}$$

(23)

The stationarity equations read:

$$\begin{array}{cc}\hfill \frac{\partial }{\partial {\mu }_s}{\mathsf{\Pi }}^I& =0\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial \tau }{\mathsf{\Pi }}^I& =0\hfill \end{array}$$

(25)

After solving for the stationarity equations, the homogenized dissipation potential reads:

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}^I& =S^I{D}_v+\hfill \\ \hfill & \mathtt{sgn}\left({\chi }^I\right)\sqrt{{\left(A^I\right)}^2{D}_v^2+2\mathtt{sgn}\left({\chi }^I\right){\left(B^I\right)}^2{D}_d^2}\hfill \end{array}$$

(26)

The functionals ${\chi }^I$, $S^I$, $A^I$, and $B^I$ are defined by:

$$\begin{array}{cc}\hfill {\chi }^I& =\eta \frac{{\alpha }_0^2}{{\alpha }_s^2}-{\kappa }^I\left(1+\frac{{\alpha }_0^2}{{\alpha }_s^2}\right)\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill S^I& =-\frac{\eta {\kappa }^I}{{\chi }^I}\frac{c_s}{{\alpha }_s}\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill A^I& =c_s\frac{{\alpha }_0}{{\alpha }_s}\frac{\eta }{|{\chi }^I|}\sqrt{\left(\eta -{\kappa }^I\right)\frac{{\kappa }^I}{{\alpha }_s^2}}\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill B^I& =c_s\frac{{\alpha }_0}{{\alpha }_s}\sqrt{\frac{\eta \left(\eta -{\kappa }^I\right){\mathcal{M}}^I}{|{\chi }^I|}}\hfill \end{array}$$

(30)

We can then derive the homogenized yield criterion at Level I as:

$$\begin{array}{c}\hfill {\left(\frac{{\mathsf{\Sigma }}_m+(p_I-p_{I_0})-S^I}{A^I}\right)}^2+\mathtt{sgn}\left({\chi }^I\right){\left(\frac{{\mathsf{\Sigma }}_d}{\sqrt{2}{B}^I}\right)}^2=1\end{array}$$

(31)

where ${\mathsf{\Sigma }}_m=\frac{1}{3}\mathbf{1}:\mathbf{\sigma }$ is the global volumetric stress, and ${\mathsf{\Sigma }}_d=\sqrt{(\mathsf{\Sigma }-{\mathsf{\Sigma }}_m\mathbf{1}):(\mathsf{\Sigma }-{\mathsf{\Sigma }}_m\mathbf{1})}$ is the deviatoric stress. $\mathtt{sgn}$ is the sign function. $p_I$ is the pressure in micropores and $p_{I_0}$ is the initial pressure in micropores. The macroscopic yield criterion is elliptical for ${\chi }^I>0$ and hyperbolic for $\chi <0$. The resulting degree of non-linearity function reads:

$$\begin{array}{c}\hfill {\nu }^I=\frac{1}{2}\frac{1}{{\mu }^I}\frac{{(S^I-\tau )}^2}{\frac{k^I}{{\mu }^I}-\mathtt{sgn}\left({\chi }^I\right)\frac{{\left(A^I\right)}^2}{{\left(B^I\right)}^2}}+\frac{\mathtt{sgn}\left({\chi }^I\right)}{2}\frac{{\left(B^I\right)}^2}{{\mu }^I}\end{array}$$

(32)

3.3.2. Pore Level (I) → Core Level (II)

At Level II, three phases coexist: macropores, microporous feldspar, and quartz. We model quartz as a cohesive frictionless phase, which is characterized by a volume fraction ${\eta }_2$, a bulk modulus ${\kappa }_2$, a shear modulus ${\mu }_2$, and a cohesion $c_2$. In turn, the volume fraction of macropores is $\varphi$, and the packing density of the microporous feldspar matrix phase is ${\eta }^I=1-\varphi -{\eta }_2$. The generic expression of the macroscopic elastic strain rate energy ${\tilde{W}}_0^{II}$ reads:

$$\begin{array}{c}\hfill {\tilde{W}}_0^{II}=\frac{1}{2}{\kappa }^{II}{D}_v^2+{\mu }^{II}{D}_d^2+{\beta }^{II}{\tau }^I{D}_v-\frac{{\gamma }^{II}}{2}\frac{{{\tau }^I}^2}{{\mu }^I}\end{array}$$

(33)

where ${\tau }^I$ is the prestress in the microporous feldspar matrix and the functions ${\beta }^{II}$ and ${\gamma }^{II}$ are given by:

$$\begin{array}{cc}\hfill {\beta }^{II}& =\frac{{\kappa }^{II}-{\kappa }_2}{{\kappa }^I-{\kappa }_2}\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill {\gamma }^{II}& =\frac{{\mu }^I}{{\kappa }^I-{\kappa }_2}\left({\eta }^I-\frac{{\kappa }^{II}-{\kappa }_2}{{\kappa }^I-{\kappa }_2}\right)\hfill \end{array}$$

(35)

The effective plastic dissipation potential is ${\mathsf{\Pi }}^{II}={\tilde{W}}_0^{II}+{\eta }_1{\nu }^I+{\eta }_2{\nu }_2$. ${\nu }^I$ is the nonlinearity function of the microporous feldspar phase and is given by Equation (32). ${\nu }_2=\frac{c_2^2}{2{\mu }_2}$ is the nonlinearity function of the quartz phase. The effective plastic dissipation potential obeys stationarity conditions with respect to ${\mu }^I$, and with respect to ${\tau }^I$. By enforcing the stationarity conditions, we can show that the generic expression for the effective potential and of the macroscopic failure criterion is:

$$\begin{array}{c}\hfill {\left(\frac{{\mathsf{\Sigma }}_m+(p_{II}-p_{{II}_0})-S^{II}}{A^{II}}\right)}^2+\mathtt{sgn}\left({\chi }^{II}\right){\left(\frac{{\mathsf{\Sigma }}_d}{\sqrt{2}{B}^{II}}\right)}^2=1\end{array}$$

(36)

where the functions $A^{II}$, $S^{II}$, $B^{II}$, and ${\chi }^{II}$ are given by:

$$\begin{array}{cc}\hfill {\chi }^{II}& =\frac{{\kappa }^{\mathcal{II}}}{{{\alpha }_0^I}^2}-\frac{{{\beta }^{II}}^2}{\frac{{\eta }^I}{\zeta }-{\gamma }^{II}}\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill S^{II}& =\frac{\frac{{\eta }^I}{\zeta }{\beta }^{II}}{\frac{{\eta }^I}{\zeta }-{\gamma }^{II}}{S}^I\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill A^{II}& =\sqrt{{\eta }_2\frac{{\mu }^I}{{\mu }_2}{c}_2^2+{\eta }^I\mathtt{sgn}\left({\chi }^I\right){B^I}^2-\frac{\frac{{\eta }^I}{\zeta }{\gamma }^{II}}{\frac{{\eta }^I}{\zeta }-{\gamma }^{II}}{S^I}^2}\sqrt{{\chi }^{II}}\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill B^{II}& =\sqrt{{\eta }_2\frac{{\mu }^I}{{\mu }_2}{c}_2^2+{\eta }^I\mathtt{sgn}\left({\chi }^I\right){B^I}^2-\frac{\frac{{\eta }^I}{\zeta }{\gamma }^{II}}{\frac{{\eta }^I}{\zeta }-{\gamma }^{II}}{S^I}^2}\sqrt{{\mathcal{M}}^{II}}\hfill \end{array}$$

(40)

with $\zeta =\frac{k^I}{{\mu }^I}-\mathtt{sgn}\left({\chi }^I\right)\frac{{\left(A^I\right)}^2}{{\left(B^I\right)}^2}$, ${\kappa }^{II}=\frac{{\kappa }^{II}}{{\kappa }^I}$, ${\mathcal{M}}^{II}=\frac{{\mu }^{II}}{{\mu }^I}$, and ${\alpha }_0^I=\sqrt{\frac{{\mu }^I}{{\kappa }^I}}$. Finally, ${\alpha }^{II}=B^{II}/A^{II}$ is the effective friction coefficient at Level II.

3.3.3. Hardness-Packing Density Relations at Pore Level (I)

At Level I, the material consists of feldspar grains intermixed with micropores. At this level, nanoindentation is employed to characterize the mechanical response. An important task is to correlate the indentation modulus M and indentation hardness H to the microporosity $\phi$ and the elasto-plastic constants of the solid in the feldspar skeleton $(m_s,{\nu }_s,c_s,{\alpha }_s)$. Therefore, we are looking for the functionals $\tilde{\mathcal{M}}$ and $\mathcal{H}$ so that:

$$\begin{array}{c}\hfill M=m_s\tilde{\mathcal{M}}\left({\nu }_s\right);H=c_s\mathcal{H}({\nu }_s,{\alpha }_s)\end{array}$$

(41)

To this end, we carried out 2D finite element simulations of indentation tests in a drained porous material using Abaqus. The initial and the deformed configurations are shown in Figure 7a,b: a Berkovich diamond idealized by a rigid cone of half-apex angle $\theta ={65.27}^{\circ }$ was pushed downward against a cylinder with a prescribed vertical displacement of maximum amplitude 200 nm in a loading/unloading sequence. Symmetric boundary conditions were applied at the left of the material volume, whereas the vertical displacement was restricted at the bottom. Low-friction contact conditions were prescribed between the indenter and the material, with a friction coefficient of 0.01. We conducted a nonlinear geometry analysis to account for potential finite strains.

The mesh consisted of four-node bilinear quadrilateral axisymmetric elements. The mesh was refined right below the probe to capture the gradient of stresses. Meanwhile, a coarse mesh was used away from the probe, to reduce the computational time. The constitutive material behavior consisted of the Level I drained bulk and shear moduli, and of the Level I yield strength criterion (for drained conditions).

We carried out a series of numerical simulations corresponding to values of the packing density $\eta$ between 0.5 and 1, and values of the coefficient of internal friction, ${\alpha }_s$, between 0.2 and 0.7. For each simulation, the computational load–displacement curve was analyzed using the Oliver and Pharr method to extract a theoretical value of the indentation modulus and hardness. Figure 7c,d displays the resulting evolution of the normalized indentation modulus $M/m_s$ and indentation hardness $H/c_s$ as a function of the packing density for ${\nu }_s=0.25$, assuming a self-consistent scheme. As expected, the normalized indentation modulus does not depend on ${\alpha }_s$ and is linear with respect to the packing density $\eta$ with a percolation threshold of 0.5. Meanwhile, the normalized indentation hardness displays a concave evolution with respect to the packing density. $H/cs$ increases as ${\alpha }_s$ increases.

4. Results

4.1. CO${}_2$-Induced Micro-Structural Changes at Pore Level (I) from SEM Imaging

We explored potential microstructural changes, following incubation in CO${}_2$-saturated brine, using scanning electron microscopy and our theoretical model. Figure 8 displays backscattered scanning electron microscopy (SEM) images of Mt. Simon sandstone 6925-AS1, which results from incubation of Mt. Simon sandstone in brine-saturated CO${}_2$ at 17.2 MPa and 50 ${}^{\circ }$C for one week. The SEM images were captured using a FEI Quanta 650 environmental scanning electron microscope under low vacuum with a walking distance of 12.0 mm, an accelerating voltage of 20–25 kV, and magnification levels of 55× to 283×. As shown in Figure 8, a porous granular structure is observed with quartz grains in dark gray and feldspar grains in light gray.

We observe significant microstructural changes by comparing backscattered SEM images of unaltered (6925-U) and altered (6925-AS1) Mt. Simon sandstone specimens. For instance, Figure 4a of unaltered Mt. Simon sandstone 6925-U and Figure 8 displaying altered Mt. Simon sandstone 6925-AS1 were captured with the same equipment, environmental SEM FEI quanta 650, under low vacuum with the same magnification level, 55×, and close values of the walking distance, 12.0–12.6 mm, and accelerating voltage, 20–25 kV. Through a simple comparison, we observed an increase in porosity and a decrease in feldspar grains (in light gray) in the altered sample. Moreover, the increase in porosity following incubation in CO${}_2$-saturated brine is non-uniform, with the formation of channels shown using red arrows (see Figure 8a,b). Finally, there is an increase in the pore throat size from 30–40 $\mathsf{\mu }$m for unaltered Mt. Simon sandstone 6925-U to 70–80 $\mathsf{\mu }$m for altered Mt. Simon sandstone 6925-AS1. Thus, our experimental observations show a change in composition and pore structure following incubation in CO${}_2$-saturated brine.

The observed increase in pore sizes is due to local non-equilibrium flow conditions despite overall static flow incubation conditions without brine recirculation. Consider locally a macropore or micropore filled with acidic brine. Initially, the air void is deformed due to the poro-elastic response. As CO${}_2$-induced dissolution reactions occur, both the air void size and the amount of available reactive surface area increase. This increase in both the air void size and the amount of reactive surface area accessible will lead to an acceleration of the geochemical dissolution reaction while moving the reactive front forward. For more details on pore creation under reactive flow, the reader is referred to [31,32].

4.2. Calibration and Validation of the Theoretical Model

We validate our theoretical model, at Level I, using the nanoindentation tests conducted on the unaltered Mt. Simon sandstone sample, 6925-U. The resulting distribution of indentation modulus and indentation hardness, ${(M_i,H_i)}_{1\le i\le N}$, was analyzed using statistical deconvolution to identify mechanical phases (see Figure A1 and Figure A2). Afterward, the feldspar and microporous phases were isolated. Additional details regarding the statistical deconvolution analysis are provided in Appendix A.

To quantify the local microporosity, nonlinear optimization was performed, using the programming language Python, by minimizing the quantity:

$$\begin{array}{c}\hfill min\sum _i{\left(\frac{H(c_s,m_s,{\nu }_s,{\alpha }_s,{\phi }_i)-H_i}{H_0}\right)}^2+\\ \hfill {\left(\frac{M(c_s,m_s,{\nu }_s,{\alpha }_s,{\phi }_i)-M_i}{M_0}\right)}^2\end{array}$$

(42)

where $M_0$ and $H_0$ are scaling coefficients chosen as the median values of the indentation modulus and indentation hardness, respectively, from indentation tests on macropores and on microporous feldspar.

In Figure 9, we show the calibration of the theoretical model using nanoindentation tests on unaltered Mt. Simon sandstone 6925-U with $N=441$ indentation tests. Figure 9a compares the theoretical and experimental values of the indentation modulus, whereas Figure 9b compares the theoretical predictions with the experimental measurements for the indentation hardness. A good agreement exists between theory and experiments. For the prediction of the indentation modulus, the coefficient of determination is $R^2=0.9260$, and the root mean square error (RMSE) is $\mathtt{RMSE}=6.4815$. Meanwhile, for the prediction of the indentation hardness, the coefficient of determination is $R^2=0.9928$, and the root mean square error is $\mathtt{RMSE}=0.2276$. The optimized elasto-plastic coefficients for the feldspar skeleton solid at the molecular level are: $m_s=70.95$ GPa, ${\nu }_s=0.25$, $c_s=0.37$ GPa, and ${\alpha }_s=0.12$.

Figure 10, Figure 11 and Figure 12 display the validation of our theoretical model for indentation tests carried out on Mt. Simon sandstone samples 6925-AS1, 6927.55-U, and 6927.55-AS1, respectively. For both the indentation modulus and the indentation hardness, a good agreement exists between experiments and theory, as evidenced by the high values of $R^2$ and the low values of the RMSE. Overall, a better agreement exists between theory and experiments for the indentation hardness than for the indentation modulus.

We note that the scatter in the indentation modulus remains small compared to the average value of the indentation modulus. $RMSE/<M>=$7.8%, 8.7%, 8.3%, and 17.00% respectively for Mt. Simon sandstone 6925-U, 6925-AS1, 6927.55-U, and 6927.55-AS1. Thus, the relative RMSE remains low, especially given that for rocks, the uncertainty on the mechanical characterization is usually on the order of 15%. The scatter is due to the strong influence of the structural disorder on the Young’s modulus. Laubie et al. [33] demonstrated that standard micromechanics solutions do not apply when the microstructure significantly deviates from a statistically disordered configuration. In our experiments, however, we note higher values of the RMSE and $RMSE/<M>$ for altered specimens. As shown in Figure 8, CO${}_2$-induced geochemical alterations introduce localized porosity increases, such as channels and wormholes. These alterations lead to a microstructure that is less and less statistically disordered. Nevertheless, our multiscale approach yields a good representation of the mechanical response of the Mt. Simon sandstone. In the next sections, we investigate CO${}_2$-induced changes of the microstructure and macroscopic behavior.

4.3. CO${}_2$-Induced Micro-Structural Changes at Pore Level (I) Using Nanoindentation Testing Integrated with Micromechanical Modeling

The nanoindentation results point toward a change in the pore structure following incubation in CO${}_2$-saturated brine. Figure 13 displays the probability distribution as a function of the local packing density in Mt. Simon sandstone 6925 ft. before (6925-U in Figure 13a and after (6925-AS1 in Figure 13b incubation in CO${}_2$-saturated brine. We calculated the local packing density values using an inverse analysis, as described in Equation (42). We define the macropores by a threshold in the packing density, ${\eta }_{\mathrm{thr}}=0.78$, because macropores are filled with epoxy resin. (The first step in the specimen preparation procedure is to embed samples in epoxy resin under vacuum. The role of the vacuum embedding steps is to protect the specimen from subsequent mechanical damage throughout the grinding and polishing procedure and preserve the microstructure.) We observed a shift of the curve to lower values of the packing density after incubation. We also observed a net increase in macroporosity, as the area under the probability distribution curve between $\eta =0.5$ and $\eta =0.78$ is greater for the altered specimen (see Figure 13b) than for the unaltered specimen (see Figure 13a). Similar observations apply to Mt. Simon sandstone at 6927.55 ft. By comparing the altered state (Figure 14b), and the unaltered state (Figure 14a), we observed a shift of the probability distribution function for the local packing density towards lower values of the local packing density. Furthermore, the area under the probability distribution function of the local packing density, between $\eta =0.5$ and $\eta =0.78$, is greater for the altered specimen than for the unaltered one, indicating an increase in macroporosity.

2

Table 3. Measured alterations in pore structure in Mt. Simon sandstone following incubation in CO${}_2$-saturated brine.

Material	Macroporosity $\mathit{\varphi }$	Microporosity $\mathit{\phi }$
6925-U	0.181	0.083
6925-AS1	0.224	0.096
6927.55-U	0.147	0.105
6927.55-AS1	0.245	0.111

reports the values of the macroporosity and the nano-porosity for Mt. Simon sandstone samples at 6925 and 6927.55 ft., both in the unaltered and altered states. The value of the porosity for unaltered Mt. Simon sandstone 6925-U and 6927.55-U are 18% and 15%, respectively, and are in agreement with well density logs [34].

The nanoindentation tests show an increase in both the macroporosity and the microporosity in host rocks after incubation with CO${}_2$-saturated brine. For the Mt. Simon sandstone samples at 6925 ft., we observe a 24% increase in the macroporosity and a 4% increase in the microporosity. For the Mt. Simon sandstone sample at 6927.55 ft., we observe a 66% increase in the macroporosity and a 6% increase in the microporosity after incubation in CO${}_2$-saturated brine. The observed increase in the macroporosity and the microporosity is higher for the sample at 6927.55 ft. than for the sample at 6925 ft. due to compositional differences. Our results agree with prior studies that reported an increase in porosity in host rocks due to interactions with CO${}_2$ under relevant geological conditions [35,36]. Here, we show that void formation occurs at both the nanometer and the micrometer lengthscales.

4.4. Changes in Cohesive-Frictional Response at the Core Level (II)

Incubation in CO${}_2$-saturated brine results in a change in the yield surface. Figure 15 displays the predicted change in the macroscopic failure envelope under drained conditions for Mt. Simon sandstone samples at 6925 and 6927.55 ft. after incubation with CO${}_2$-saturated brine. For samples at 6925 and 6927.55 ft., in both altered and unaltered states, the predicted failure envelope is an ellipse with an offset. Our theoretical model predicts that incubation in CO${}_2$-saturated brine will shrink the yield surface. This reduction in size of the yield surface, following CO${}_2$-induced alteration, represents a reduction in the magnitude of the statically admissible stresses that the material can sustain.

2 Our theoretical model predicts that, at the macroscopic scale, there will be a mechanical weakening following incubation in CO${}_2$-saturated brine, with a strong impact on the stiffness, the friction coefficient, and the cohesion of the material. The macroscopic constitutive equations were computed by application of the theoretical models derived in Section 3.2 and Section 3.3. The particle aspect ratio $\delta$ was adjusted to yield an agreement between, the theoretical macroscopic Young’s modulus of unaltered Mt. Simon sandstone—6925 ft. and 6927.55 ft.—and shear and compression wave velocity readings [2]. The remaining model parameters (see Table 2) were computed through a statistical deconvolution analysis of grid indentation data, as described in Section 4.2.

Table 4. Model parameters. For each sample—Mt. Simon sandstone 6925 ft. altered and unaltered and Mt. Simon sandstone 6927.55 altered and unaltered—the phase distribution $\left({\varphi }_i\right)$ is given by application of Gaussian mixture modeling to the grid indentation data; meanwhile, the macroporosity $\varphi$ and the microporosity $\phi$ is given by application of our nonlinear micromechanics model. The values of the macroporosity and microporosity are displayed in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13 and summarized in Table 3.

Level	Model Parameters
0	$(m_s=70.95\mathrm{GPa},{\nu }_s=0.25,c_s=0.37\mathrm{GPa},{\alpha }_s=0.12)$
I	$\phi$
II	$(\varphi ,{\varphi }_i,\delta =0.022)$

summarizes the values of the microstructural parameters. We found that the predicted values of the friction coefficient for unaltered Mt. Simon sandstone (see

Table 5. Predicted alterations in drained macroscopic elasto-plastic constants in Mt. Simon sandstone following incubation in CO${}_2$-saturated brine. $S^{II}$ is the offset. ${\alpha }^{II}=B^{II}/A^{II}$ is the macroscopic drained friction coefficient and ${\theta }^{II}=\mathtt{arctan}\left({\alpha }^{II}\sqrt{3}\right)$ is the friction angle. $A^{II}$ (respectively, $B^{II}$) is the long axis (respectively, short axis) of the macroscopic yield curve for drained conditions.

Material	${\mathit{E}}^{\mathit{II}}$ (GPa)	${\mathit{S}}^{\mathit{II}}$ (MPa)	${\mathit{\alpha }}^{\mathit{II}}$	${\mathit{\theta }}^{\mathit{II}}$ (${}^{\circ }$)	${\mathit{B}}^{\mathit{II}}$ (MPa)
6925-U	23.89	−3.45	0.50	40.76	338.57
6925-AS1	20.32	−4.37	0.45	37.71	297.91
6927.55-U	30.34	−4.04	0.54	42.95	379.84
6927.55-AS1	15.93	−10.18	0.39	33.79	229.10

) are congruent with macroscopic triaxial test measurements reported for clay-bearing rocks [37] in general and for Mt. Simon sandstone [2] in particular.

Table 5 report the theoretical predictions for the macroscopic elasto-plastic constants of Mt. Simon sandstone samples at 6925 and 6927.55 ft., before and after incubation in CO${}_2$-saturated brine, as obtained by assuming drained conditions. Figure 15 displays the failure surface, under drained conditions, for both sample depths in unaltered and in an altered state. Incubation with CO${}_2$-saturated brine leads to a reduction of the failure surface for both samples. Following incubation with CO${}_2$-saturated brine for one week at 50 ${}^{\circ }$C and 17.2 MPa, the sample at 6925 ft. experiences a 15% decrease in the macroscopic drained Young’s modulus, an 11% decrease in drained overall friction coefficient, and a 12% decrease in $B^{II}$ (the short axis of the failure envelope).

Thus, our model predicts a higher reduction in strength resistance than in stiffness. A similar trend is observed for the sample at 6927.55 ft. The macroscopic drained Young’s modulus is decreased by 48%, the macroscopic drained friction coefficient is reduced by 28%, and the value of the short axis of the macroscopic drained failure envelope is reduced by 40%, following incubation with CO${}_2$-saturated brine. A higher level of mechanical weakening is predicted for the Mt. Simon sandstone sample at 6927.55 ft. compared to the sample at 6925 ft. This relative increase in mechanical weakening mirrors the relative increase in both the macro- and microporosity following incubation with CO${}_2$-saturated brine. Moreover, our mechanistic model correlates the decrease in the friction coefficient to an increase in macro- and microporosity.

Our predictions conform with the experimental observations of Rathnaweera et al. [38] on quartz-rich sandstone exposed to CO${}_2$-saturated brine at a pressure of 8 MPa and a temperature of 32 ${}^{\circ }$C. Using uniaxial compression tests, they recorded a significant decline in the Young’s modulus and the uniaxial compression strength. Furthermore, SEM and XRD analysis revealed dissolution of quartz and carbonate minerals. Our conclusions also agree with the experimental measurements of Xie et al. [39] on limestone exposed to an acidic solution, resulting in an increase in the porosity. Furthermore, using hydrostatic and triaxial compression tests, they recorded a shrinkage of the yield surface between virgin and altered specimens. Accordingly, our model effectively connects CO${}_2$-induced morphological and compositional changes to alterations in the poro-elasto-plastic response.

5. Discussion

5.1. Origins of Mechanical Strength in Mt. Simon Sandstone

Our model points to the solid skeleton cohesion, the packing density, and the friction between grains as the primary contributors to the strength of Mt. Simon sandstone, in general. An important consequence is that the CO${}_2$-induced geochemical reactions influence rock geomechanical characteristics by altering the packing density locally, and by affecting the coefficient of intergranular friction. More specifically, the increase in micro- and macroporosity is at the root of the mechanical weakening. This increase in air voids, at the micrometer and nanometer levels, results in a decrease in the friction coefficient, a decrease in cohesion values, and a reduction in the yield surface’s size at the macroscopic level.

5.2. CO${}_2$-Induced Geochemical Reactions

The increase in porosity and microporosity, observed for Mt. Simon sandstone, and the predicted mechanical weakening, at the macroscopic level, result from CO${}_2$-induced mineral dissolution reactions. As supercritical CO${}_2$ dissolves into pore waters, the formation of a weak carbonic acid occurs. The increased acidification of the pore waters leads to mineral dissolution and/or precipitation reactions. For the Mt. Simon sandstone, feldspar and quartz are two reactive minerals that are susceptible to undergo dissolution reactions. More specifically, Shi et al. [40] reported illite and feldspar dissolution reactions in Mt. Simon sandstone following incubation with CO${}_2$-saturated brine for 1–2 weeks. Yoksoulian et al. [41] observed clay and quartz dissolution for Mt. Simon sandstone exposed to CO${}_2$-saturated brine under geologic reservoir conditions for six months. Although secondary precipitation has been reported in other host rocks exposed to CO${}_2$-saturated brine under geological conditions, secondary precipitation chemical reactions have not yet been reported for Mt. Simon sandstone. Moreover, our experimental cleaning protocol, which involves a bath in N-decane under ultrasonic vibrations, will most likely remove loosely bonded particles resulting from potential secondary precipitation reactions. Nevertheless, we minimized the impact of the cleaning protocol by using an inert solvent, N-decane, and by using short cleansing times (2 min). Thus, feldspar and quartz dissolution reactions are the primary reason for the observed increase in microporosity and macroporosity for Mt. Simon sandstone samples at 6925 and 6927.55 ft. following incubation with CO${}_2$-saturated brine for one week at 50 ${}^{\circ }$C and 17.2 MPa.

5.3. Implications for CO${}_2$-Induced Microseismicity

The decrease in the friction coefficient and in the cohesion, as observed for Mt. Simon sandstone, following incubation with CO${}_2$-saturated brine, suggests that the likelihood of microseismic events for the host rock will increase after injection of CO${}_2$. The risk of slip, ${\mathcal{P}}_{\mathrm{slip}}$, is described by the applied shear stress, $\tau$, relative to the critical shear stress [42,43]:

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{slip}}=\frac{\tau }{\alpha \left({\sigma }_n-P_f\right)}\end{array}$$

(43)

where $\alpha$ is the host rock friction coefficient, $\tau =\frac{1}{2}\left({\sigma }_1-{\sigma }_3\right)sin2\theta$ is the geological shear stress, ${\sigma }_n=\frac{1}{2}\left({\sigma }_1+{\sigma }_3\right)-\frac{1}{2}\left({\sigma }_1-{\sigma }_3\right)cos2\theta$ is the geological normal stress, $\theta$ is the fault orientation, $P_f$ is the pore fluid pressure, and $({\sigma }_1,{\sigma }_2,{\sigma }_3)$ are the principal stresses of the in-situ geological stress field. For the Mt. Simon sandstone formation, at depths of 6925–6927.55 ft., the in-situ geological stresses are given by [2]: $({\sigma }_1=86.71\mathrm{MPa},{\sigma }_2=58.1\mathrm{MPa},{\sigma }_3=52.03\mathrm{MPa})$. Thus, the threshold pore fluid pressure to induce failure is given by:

$$\begin{array}{c}\hfill P_{f_{\mathrm{thr}}}=min\left({\sigma }_n\pm \frac{\tau }{\alpha }\right)\end{array}$$

(44)

Table 6. Predicted threshold pore fluid pressure $P_{f_{\mathrm{thr}}}$ and critical cohesionless fault orientation for Mt. Simon sandstone in both unaltered and altered states.

Material	${\mathit{P}}_{{\mathit{f}}_{\mathbf{thr}}}$ (MPa)	${\mathit{P}}_{{\mathit{f}}_{\mathbf{thr}}}$ (psi)	Critical Fault Orientation (${}^{\circ }$)
6925-U	26.31	3815	±33
6925-AS1	23.65	3430	± 34
6927.55-U	27.87	4042	± 33
6927.55-AS1	19.39	2812	± 35

reports the predicted values of the threshold pore fluid pressure, $P_{f_{\mathrm{thr}}}$, for Mt. Simon sandstone samples at 6925 and 6927.55 ft. in both altered and unaltered states. Figure 16 displays the evolution of the quantity ${\sigma }_n\pm \frac{\tau }{\alpha }$ as a function of the fault orientation $\theta$. The value of the friction coefficient is given by the macroscopic friction coefficient ${\alpha }^{II}=B^{II}/A^{II}$ as calculated by Equations (39) and (40). Incubation in CO${}_2$-saturated brine mildly affects the critical orientation of the activated fault upon slip. More importantly, incubation in CO${}_2$-saturated brine results in a 10% reduction in the value of the threshold pore-fluid pressure $P_{f_{\mathrm{thr}}}$ for Mt. Simon sandstone at 6925 ft., and in a 30% reduction for Mt. Simon sandstone at 6927.55 ft. This sharp reduction in the threshold pore fluid pressure is mainly driven by the drop in the friction coefficient, which in turn is caused by the CO${}_2$-induced increase in nano- and macroporosity. Thus, the presented workflow enables us to correlate CO${}_2$-induced microstructural changes to a reduction in the threshold pore-fluid pressure.

5.4. Induced Microseismicity Within the Illinois-Basin Decatur Project

Based on our theoretical approach, we carefully examine the case of the US DOE-funded Illinois Basin Decatur Project (IBDP). During three years, from November 2011 to November 2014, 1 million tons of CO${}_2$ were injected into the Mt. Simon sandstone formations at depths of 7025–7050 ft [45]. The frequency of microseismic events was continuously monitored pre- and post CO${}_2$ injection using geophones [46]. A surge in microseismic activity was recorded during the injection phase. Furthermore, after the permanent shut-in in November 2014, a significant number of microseismic events was still recorded, with most post-injection events located closer to the injection well and the projected CO${}_2$ plume [2,31,44].

Figure 16 explores the impact of an increase in the pore-fluid pressure, following CO${}_2$ injection, on the risk of a fault slip. The material in question is the fault gouge that is severely sheared, very porous, and significantly weaker and more compliant relative to the intact host rock [47]. Therefore, we provide a very conservative analysis by using the properties of an intact host rock as the basis for our analysis.

In Zone 3 of the injection well, corresponding to a depth of 6946 ft., the recorded pore-fluid pressures were 3165 psi (21.82 MPa) before CO${}_2$ injection and 3340 psi (23.03 MPa) at the peak of the CO${}_2$ injection [44]. Both values of the pore fluid pressure are well below the threshold pressure required to induce fault slip for any fault located in the Mt. Simon sandstone at depths of 6925 and 6927.55 ft., in the absence of CO${}_2$-induced geo-chemo-mechanical alterations. However, CO${}_2$-induced geo-chemo-mechanical alterations drastically lower the threshold pore-fluid pressure. As a result, for the sample at 6927.55 in the altered state, there is a window of fault orientations for which a fault slip is likely to occur, at pore fluid pressures of both 3340 psi (23.03 MPa) and 3165 psi (21.82 MPa). As for the sample at 6925 ft. in the altered state, it is only during the peak pore fluid pressure of 3340 psi (23.03 MPa) that a window of fault orientations exists that can lead to a fault slip.

Figure 17 follows the behavior of cohesionless faults in four hypothetical sites within the Mt. Simon sandstone location: 6925 ft. in an unaltered state, 6925 ft. in a geochemically altered state, 6927.55 ft. in an unaltered state, and 6927.55 ft. in a geochemically altered state. We focus on three plausible points in time: (1) pre-injection with a reading in October 2011; (2) during CO${}_2$ injection where we select two readings from February 2012 (at the peak level for the pore-fluid pressure) and from January 2013; and (3) post-injection in June 2015, seven months after the end of the CO${}_2$ injection.

Figure 17 reports the risk of slip, ${\mathcal{P}}_{\mathrm{slip}}$, for the different scenarios. Our model predicts that faults in the unaltered regions of the Mt. Simon sandstone (6925 and 6927.55 ft.) would not be activated before, during, and after CO${}_2$ injection: this is expected as the increase in pore-fluid pressure, which is less than 5% of the initial pore-fluid pressure, remains below the threshold pressure of the material. However, due to CO${}_2$-induced microstructural changes, a mechanical weakening occurs. In Figure 17b,c, CO${}_2$ injection is predicted to result in a high risk of fault slip in geochemically altered areas (6925-AS1 and 6927.55-AS1). Our model predicts that after the end of CO${}_2$ injection, some faults that had been activated during injection due to a CO${}_2$-induced mechanical weakening would still remain active (6927.55-AS1). These predictions agree with field observations of post-injection microseismic events. The location of these post-injection microseismic events, near the injection well and close to the projected CO${}_2$ plume, also agrees with our predictions of mechanical weakening due to CO${}_2$ induced microstructural changes. In future work, we will integrate the coupling between geological stress and pore fluid pressure. Nevertheless, the novelty of our work is that our study provides a mechanics-based connection between CO${}_2$-induced microstructural changes and the observed microseismic events during and after CO${}_2$ injection.

6. Conclusions

We investigated the influence of CO${}_2$-induced geochemical alterations on the cohesive frictional behavior of Mt. Simon sandstone within the context of carbon dioxide geological sequestration. To this end, we carried out nanoscale mechanical characterization and microstructural characterization of Mt. Simon sandstone specimens in their native state and following incubation with CO${}_2$-saturated brine for one week under geological conditions. Furthermore, we formulated a multiscale microporomechanical model to correlate the nanostructure, microstructure, and mineralogy of Mt. Simon sandstone to the elasto-plastic parameters at the macroscopic level. Our theoretical model was validated by indentation tests, SEM analysis, sonic velocity logs, and well density logs. Our main findings include:

• Incubation with CO${}_2$-saturated brine, under geological conditions, results in a significant increase in the microporosity and macroporosity.
• CO${}_2$-induced microstructural changes result in a sharp decrease in the size of the yield surface, and in the value of the friction coefficient.
• The reduction in the friction coefficient, following CO${}_2$-induced geochemical alterations, leads to a drop in the threshold pore-fluid pressure and an increase in the risk of fault slip.
• Using the Illinois Basin Decatur Project as a case study, our investigation showed that CO${}_2$-induced geochemical reactions, causing changes in nano- and macroporosity, are a plausible explanation for the surge in microseisimic activity observed during CO${}_2$ injection in the Mt. Simon formation and for the observed residual microseismicity after cessation of the CO${}_2$ injection activities.

In future studies, we will integrate our multiscale nonlinear host rock constitutive model with geomechanical finite element models and reactive transport models [48] to provide a fully coupled analysis of the effect of CO${}_2$-induced geochemical reactions on the risk of induced microseismicity.