Image-based modeling of gas adsorption and deformation in porous media

Understanding adsorption of CO2 in porous formations is crucial to its sequestration in geological formations. We describe a model for adsorption of CO2 and the deformation that it induces in a sandstone formation over wide ranges of temperature and pressure. The model couples the thermodynamics of sorption with elastic deformation of the solid. Finite-element computations are then used in order to compute CO2 adsorption isotherms along with the induced strain in the formation. We also compute the Darcy permeability of the porous medium using the lattice-Boltzmann method. All the computations are carried out with a three-dimensional image of a core sample from Mt. Simon sandstone, the target porous formation for a pilot CO2 sequestration project that is currently being carried out by Illinois State Geological Survey. Thus, no assumptions are made regarding the shape and sizes of the pore throats and pore bodies. The computed CO2 sorption isotherm at 195 K is in excellent agreement with our experimental data. The computed permeability is also in good agreement with the measurement. As a further test we also compute the sorption isotherm of N2 in the same formation at 77.3 K, and show that it is also in good agreement with our experimental data. The model is capable of predicting adsorption of CO2 (or any other gas for that matter) in porous formations at high pressures and temperatures. Thus, it is used to study the effect of hydrostatic pressure on adsorption and deformation of the porous formation under various conditions. We find that the effect of the confining pressure is more prominent at higher temperatures. Also computed is the depth-dependence of the capacity of the formation for CO2 adsorption, along with the induced volumetric strain.

Adsorption-induced swelling of porous formations is a main factor that contributes to the change in the porosity, permeability, diffusivity, and surface area of a pore space 39 .
The main goal of the present paper is to address the aforementioned shortcomings. We describe a statistical mechanical model of gas sorption and the associated deformation in a porous medium. The model is utilized to compute the sorption isotherms and the deformation that it induces in a core sample from a sandstone reservoirthe Mount (Mt.) Simon formation. The formation is generally referred to the basal Cambrian formation in the upper Mississippi Valley and southern Great Lakes areas. Numerical simulation of the theoretical model is carried out using a three-dimensional (3D) image of a core sample from the formation. We then compare the results with our own experimental sorption isotherms and the permeability, measured with the same sample porous medium whose image we utilize in the simulations.

