## Abstract

CO_{2} storage in salt rock is simulated with the finite element method (FEM), assuming constant gas pressure. The initial state is determined by simulating cavity excavation with a continuum damage mechanics (CDM) model. A micro–macro healing mechanics model is proposed to understand the time-dependent behaviour of halite during the storage phase. Salt is viewed as an assembly of porous spherical inclusions that contain three orthogonal planes of discontinuity. Eshelby's self-consistent theory is employed to homogenize the distribution of stresses and strains of the inclusions at the scale of a representative elementary volume (REV). Pressure solution results in inclusion deformation, considered as eigenstrain, and in inclusion stiffness changes. The micro–macro healing model is calibrated against Spiers’ oedometer test results, with uniformly distributed contact plane orientations. FEM simulations show that independent of salt diffusion properties, healing is limited by stress redistributions that occur around the cavity during pressure solution. In standard geological storage conditions, the displacements at the cavity wall occur within the first 5 days of storage and the damage is reduced by only 2%. These conclusions still need to be confirmed by simulations that account for changes in gas temperature and pressure over time. For now, the proposed modelling framework can be applied to optimize crushed salt back-filling materials and can be extended to other self-healing materials.

The low permeability and self-healing potential of salt rock (halite) makes it a good potential host material for long-term geological storage, especially in wet conditions. So far, the actual use of salt formations for geological storage has been limited by technological challenges related to retrievability. However, in favourable stress, pore pressure and temperature conditions, salt stiffness and strength can indeed increase due to pressure solution, solid diffusion and recrystallization (Chan *et al.* 1998; Tsang *et al.* 2005; Zhu & Arson 2015). Pressure solution is a very effective healing mechanism, common in crystalline media that contain water films, especially in halite. Salt minerals dissolve at contacts that are under high stress, diffuse along grain boundaries and reprecipitate at pore walls, which are under lower stress (Paterson 1973; Raj 1982; Rutter 1983). Based on thermodynamic equations established at the grain scale, experiments were conducted on granular salt and phenomenological models were proposed to predict healing in halite (Spiers *et al.* 1990; Yang *et al.* 1999; Houben *et al.* 2013). In this paper, we use a homogenization approach to understand the macroscopic effect of local pressure-solution mechanisms at the scale of a representative elementary volume (REV).

Eshelby's theory allows the calculation of the stress and strain fields of a REV made of an elastic matrix that contains ellipsoidal heterogeneities (Eshelby 1957). Based on this theory, Mori & Tanaka (1973) proposed an explicit method to predict the homogenized stiffness of the REV, in which the interaction between the matrix and the inclusions is accounted for. However, when it is impossible to identify a dominating phase that can be considered as the matrix, as in polycrystalline materials, it is necessary to use a so-called self-consistent method (Kröner 1961; Hill 1965), in which the REV is seen as a juxtaposition of inclusions. Each inhomogeneity is seen as an inclusion embedded in a matrix that has the yet-unknown homogenized properties of the REV. Hence, the properties of the matrix are not known *a priori*, which makes the model implicit. In what follows, we use a self-consistent method to upscale stresses and strains induced by local pressure solution. Since the time-dependent strains of the inclusions depend on thermodynamic processes that cannot be predicted from the far-field stresses alone, we model local pressure-solution strains as eigenstrains (e.g. Pichler & Hellmich 2010). We then use our micro–macro model to simulate healing around salt cavities used for CO_{2} storage.

We first explain the pressure-solution phenomenon and present the corresponding thermodynamic equations at the inclusion scale. Secondly, we formulate the micro–macro healing model based on a self-consistent homogenization scheme. We then calibrate the model against published experimental data. Lastly, we present finite element simulations of salt-cavity healing during CO_{2} storage.

## Pressure-solution model