Results and Discussion
Estimating the parameters of the model. As Eqs (5) and (7) in the Methods section indicate, the Hamiltonian (total energy) of the system contains three parameters, namely, the fluid-fluid interaction energy parameter ff  , the fluid-solid interaction energy parameter fs  , and the coupling coefficient λ that links the deformation to sorption. One way of estimating ff  is through the mean-field theory that relates the bulk critical temperature T c to ff  by, where k B is the Boltzmann's constant. As is well-known, however, mean-field theories usually over-or underestimate various properties of fluid-solid systems. A more accurate way of estimating ff  and  fs is assuming that they are the energy interaction parameters that are used in molecular dynamics (MD) simulations. As Table 1 indicates, x-ray diffraction of the core samples of Mt. Simon shows that quartz -silica -constitutes about three-fourth of the sandstone's composition. In fact, the quartz content of several core samples from the same location in Mt. Simon ranges from 73 percent to over 90 percent. Therefore, we use the Lennard-Jones (LJ) energy parameters for CO 2 and CO 2 -silica 40 K. In addition, since we also compute adsorption isotherm of nitrogen in the same core sample, we use the LJ parameters 40 K. In principle, the parameter λ can also be determined by MD simulations 21 . Here, however, we treat it as an adjustable parameter in order to fit the model to the data. Thus, λ, the model's single fitting parameter, was estimated by fitting the calculated sorption isotherms to our experimental data for CO 2 and N 2 adsorption in the core samples at 195 K and 77.29 K, respectively. The optimal values of  λ λ = ⁎ / ff were determined to be −40 and −30 for CO 2 and N 2 , respectively. Negative values of λ * imply that the porous medium swells as a result of gas adsorption.
As described in the Methods section below, we use the mean-field density-functional theory (DFT) to compute the contribution of the gas to the total energy of the system. As in any DFT formulation, one must estimate the maximum adsorption capacity, n max v p , with n max and v p being, respectively, the maximum density of the adsorbed phase, and the sorption pore volume. Their estimates are necessary because the gas density and the adsorbed amounts predicted by the DFT are in units of lattice site occupancy and, thus, must be converted to physical units. The minimum value of n max is equal to reciprocal of the van der Waals co-volume, which for CO 2 is 1.03 g/cm 3 42 . The maximum value of n max is the number density of a close packing of gas particles, assumed to be spherical, which is σ 2 / 3 with σ being the molecular diameter of the gas, which for CO 2 is 1.74 g/cm 3 . Since the maximum theoretical value of n max agrees with the experimental data, we used it in our calculations. The corresponding value for N 2 is 0.723 g/cm 3 . The sorption pore volume v p is that of the core sample, computed based on the experimental values of the density and porosity of the sample whose sorption isotherm we have measured. The pore volume of the sample is 0.094 cm 3 /gr. The maximum adsorption capacity for CO 2 and N 2 are, respectively, 0.164 gr/gr and 0.068 gr/gr. The permeability of the core sample. Measurement of the gas permeability k of the core sample yielded, k = 325 millidarcy (mD). As described in the Modeling section, we computed the permeability of the core sample where M is the total number of grid sites in the pore space, and ρ b is the density of the bulk fluid. The computations were carried out as a function of the pressure P. The bulk density ρ b used in the DFT is in units of lattice site occupancy and is related to the physical bulk density ρ p through, ρ b = ρ p /n max . To estimate ρ p , one may use either the experimental data or an equation of state. Since the Peng-Robinson equation of state 43 provides 44 accurate predictions for the density of gases over a wide range of pressure, we utilized it for calculating ρ p . The bulk chemical potential μ b is given by the standard mean-field equation 45 μ ρ where Z is the coordination number. In the model Z = 4, as we used tetrahedral grid blocks. A series of bulk pressures P were selected and the corresponding bulk density ρ p were calculated using the Peng-Robinson equation of state. Then, the corresponding bulk DFT densities, ρ b = ρ p /n max were computed, and Eq. (3) was used to compute the bulk chemical potential μ b . The results were then used in Eqs (8) and (15) (see the Methods section below) in order to compute the sorption isotherms along with the mechanical response of the porous solid. The bulk pressure does, of course, exert an equivalent external force on the boundaries of the porous sample, which is utilized in the calculation of the volumetric strain via the finite-element formulation that is described below in the Methods section. Figure 1 presents the 3D image of the porous medium and its grid structure used in the computations. The shear modulus G and the Poisson's ratio ν of the solid matrix of the core sample were taken to be 2.73 GPa and 0.26, respectively, which are the typical experimental values reported for sandstone.
The computed sorption isotherms are shown in Fig. 2, where they are compared with our experimental data. It is clear that the data and the computed isotherms agree closely. Also shown are the computed desorption isotherms in the same core sample. The difference between adsorption and desorption isotherms for both gases is small, although it is a bit larger in the case of N 2 . According to the classification of the International Union of Pure and Applied Chemistry 46 (IUPAC), the sorption isotherm of N 2 is of Type III, which corresponds to a case in which attractive adsorbent-adsorbate interactions are weak, but molecular interactions between the adsorbates themselves are strong. On the other hand, although Fig. 2 does not indicate it because the calculations and measurements shown in the figure are restricted to low pressures, in the IUPAC classification the CO 2 sorption isotherm corresponds to Type I, which occurs when adsorption is limited to at most a few molecular layers, so that with increasing pressure the amount adsorbed reaches a constant value and does not change at still higher pressures. This is demonstrated shortly.
The volumetric response of the solid matrix of the porous medium represents its adsorption-induced swelling and elastic compression arising due to the fluid pressure. Interestingly, the volumetric strain "isotherms" in the solid matrix of the core sample, calculated as the sum of the principal strains, corresponds closely to the sorption isotherms of the two gases. This is also shown in Fig. 2, where we present the pressure-dependence of the strain during sorption of the two gases. The deformation of the rock due to adsorption of the gases is small, because the temperature is very low and the core sample is mostly quartz, which is very difficult to deform.
Since the focus of the work is sequestration of CO 2 in rock, hereafter we present the results only for CO 2 . Thus, we simulated sorption of CO 2 in the core sample at relatively high temperatures of 313 K and 350 K and a wide range of pressures that are relevant to the geological conditions for CO 2 sequestration in Mt. Simon. Figure 3 presents the absolute and excess adsorption uptakes of CO 2 at the two temperatures, demonstrating that the isotherms are indeed of the aforementioned Type I. The excess adsorption attains a maximum, followed by a decreasing trend at higher pressures. The reason for the trend is that, at pressures below the maximum the gas densities in the pore space and in the bulk both increase with the same rate. At pressures higher than the pressure at which the maxima are attained, however, the pore space is saturated and, thus, the adsorbed amount no longer change 47,48 . The smaller magnitude of the maximum at the higher temperature and its shift to a higher pressure are, of course, expected. At low pressures the differences between the two isotherms are negligible, consistent with the previous experimental and theoretical studies of CO 2 adsorption in shale and coal formations 22,49 .
The volumetric strain V  and porosity change, induced by CO 2 adsorption in the porous formation at 313 K and 350 K, which correspond to Fig. 3 are depicted in Fig. 4. V  is dependent upon the temperature and is higher at lower temperatures, because the amount of gas adsorbed is higher at low T. The strain "isotherm" follows a pattern similar to sorption, namely, it increases rapidly for pressures up to 10 MPa at 313 K and 15 MPa at 350 K, whereas at still higher pressures the strain changes only slightly, because the pore space has been saturated by CO 2 . Similar trends were reported for CO 2 adsorption in coalbeds 50,51 . The maximum  V due to the adsorption is 0.7% and 0.68% for 313 K and 350 K, respectively, in line with the experimental data for CO 2 adsorption on coalbeds 49 . Note that the pressure at which the strain reaches its maximum corresponds to that at which maximum adsorption is obtained. On the other hand, the porosity increase of the porous formation induced by CO 2 adsorption is only about 1.15% and 1.11% at 313 K and 350 K, respectively. The small change in the porosity is, once again, reflective of the fact that the core sample is largely quartz and very difficult to deform.
Effect of the hydrostatic stress. Since the optimal sites for CO 2 sequestration are located deep within an underground porous formation, it is essential to consider the effect of external confining stress on the adsorption and deformation. Thus, we computed adsorption isotherms and volumetric strains for confining pressure of 50 and 100 MPa. In other word, in addition to the pressure of the gas in the bulk phase, the porous medium experiences an excess stress from the hydrostatic pressure. The results for the hydrostatic pressure of 50 MPa are shown in Fig. 5. The maximum values of the excess adsorption at 313 K and 350 K are, respectively, 1.11 mmol/gr and 0.95 mmol/gr, while the corresponding volumetric strains are 0.69% and 0.67%. These should be compared with the results shown in Fig. 4. The inset of Fig. 5(c) indicates negative strain, or solid contraction, at low pressures, followed by increasing positive strain, indicating swelling. Thus, at low pressures the stress arising from the confining pressure dominates the effect of swelling induced by adsorption. The effect is more prominent at higher temperatures. But, at high pressures, it is adsorption-induced swelling that is more effective than the hydrostatic compression.
The results for a hydrostatic pressure of 100 MPa are shown in Fig. 6. They demonstrate patterns similar to those shown in Fig. 5, with only slight changes in the numerical values of the maximum excess adsorption and the volumetric strain, and with the difference that the contraction-dominated stage occurs over a broader pressure range and higher compressive strains. To better understand the effect of the external pressure, we compare in Fig. 7 the strain-dependence of the absolute and excess amounts of sorbed CO 2 versus volumetric strain at 313 K and 350 K and for three hydrostatic pressures.
Adsorption under geological conditions. Due to the variations of temperature and pressure with depth in geological formations, their effect on the adsorbed amounts of gases cannot be neglected. To this end, we assumed that the pressure and temperature are linear functions of the depth, using a pressure gradient of 15 MPa/km  and temperature gradient of 27.3 °C/km 52 . We then computed the depth-dependence of the maximum adsorption capacity and volumetric strain for CO 2 . The results for sorption are presented in Fig. 8. In shallower depths, up to about 200 m, the difference between the excess and absolute adsorption is negligible. The shapes of the two curves are also similar to what we have presented so far, namely, that the excess amount of adsorbed CO 2 rapidly increases to a maximum value at small depths, followed by a gradual decrease at larger depths. The depth-dependence of the absolute adsorption is also similar to what we have presented so far: It increases sharply at small depths, and then remains essentially unchanged in the samples that are at larger depths. Both sets of results are consistent with the available experimental data for CO 2 sorption in shale and coal 26,52 . The results shown in Fig. 8 and their comparison with the results presented so far suggest that, although high temperatures in deep underground formations adversely affect the maximum adsorption capacity, the effect of high pressures at the same depth is strong enough to offset the unfavorable thermal effect.