Let us consider two halite crystals separated by a thin fluid film. A stress increase normal to the crystal boundary (called the contact plane) leads to an increase in the chemical potential (Δ*μ*) in the solid constituent in reference to the solute:*σ*_{e} is the effective stress normal to the contact plane (defined as the difference between the normal stress and the fluid pressure) and *Ω* is the molar volume of NaCl (2.7 × 10^{−5} m^{3} mol^{−1}). The drop in the chemical potential in the solute can be expressed as a function of mineral concentration:*R** is the gas constant, *T* is the Kelvin temperature, and Δ*C* is the difference between the ionic concentration in the fluid films located at contact planes and that in the pores. *C*_{o} is the reference concentration in the pores, located at crystals’ edges (considered as sinks). Due to the difference of solute concentration between the contact planes and the pores, salt ions diffuse from high-concentration sites (contact planes) to low-concentration sites (pores). Figure 1 shows a schematic representation of the diffusion path around a pore. We define the elementary heterogeneity of salt rock as a hollow spherical inclusion that contains a pore located at the intersection of three orthogonal crystal contact planes (plane *XY*, plane *YZ* and plane *XZ* in Fig. 1). The salt crystals are assumed to be cubic with dimensions of 2*r*_{g} (Overton & Swim 1951; Zhang & Spiers 2005). The diameters of the grains of halite range between 0.01 mm and several decimetres (Urai & Spiers 2007; Urai *et al.* 2008). In this paper, we work with a grain dimension *r*_{g} of the order of 10^{−4} m, which is the typical size of salt crystals studied in the laboratory. This will make it possible to calibrate our model against reported experimental data in the sequel. The thickness of the solid wall around the pore is assumed to be uniform, so that the radius of an inclusion is equal to *r*_{g}. Noting *W*, the thickness of the shell around the pore, the radius of the pore is *r*_{g}−*W*.

The pressure-solution mechanism is illustrated in Figure 2. Salt dissolves at contact planes that are under high stress. We consider that pores act as sinks, so that ions diffuse towards the central pore and uniformly precipitate on the pore wall. The increase in the pore-wall thickness is noted as δ*W*.

According to Fick's first law, the radial diffusion flux, *J*(*w*), along the contact plane is related to the diffusion coefficient, *D*, and the mineral concentration, *C*, as follows:*w* is the distance between the periphery of the contact surface and the dissolution front (i.e. the wall thickness at a given time). The shape of the pore–film interface is an annulus, the area of which can be calculated as the product of the film thickness, *S*, by the perimeter of the pore wall, *V*_{c} is the dissolution velocity. Pressure solution induces energy dissipation by diffusion. The energy dissipation per unit volume (*C*(*w*), is equal to the solute concentration in the pores (*C*_{o}) (Rutter 1983; Schutjens & Spiers 1999; Pluymakers & Spiers 2015). The total energy dissipated along the contact plane *XY* is obtained as:*V _{XY}* is the dissolution velocity on the plane

*XY*. We assume that the total input work at a contact plane is entirely dissipated by diffusion. The velocity

*V*can thus be calculated as:

_{XY}*XY*. With the expression of the dissolution velocity, we obtain the inclusion strain rate and the change in pore-wall thickness, as follows:

*A*

_{s}is the surface area of the pore in the centre of the inclusion, and

*γ*is the chemical strain.

_{z}## Homogenization scheme

In order to predict the stiffness, deformation and stress of salt rock at the REV scale, a self-consistent homogenization scheme is adopted. The strain of a spherical inclusion can be expressed as (Dvorak 1992):*ε ^{i}* and

*A*are the strain and concentration tensors of inclusion

^{i}*i*, respectively;

*n*is the total number of inclusions in the REV (characterized by their porosity and plane orientations);

*D*is the influence tensor of inclusion

^{ij}*j*on inclusion

*i*;

*γ*is the eigenstrain of inclusion

^{j}*j*(the chemical strain induced by pressure solution); and