Summary
We presented a model for sorption-induced deformation of porous formations that is capable of predicting the sorption isotherms in the formations. The model was used to study CO 2 adsorption in sandstone, along with the associated deformation that it induces in the porous sample. The computations utilized the 3D image of the core sample from Mt. Simon sandstone and, therefore, no assumption was made regarding the shapes and sizes of the pore bodies and pore throats of the pore space. The model's predictions for low-temperature sorption isotherms of CO 2 and N 2 are in good agreement with our experimental data. A study of CO 2 adsorption in the porous medium under a hydrostatic pressure indicated stress-induced compression in the solid matrix at low pressures, followed by swelling at high pressures. T The model is capable of predicting sorption isotherms and the resulting deformation in deep subsurface porous formations, and indicates that in such formations the excess adsorption quickly reaches a maximum, beyond which it gradually decreases, with absolute adsorption and volumetric strain following the same trends with the depth.

Methods
Measurement of the sorption isotherm. Gas adsorption in two cylindrical core samples from Mt. Simon formation were measured. The samples are both horizontal plugs from a larger section that had been extracted along a well in the formation. More specifically, sample 1 was from a depth of 6979.57 ft, whereas sample 2 was from a depth of 6985.85 ft. Their dimensions and weights are listed in Table 2. The mineralogy of the parent core  sample from which the two cores were extracted is presented in Table 1. Prior to the sorption measurements, the samples' porosities were measured via the helium expansion method. They were virtually identical, and are listed in Table 2.
Before sorption experiments began, the samples were dried at 110 °C under vacuum (1 mm Hg) for 24 hours. The N 2 and CO 2 adsorption isotherms were measured using a Micromeritics ASAP 2010 Instrument. The N 2 sorption test was carried out at the nominal temperature of 77.3 K, which was maintained by immersing the specialty-made sample-holder containing the core sample in a liquid-nitrogen bath. Note, however, that the true adsorption temperature varies a bit during the experiment. The instrument, however, measures the true experimental temperature and calculates the corresponding N 2 saturation pressure. The CO 2 adsorption measurements were carried out at 195 K, which was maintained by immersing the sample-holder containing the core in a bath containing a mixture of dry ice and acetone. To generate the adsorption isotherm, the pressure in the sorption cell was increased in a step-wise manner, and the adsorbed gas volume was measured and recorded. The data will be compared with the predictions of the model (see the Results and Discussion section).
The model. Assuming that the sorption-induced swelling is slow and quasi-static, we use equilibrium statistical mechanics to formulate a model that describes sorption and its effect on the deformation and swelling of the pore space. The total energy E of the system -the porous formation and the gas -is given by where E f , E s , and E fs represent, respectively, the energy due to the fluid, the solid, and the interactions between the two. The energy of the fluid is expressed by a mean-field DFT 21,53-55 , where ρ i is the average fluid density (average occupancy) at lattice site i, μ is the chemical potential of the bulk fluid, and ff  is the interaction parameter between neighboring fluid sites. Equation (5) arises from a lattice-gas model in which only the nearest-neighbor interactions are taken into account. The density of each lattice site i is determined by considering two indicator functions I i and  i for each site, such that = I 0 i f ( ) and 1 and = I 1 i s ( ) and 0 denote, respectively, the fluid and the solid indicators. Then, the density at each site is given by, where 〈·〉 indicates the average value. We assume that the solid matrix of the porous formation is linearly elastic (an assumption that can be relaxed, provided that one has a constitutive equation that expresses the elastic properties of the solid 56,57 ). Then, E s , the energy of the solid is represented by  where K i is the elastic stiffness, and σ i and ∇ ⋅ u ( ) i are the stress and strain in element i, respectively. The energy E fs associated with the fluid-solid interactions is the sum of wetting energy E w at the fluid-solid interface, plus the energy change E d in the solid as a result of its sorption-induced deformation/swelling 21,58 : Here, u i is the displacement vector at i,  fs is the wettability coefficient that measures the strength of the interaction between fluid and solid, and λ is the coupling constant that relates the variation in the strain energy of the solid (i.e., its deformation) due to the fluid and solid phases. The sum runs over all the lattice sites. To obtain the equilibrium state of the porous medium and its fluid content, which is in contact with a bulk reservoir at temperature T and chemical potential μ, we minimize the total energy with respect to the fluid density and strain u i of the solid. This completes the formulation of the model.