*E*is the strain of the REV. Based on Eshelby's inclusion-matrix theory, the concentration tensor

*A*can be calculated as (Dormieux

^{i}*et al*. 2006):

*φ*

^{j}is the volume fraction of inclusion

*j*in the REV,

*P*

^{i}is the P tensor (also called Hill's tensor in some references),

*C*is the stiffness of inclusion

^{i}*i*and

*C*is the homogenized stiffness of the REV.

^{h}*P*

^{i}is a fourth-order tensor that depends on the shape of the inclusion and on the stiffness of the REV. The explicit expression of the P tensor was given by Mura (1987).

*C*is calculated based on a damage model, to be determined. Two examples of such damage models will be given, one at laboratory scale for calibration and verification purposes, and one at cavity scale for numerical simulations. The homogenized stiffness

^{i}*C*can be calculated implicitly, as follows:

^{h}*D*are expressed as follows (Pichler & Hellmich 2010):

^{ij}*t*, the strain of the REV (

*E*), the chemical strain (

_{t}## Calibration against oedometer test results

We first applied our model at laboratory scale. We closed the formulation by introducing a damage model at inclusion scale. In granular materials, the shear modulus depends on porosity (Kováčik 2008) and increases with the hydrostatic stress (Digby 1981). Accordingly, for each inclusion, we assumed that the shear modulus, *μ**, was related to the porosity of the inclusion, *ϕ*, and to the effective stress, *σ*_{e}, as follows:*η*, *m* and *n* are constants that need to be calibrated. Note that *η*, *m* and *n* are assumed to be the same for all inclusions. *μ** and *ϕ* are inclusion-specific (but the *i* index was dropped in equation 16 for clarity). The stiffness of inclusion *i* is expressed as:** I** is the second-rank identity tensor and

*δ*is the symmetrical part of the fourth-rank identity tensor.

We simulated oedometer tests conducted by Spiers’ group on saturated granular salt (Spiers & Brzesowsky 1993). The micro–macro model presented above was implemented in a finite-element method (FEM) package and simulations were done with only one element, to reflect the behaviour at the material point. In conformity with the experiments, the crystal and void size distributions were assumed uniform in the REV, with an initial porosity of 42% and a crystal size of *r*_{g} = 0.1375 mm (equal to the inclusion size in our microstructure model). We assumed that the distribution of contact plane orientations was uniform. Numerical creep curves were fitted to the experimental results obtained under compressions of −8.0, −4.2 and −1.1 MPa (Fig. 3). The corresponding calibrated parameters are reported in Table 1. We define the prediction error as the difference between the area below the experimental creep curve and the area below the numerical creep curve, normalized by the area below the experimental creep curve. The errors made in the numerical predictions are 5.7, 4.0 and 3.7% for oedometer tests under −8.0, −4.2 and −1.1 MPa, respectively.

During the oedometer tests, the creep rate decreases due to: (i) the increase of the diffusion path length in each inclusion: precipitation at the pore wall increases the migration distance of the ions from the periphery of the inclusion to the pore wall; and (ii) the decrease in compressive stress at inclusion and REV scales: during pressure solution, the dimensions of the inclusions decrease in the directions normal to the contact planes where dissolution occurs; according to equations (7) and (8), compressive stress at contact planes is relaxed, which leads to smaller chemical potential differences between the solid and the solute, and hence less dissolution.

## Simulation of CO_{2} storage in a salt cavern

Secondly, we applied our model at cavity scale. We closed the formulation by introducing a damage model at REV scale. The FEM is used to simulate CO_{2} storage in a salt cavern. Based on the work of Dusseault *et al*. (2004), the salt cavern was modelled as an oblate spheroid with a vertical axis of 100 m and horizontal axes of 150 m. The centroid of the cavern was located at a depth of 1200 m and halite density was taken as equal to 2400 kg m^{−3}. After excavation, the cavern was sealed at a pressure of 15 MPa. The FEM model's dimensions were 1650 × 375 × 375 m (Fig. 4). The lateral displacement of salt rock at the external boundary was fixed. We used a continuum damage mechanics (CDM) model (described below) to simulate the excavation phase. In a second stage, the micro–macro healing model presented above was used to simulate the storage phase, by applying a 15 MPa pressure at the cavern wall. Note that the incremental relationship is obtained by a local homogenization and upscaling method, and not through a macroscopic law previously established. The initial value of porosity in the healing simulation phase was set equal to the value of the damage variable, *Ω**, calculated after the excavation phase. During the storage phase, porosity decreases due to pressure solution, which is expected to increase halite stiffness – a process referred to as ‘self-healing’. We first performed a mesh sensitivity analysis. Results confirmed that mesh size only influenced the damage value after the excavation phase, but has very little effects on the predicting of healing.

A thermodynamics-based CDM model is used to predict excavation damage. At REV scale, the expression of Helmholtz free energy, Ψ_{s}, is given as (Halm & Dragon 1998; Zhu & Arson 2015):*λ*_{o} and *μ*_{o} are Lamé constants; *α* and *β* are damage material parameters; and *Ω** is the dimensionless damage variable, considered as a scalar porosity in the present study. The damage driving force *Y*_{d} is expressed from thermodynamic conjugation relationships, as follows:*k*_{o} is the damage initiation threshold and *k*_{1} is a damage hardening parameter. During the initiation and propagation of damage, the consistency conditions hold (i.e. *f*_{d} = 0 and δ*f*_{d} = 0). The damage parameters used for the simulation, taken from a prior study (Zhu & Arson 2015), are reported in Table 2.

Figure 5 shows the distribution of damage around the salt cavern after excavation. The maximum damage observed reaches 9.8% and appears in the middle of the sidewall, due to the high compressive principal stress. Damage drops rapidly away from the cavern wall. The size of the damage zone is of the order of the cavern's radius.

In a second step, we use the micro–macro healing model with randomly orientated inclusions to simulate the storage phase. The REV porosity is calculated from the porosity of the inclusions, and used as the damage variable in the expression of REV stiffness that derives from the expression of the free energy in the CDM model above (equation 17). The minimum (compressive) principal stress distribution around the cavern is represented in Figure 6, for *DS* = 1.0 × 10^{−18} m^{3} and *r*_{g} = 0.13 mm. The corresponding evolution of damage at the cavern wall is shown in Figure 7. The evolution of cavity convergence is illustrated in Figure 8.

Figure 6 clearly shows that the compressive principal stress decreases over time during the storage phase. This is because the dissolution of salt at the contact planes triggers negative chemical strains at the inclusion scale. Since halite elements are geometrically constrained, inclusion shrinkage results in tensile stresses at the REV scale, which reduce the overall compressive stresses around the cavern. The decrease in compressive stress is particularly visible in the elements with large initial damage in the middle of the sidewall, where the compressive principal stress switches from −43 MPa after excavation to −15 MPa during storage. Figure 7 confirms that stress redistribution at the cavern wall is due to healing: that is, to the reduction of damage, defined here as porosity. Because pressure solution relaxes compressive stress at the contact planes, the dissolution rate decreases over time as salt precipitation occurs at the pore walls. In other words, healing is self-limited. As a result, the healing rate decreases rapidly during the first 3 days of storage, after which damage reaches a plateau. At the sidewall, damage amounts to 9.8% after excavation and to 9.6% after 4 days of storage. Healing thus reduced the maximum damage by 2%. Just after excavation, the vertical convergence is −0.048 m, and the horizontal convergence is +0.015 m (see Fig. 8). Pressure solution increases convergence, which reaches 1% in both the horizontal and vertical directions after 4 days of storage. Cavern deformation stabilizes when healing reaches a plateau.

Figure 9 presents a sensitivity analysis of the parameter *DS*, which is the product of the diffusion coefficient (*D*) by the thickness of the fluid films at contact planes (*S*). Clearly, the healing rate increases with *DS*, which can be seen as a diffusive efficiency parameter. Physically, pressure solution occurs faster when the diffusion coefficient is higher and/or when the intercrystalline space is larger. However, *DS* does not influence the final damage value at the cavity wall, because the rate of pressure solution is independent of the stress redistributions that occur around the cavern (Fig. 9a). Similarly, a high diffusive efficiency accelerates convergence but does not influence the final shape of the storage facility (Fig. 9b).

Figure 10 illustrates the influence of the size of the inclusions on salt healing. Inclusion size is equivalent to crystal size and indicates how long the diffusion path is from high stress dissolution sites to low stress precipitation sites. Larger inclusions imply longer diffusion paths and thus a lower healing rate. At the same initial porosity, larger inclusions also mean fewer inclusions. After 2 days of storage, healing with a crystal size of 0.26 mm is only 50% of the healing obtained with a crystal size of 0.13 mm, and the asymptotic values of damage and convergence are reached after 20 days. The healing rate is negligible with a crystal size of 1.30 mm, which shows that the healing rate is more sensitive to the crystal size than to the diffusive efficiency.

## Conclusions

In this paper, we propose a micro–macro healing mechanics model that captures the effect of pressure solution on halite porosity. Thermodynamic equations of pressure solution are established at the scale of an inclusion defined as a hollow sphere intersected by three orthogonal contact planes that contain a fluid film. Stress at the contact planes determines the dissolution rate. The pore at the centre of the inclusion is viewed as a sink: ions diffuse from the dissolution sites to the pore wall, where they precipitate. The resulting rate of deformation of the inclusion is calculated from mass-balance equations, and then defined as an eigenstrain in a self-consistent homogenization scheme. The micro–macro healing model is calibrated against published oedometer test results and implemented in a FEM package.

CO_{2} storage in a deep oblate spheroidal salt cavity is simulated. A continuum damage model is used to calculate the damage induced by excavation. The micro–macro healing model is used to simulate the storage phase under constant gas pressure. In real storage conditions, gas pressure increases due to convergence. However, in the present study, convergence was limited to an asymptotic value of 1%, which was not found to influence gas pressure. Due to the healing process, stress around the salt cavern is redistributed and becomes uniform. The largest compressive principal stress is decreased by 65%. Under these conditions, numerical calculations show that only 2% of the excavation damage can be recovered at the sidewall, where damage is the highest. Damage and convergence reach a plateau after 4 days of storage. The healing rate decreases over time because: (i) salt precipitation at the pore walls lengthens the diffusion path of the ions dissolved at high compression stress sites; and (ii) healing results in a stress redistribution around the cavity, reducing compressive stresses, thus limiting the triggering of pressure solution. Higher diffusion coefficients and thicker fluid films can accelerate healing, but cannot change its asymptotic value. A larger crystal size significantly decreases the healing rate.

Simulation results indicate limited healing potential around salt cavities used for CO_{2} storage, but it has to be noted that the healing model was calibrated at laboratory scale for a uniform size distribution of crystals, which is not representative of realistic storage conditions. In addition, temperature effects were omitted in this study. In future work, we will study the effect of grain size distribution, and formulate a more comprehensive damage and healing thermohydromechanical model. As is, the proposed modelling framework can be used to optimize some microstructure parameters of crushed salt buffers (such as porosity and fluid inclusion distribution) and be extended to other self-healing materials.

## Funding

This work was supported by the United States National Science Foundation, under Grant No. CMMI#1552368: ‘CAREER: Multiphysics Damage and Healing of Rocks for Performance Enhancement of Geo-Storage Systems - A Bottom-Up Research and Education Approach’.

- © 2019 The Author(s). Published by The Geological Society of London for GSL and EAGE. All rights reserved