Numerical simulation of sorption and deformation.
We used the finite-element method to solve the governing equations. Minimizing the total energy with respect to ρ i and u i yields two equations, one for the local fluid density ρ i , and a second one for the displacement field of the solid. The two equations are coupled and must be solved numerically. The governing equation for density of fluid ρ i is given by, In Eq. (9) e ni are the unit vectors from the center of the element i to its vertices n, and k e is the elemental stiffness matrix given by e V T where C and U are, respectively, the elastic coefficient and the strain-displacement matrices, and T denotes the transpose operation. The integration is over the volume V of the element in the computational grid. C is given by, where λ′ and μ′ are the Lamé constants that are related to the shear modulus G and the Poisson's ratio ν, The 3D computational grid that we utilize consists of tetrahedron grid cells; see Fig. 1. We first introduce three normalized coordinates (ξ, η, ζ) ≡ (f 1 , f 2 , f 3 ) that are related to the Cartesian coordinates (x, y, z) by, being the coordinates of the vertices of the tetrahedra. We then introduce a shape vector, N = (N 1 , N 2 , 4 1 , and N j = f j−1 with j = 2 − 4. In Eq. (10) U is a 3 × 4 matrix given by where D is the differentiation operator matrix, so that the entries U ij of U are given by, U ij = ∂N i /∂f j . A global stiffness matrix K is then constructed by compiling in it the stiffness matrices of the tetrahedral elements. Doing so generates the following system of linear equations that describe the displacements of elemental vertices: l nl l i ni i Equations (8) and (15) must then be simultaneously solved numerically. We compute the sorption isotherms of CO 2 and N 2 and the deformation that they induce in the 3D core sample from Mt. Simon sandstone, whose actual 3D image of the sample 59 , shown in Fig. 1(a) simulations. Due to the limitations on the resolution of the image, the porosity of the image is 10 percent. Thus, based on the correlation function of the pore sizes 60,61 we generated pores with sizes below the resolution of the image and distributed them in the image, so that its total porosity matched that of the core sample.
To carry out the computations, we must first discretize the 3D image in order to generate the computational grid. To do so, one must first determine the size of the representative elementary volume (REV) for the sample, i.e., the minimum image size such that its properties will not change if larger images and computational grids are used. In a previous study 59 in which flow of CO 2 and water in the same image was simulated, it was determined that an image with 300 × 300 × 300 voxels represents adequately the REV. Thus, we used the same image size in the present study.
The computational grid must be resolved enough to accurately represent the image's morphology, but it must also be such that the required computations do not take prohibitively long times. We discretized the domain with an adaptive grid with tetrahedral elements, shown in Fig. 1(b). After some preliminary simulations in which various grid resolutions were utilized in order to identify the most accurate grid with affordable computational time, a grid with 373,607 tetrahedral elements was determined to be accurate enough, as the porosity and computed adsorption isotherm of CO 2 did not change significantly if grids with higher resolutions were used.
The numerical procedure for determining the fluid density and displacement fields is as follows. We begin with a given chemical potential μ and initial guesses for the density ρ i (such as a uniform profile across the sample) and ∇ ⋅ u i , and substitute them in the right side of Eq. (8), which provides a new approximation for ρ i . The new approximate ρ i is then substituted in the set of equation (15) and the set is solved by the generalized minimal residual method to obtain a new approximate solution for the displacement field. The new approximations are then substituted in the right side of Eq. (8) to obtain the next approximate solution for ρ i . We then return to Eq. (15) to calculate the new displacement field. The iteration continues until constant density and displacement fields are obtained.
Numerical calculation of the effective permeability. We utilized the lattice-Boltzmann (LB) method with a single relaxation time to simulate fluid flow in the pore space of the image. The effective permeability was computed using Darcy's law. Using the standard bounce-back rule, the no-slip boundary condition was imposed on the internal solid boundaries. A constant pressure gradient was applied in the flow direction. The D3Q19 propagation scheme 62 was used in LB method in order to evolve the flow system.
The LB model is second-order accurate, with its compressibility error related to the square of the Mach number, Ma. In order to simulate an incompressible fluid, Ma is kept below 0.1. A unit relaxation parameter τ R = 1 of the LB method was used, because it leads to accurate simulation of fluid flow 62 . The simulations were carried out in the creeping flow regime with a Reynolds numbers, Re < 1. Some preliminary simulations were carried out to determine an accurate grid resolution and size. We determined that a grid of size 300 × 300 × 300 yields accurate estimate of the permeability with reasonable computation time.