Porosity Evolution during Chemo-Mechanical Compaction

Abstract

This chapter presents the constitutive equations necessary to interpret laboratory and field data when both solid and pore volume evolve through time due to chemical and mechanical processes. The equations for the porosity evolution that are developed are generic, but the examples presented are acquired from chalk core studies. The processes at play when porosity is subject to change due to volumetric compaction and fluid-rock interactions when porous chalks are continuously flooded are presented here. As the overall solid mass is a conserved quantity, the void space is not. Constitutive equations are therefore required to estimate the time-evolution of the porosity. Laboratory triaxial tests were performed on high-porosity outcrop chalks from Obourg, Liegè, and Mons (Belgium). These tests are being compacted and continuously flooded with MgCl2 brine at elevated temperature and at high stresses. As calcite is replaced by magnesite, the overall mass and solid density change, thereby changing the volume of the solid. At the same time, the bulk volume is changing. Taking both effects into consideration, the pore volume evolution can be determined. We find that the porosity changes in nonintuitive ways as the relative importance of bulk compaction and chemical interaction may vary over time.

1. Introduction

Understanding how chemically reactive and mechanically deformable rock-fluid systems consisting of solids and voids evolve is very important to several fields in the Earth sciences. Examples include the lithification of sedimentary strata [1] and long-term creep behavior of crustal rocks [2]. In addition, a range of industrial processes are affected by chemo-mechanical interactions, including, e.g., pharmaceutical and food processing industries, and geotechnical engineering concerning roadwork construction, mass transportation, and slope stability.

The pore volume fraction, the pore size distribution, and the mineral surfaces are key parameters to ensure safe disposal of radioactive waste and captured CO₂, and to understand how ores’ deposit evolves with time. In petroleum sciences, chemo-mechanical processes are important to accurately predict the porosity since it is inside the pores where hydrocarbons are stored, and it is through the pores, the hydrocarbons are being produced by miscible and immiscible fluid migration across reactive mineral surfaces that again are subject to change. Both pore volume and production rate are crucial to determine the recoverable hydrocarbon potential.

Reactive pore fluids in nonequilibrium with their host rocks lead to dissolution and precipitation transforming the mineral assembly into another, see for example, [3, 4, 5]. Dissolution and precipitation lead to changes to grain texture and morphology, and the mineral surface’s affinity to oil and water. These are factors, together with porosity, that dictate not only the flow property of the porous rock but also the mechanical parameters that control the stiffness, strength, and the rate at which compaction by grain reorganization and pore collapse occurs. The general processes that are described here is applied to understand how the porosity of chalks develops (Figure 1a, b display unaltered and altered chalks) dynamically in a controlled triaxial cell experiments (Figure 2a), with control of temperature, pore pressure, side stress, and overburden stress of cylindrical samples (Figure 2b).

Figure 1.

(a) SEM image of an unaltered chalk (Liegè, Belgium [3]). Calcite grains partially organized in coccolith rings and foraminifers. (b) Reworked Liegè chalk from the same core as (a) after 1090 days of continuous mechanical compaction and flow of reactive 0.219 M MgCl₂ brine at 130°C (Table 1).

Figure 2.

(a) Triaxial cell setup controlling axial and radial stress, the pore pressure, flow rate, and temperature. (b) Additive partitioning of the total bulk strain into a pore and solid volume component. Here, uniaxial strain is assumed (constant diameter) such that lengths relate to volumes.

This chapter deals with some of the constitutive relations that are used to describe the evolution of porous bodies. We incorporate a discussion of how rock-fluid chemistry may impact the grain volume, and review the ways in which total volume reduction may be facilitated in compressive systems. The discussion summarized in the development of a porosity evolution equation in which all effects are included. The usage of the porosity evolution equation is exemplified with references to already published experimental results.

2. Constitutive equations for porosity evolution

The basic equations that are used to quantify the porosity evolution through time are presented. The analysis is based on the work presented in Nermoen, et al. [3]. The overall bulk volume of a bi-phase material equals the sum of the solid volume and pore volume

$${\mathrm{V}}_{\mathrm{b}}={\mathrm{V}}_{\mathrm{s}}+{\mathrm{V}}_{\mathrm{p}}$$E1

Any changes in solid volume and pore volume lead to changes in the bulk

$${\mathrm{\Delta }V}_b={\mathrm{\Delta }V}_s+{\mathrm{\Delta }V}_p$$E2

The pore volume, and hence the porosity, itself is not a conserved quantity. In that case, the bulk volume (size of the object of study) and the solid volume evolution have to be used. Since the volumes are additive by nature, the changes in pore volume can be calculated

$${\mathrm{\Delta }V}_p={\mathrm{\Delta }V}_b+{\mathrm{\Delta }V}_s$$E3

At any given time through dynamic porosity evolution, the porosity is given by

$$\varphi =\frac{V_p}{V_b}=1-\frac{V_s}{V_b}$$E4

When both the bulk volume and the pore volume change dynamically from known measurements before the experiment starts (${\mathrm{V}}_{\mathrm{b},0}$ and ${\mathrm{V}}_{\mathrm{p},0}$ are known), then the time-evolution of the porosity is given by

$$\varphi \left(t\right)=\frac{V_{p,0}+\mathrm{\Delta }{V}_p\left(t\right)}{V_{b,0}+{\mathrm{\Delta }V}_b\left(t\right)}$$E5

Using Eq. (3) enables the determination of the porosity from known quantities

$$\varphi \left(t\right)=\frac{V_{p,0}+\mathrm{\Delta }{V}_b\left(t\right)-\mathrm{\Delta }{V}_s\left(t\right)}{V_{b,0}+{\mathrm{\Delta }V}_b\left(t\right)}$$E6

This equation is useful when determining pore volume evolution when considering mechanical and chemical processes that occur at reactive rock-fluid systems exposed to elevated stresses. To simplify the porosity evolution equation further, the volumetric strain and the initial porosity before chemo-mechanical processes occur are introduced

$${\epsilon }_{\mathit{vol}}\left(t\right)=-\frac{\mathrm{\Delta }{V}_b\left(t\right)}{V_{b,0}}\text{and}{\varphi }_0=\frac{V_{p,0}}{V_{b,0}}$$E7

The minus sign in the volumetric strain here are in line with the definition in geotechnical engineering that inward deformation is positive, often different from other fields of sciences. Dividing by the initial bulk volume and employing the definitions Eq. (6) become

$$\varphi \left(t\right)=\frac{{\varphi }_0-{\mathrm{\epsilon }}_{\mathit{vol}}-\mathrm{\Delta }{V}_s\left(t\right)/V_{b,0}}{1-{\epsilon }_{\mathit{vol}}}$$E8

Eq. 8 is used to analyze how the pore volume fraction changes as the overall volume and the solid volume changes through time. Typically, it is easier to quantify the changes in the solid volume and total volume because of conservation of mass, but this does not generally apply. In other cases, when the pore volume and solid volume are known, the porosity can be calculated from

$$\varphi \left(t\right)=\frac{V_{p,0}+\mathrm{\Delta }{V}_p\left(t\right)}{V_{p,0}+V_{s,0}+{\mathrm{\Delta }V}_p\left(t\right)+{\mathrm{\Delta }V}_s\left(t\right)}.$$E9

This equation could be used when the volumes of injected and produced fluid volumes are monitored and solid volume change can be back-calculated from ion chromatography (IC) of produced fluids. If, however, the bulk volume (e.g., 4D seismic) and the pore volume were obtained from monitoring the injected and produced fluid volumes, the porosity is as follows:

$$\varphi \left(t\right)=\frac{V_{p,0}+\mathrm{\Delta }{V}_p\left(t\right)}{V_{b,0}+{\mathrm{\Delta }V}_b\left(t\right)}.$$E10

3. Volumetric strain by imposed stress

In compressive hydrostatic systems, the porous rocks deform by reducing its bulk volume. This may affect the porosity through, for example, Eq. (8). In closed systems, in which the mass and density of the minerals are conserved, the bulk volume reduction equals the pore volume reduction. This is facilitated by grains moving relative to each other, and/or by pressure solution (dissolution of stressed grain contacts and precipitation in unstressed parts of the mineral framework [6]). In open systems subjected reactive flow, both mass and density of the core material may change because of mineral reactions. To evaluate the relative importance of how evolution mechanisms of the solid volume and pore volume dictate the porosity in real systems, a rigorous definition of stresses and strains are required. The aim is to pave the way for quantitative analyses of how stresses impact strains, and how strains impact the porosity chemo-mechanical compaction.

3.1. The stress tensor in porous materials

The stress tensor describes the stresses (force per unit area) in a solid porous body. For cylindrical core plug, it is convenient to express the stress tensor $\sigma$ as

$${\sigma }_{\mathit{ij}}=\left[\begin{array}{ccc}{\sigma }_{\mathit{zz}}& {\tau }_{\mathit{zr}}& {\tau }_{\mathit{z\theta }}\\ {\tau }_{\mathit{rz}}& {\sigma }_{\mathit{rr}}& {\tau }_{\mathit{r\theta }}\\ {\tau }_{\mathit{\theta z}}& {\tau }_{\mathit{\theta r}}& {\sigma }_{\mathit{\theta \theta }}\end{array}\right].$$E11

Shear and normal components are abbreviated ${\tau }_{\mathit{ij}}$ and ${\sigma }_{\mathit{ij}}$, respectively, with $\mathit{ij}=\left\{z,r,\theta \right\}$ denoting the axial ($z$), radial ($r$), and tangential ($\theta$) direction. Compressive stresses and inward deformation are defined positive. When there is no net translational or rotational force acting in the solid body (i.e., ${\tau }_{\mathit{zr}}={\tau }_{\mathit{rz}}$, ${\tau }_{\mathit{z\theta }}={\tau }_{\mathit{\theta z}}$, and ${\tau }_{\mathit{r\theta }}={\tau }_{\mathit{\theta r}}$), only six independent stress tensor components apply. For a cylindrical core plug stressed in a triaxial cell, the tangential stress equals the radial, and the principal stress directions coincide with the imposed $z$ and $r$ directed stress such that the shear stresses are zero. The stress tensor may, therefore, be expressed through the orthogonal principal stresses vector with two components

$$\left[\begin{array}{c}{\sigma }_z\\ {\sigma }_r\end{array}\right]$$E12

In reservoir systems, however, all stress components may apply, and as such, the off-diagonal elements of the stress tensor are nonzero. However, in these cases, the stress tensor can be rotated such that the principal stress notation can be obtained. It is customary procedure to arrange the first, second, and third principal stress directions as ${\sigma }_1>{\sigma }_2>{\sigma }_3$, where ${\sigma }_1$ is typically in the vertical direction (weight dominated), and, consequently, the ${\sigma }_2$ and ${\sigma }_3$ are horizontal (often abbreviated ${\sigma }_H$ as the highest horizontal stress and ${\sigma }_h$ is the least horizontal stress). ${\sigma }_h$ and ${\sigma }_H$ depend upon Poisson ratio and tectonic regional stresses.

3.2. Effective stress

In porous rocks, it is the effective stresses introduced by [7] that drive deformation. The external load applied onto a material that consists of solids and voids is balanced by the interparticle contacts in force networks (material framework) and a fraction $\alpha$ of the pore pressure. Drained conditions apply to the cases where fluids are allowed to escape to keep the pore pressure constant (hence constant effective stress), differ from undrained conditions in which the pore pressure increases because of compaction (thereby reducing the effective stress).

Simultaneously, seepage forces arising from differences in fluid pressure expose a net force onto the solid framework (see Figure 3). In partially consolidated systems, in which the cross area is given by the sum of the consolidated area (solid–solid area) and the area of the fluid-to-solid contact area ($A_{\mathit{tot}}=A_{\mathit{fs}}+A_{\mathit{ss}}$), the force from the fluid pressure differences ($\mathrm{\Delta }{P}_f$) to the solids is given as $F_{\mathit{fs}}=\mathrm{\Delta }{P}_f\alpha A_{\mathit{tot}},$ where the fraction total area is termed as the Biot coefficient, and can be expressed as $\alpha =A_{\mathit{fs}}/A_{\mathit{tot}}$. In addition, other definitions of the Biot coefficient may also apply. In weight-dominated reservoir systems of fluid saturated rocks, the solid stress increases with the lithostatic weight. The net effective stress, that is, the stress that drive deformation is given by the differences between lithostatic pressure and the fraction $\alpha$ (the Biot coefficient) of the pore pressure

$${\sigma }_v^{\prime }={\sigma }_v-\alpha P_f\text{which is greater than}{\sigma }_H^{\prime }={\sigma }_H-\alpha P_f>{\sigma }_h^{\prime }={\sigma }_h-\alpha P_f$$E13

Figure 3.

Fluid pressure differences ($\mathrm{\Delta }{\mathrm{P}}_{\mathrm{f}}$) impart forces onto the solid framework through the fluid-to-solid contact areas, which covers only a fraction $\mathrm{\alpha }$ of the cross-section (${\mathrm{A}}_{\mathrm{f}\to \mathrm{s}}={\mathrm{A}}_{\mathrm{tot}}-{\mathrm{A}}_{\text{cons}}$).

Here, the largest and smallest horizontal stress is abbreviated with an index $H$ and $h$, respectively.

In core scale experiments, the directions perpendicular to the $z$-axis are equal, ${\sigma }_r={\sigma }_{\theta }$, implying that full description of the effective stress state of a cylindrical core experiment are given by two effective stresses, ${\sigma }_r^{\prime }={\sigma }_r-\alpha P_f$ and ${\sigma }_z^{\prime }={\sigma }_z-\alpha P_f$. The stress exerted onto the core in the radial direction is in many (not all) rock mechanical experiments performed by increasing the hydraulic confining pressure of oil surrounding the core encapsulated by a rubber or plastic sleeve while a piston placed on top of the core controls the axial stress.

3.3. Defining strain

The most commonly used definition of strain, applicable to small finite deformations, is presented here. It is acknowledged that other definitions of strain also exist in the scientific literature. The strain at any time is given by the ratio of elongation divided by the initial length

$$\epsilon \left(t\right)=-\frac{L\left(t\right)-L_0}{L_0}=-\frac{\mathrm{\Delta }L}{L_0}\text{and}{\epsilon }_{\mathit{vol}}=-\frac{\mathrm{\Delta }V}{V_0}$$E14

Stresses may deform Earth materials so that two initially orthogonal directions change by an angle $\mathrm{\Psi }$. This change in angle is related to the shear strain $\mathrm{\Gamma }$ as

$$\mathrm{\Gamma }=\frac{1}{2}\text{tan}\mathrm{\Psi }.$$E15

In three dimensions (cylindrical coordinates), the pairs of shear and normal strains are organized in the strain tensor

$$\left[\begin{array}{ccc}{\epsilon }_{\mathit{zz}}& {\mathrm{\Gamma }}_{\mathit{zr}}& {\mathrm{\Gamma }}_{\mathit{z\theta }}\\ {\mathrm{\Gamma }}_{\mathit{rz}}& {\epsilon }_{\mathit{rr}}& {\mathrm{\Gamma }}_{\mathit{r\theta }}\\ {\mathrm{\Gamma }}_{\mathit{\theta z}}& {\mathrm{\Gamma }}_{\mathit{\theta r}}& {\epsilon }_{\mathit{\theta \theta }}\end{array}\right]$$E16

Similar to the stress tensor, the shear strains balance each other (${\mathrm{\Gamma }}_{\mathit{rz}}={\mathrm{\Gamma }}_{\mathit{zr}},{\mathrm{\Gamma }}_{\mathit{\theta z}}={\mathrm{\Gamma }}_{\mathit{z\theta }},{\mathrm{\Gamma }}_{\mathit{\theta r}}={\mathrm{\Gamma }}_{\mathit{r\theta }}$), thereby, reducing the number of parameters to fully describe the deformation of a volume element in 3D from nine to six parameters. In addition, for isotropic materials, the principal strains can also be found by rotating the strain matrix, the same way as the stress matrix, such that the off-diagonal elements vanish (${\mathrm{\Gamma }}_{\mathit{ij}}=0$). In addition, the radial and tangential strains are equal, such that the strain vector for cylindrical cores:

$$\left[\begin{array}{c}{\epsilon }_z\\ {\epsilon }_r\end{array}\right]$$E17

To estimate the porosity evolution, bulk volumetric strain has to be used. The volumetric strain equals the change in volume divided by the initial volume, which is the first strain invariant, given by ${\epsilon }_{\mathit{vol}}=-\mathrm{\Delta }V/V_0=\mathrm{Tr}\left({\epsilon }_{\mathit{ij}}\right)$. The volume strain remains unchanged upon coordinate change (i.e., the volume is the same irrespective of which coordinate system is used). Depending upon the geometry of the setup, the way in which strain measurements and hence the strain tensor components will vary. For cylindrical geometries, in which the volume of a cylinder is given by $V=\pi D^{2}L/4$, where $D$ is the diameter and $L$ is the length, the volumetric strain can be calculated from the radial and axial strains

$${\epsilon }_{\mathit{vol}}={\epsilon }_z+2{\epsilon }_r+2{\epsilon }_z{\epsilon }_r+{\epsilon }_r^2+{\epsilon }_z{\epsilon }_r^2$$E18

If the length and diameter of cylindrical cores are being measured continuously, then the volumetric strain can be estimated. Typically, for small strains, the second and third order terms can be omitted, hence, ${\epsilon }_{\mathit{vol}}\simeq {\epsilon }_z+2{\epsilon }_r$.

4. Partitioning time-independent and time-dependent deformation

The volumetric strain can be split into an immediate strain, occurring when the effective stress is being changed, and time-dependent deformation. The two cases are presented briefly in the following sections, even though this is a large area of research. For the time-independent case, Hooke’s law is described before nonlinear models are presented, followed by a short note on plasticity and other failure mechanisms before time-dependent models are described.

4.1. Elastic strain: linear elasticity

Hooke’s law is the simplest relation to describe the relation between the stress-strain tensors. It assumes that the deformation is immediate, linear, and reversible. In continuous media, for small stress and strain increments in the linear limit, the ${\epsilon }_{\mathit{ij}}$ and ${\sigma }_{\mathit{ij}}$ are described by the compliance (stiffness) fourth order tensor $c_{\mathit{\text{ijkl}}}$. In 3D systems, it consists of 81 real numbers, and the tensorial equation attains the compact form ${\sigma }_{\mathit{ij}}=c_{\mathit{\text{ijkl}}}{\epsilon }_{\mathit{kl}}$, where the indexes $i,j,k,l$ represent the three spatial dimensions $\left[x,y,z\right]$ in Carthesian co-ordinate systems, and $\left[z,r,\theta \right]$ in cylindrical systems. In the case, when the rotational forces balances, which applies to most continuum mechanical cases, the number of stiffness parameters describing the stress-strain relation reduces to 27. In the case of isotropic materials, the number of elastic parameters that describe the stress-strain relation of a volume element along the principal directions is further reduced to the Young’s modulus ($E$) and Poisson’s ratio ($\nu$) via the matrix equation

$$\left[\begin{array}{c}{\epsilon }_z\\ {\epsilon }_r\\ {\epsilon }_{\theta }\end{array}\right]=\frac{1}{E}\left[\begin{array}{ccc}1& -\nu & -\nu \\ -\nu & 1& -\nu \\ -\nu & -\nu & 1\end{array}\right]\left[\begin{array}{c}{\sigma }_z^{\prime }\\ {\sigma }_r^{\prime }\\ {\sigma }_{\theta }^{\prime }\end{array}\right]$$E19

By adding up the three equations expressed in the matrix form earlier

$$\left({\epsilon }_z+{\epsilon }_r+{\epsilon }_{\theta }\right)E=\left(1-2\nu \right)\left({\sigma }_z^{\prime }+{\sigma }_r^{\prime }+{\sigma }_{\theta }^{\prime }\right),$$E20

we may use this equation to define the bulk modulus in hydrostatic tests. When omitting higher order terms in the volumetric strain (Eq. (18)), the left hand side of Eq. (20) equals the volumetric strain. For hydrostatic tests, in which the stresses in all spatial directions equal, ${\sigma }_z^{\prime }={\sigma }_r^{\prime }={\sigma }_{\theta }^{\prime }={\sigma }_p^{\prime }$, then Eq. (20) simplifies to

$$\frac{E}{1-2\nu }{\epsilon }_{\mathit{vol}}={\sigma }_p^{\prime }\to K{\epsilon }_{\mathit{vol}}={\sigma }_p^{\prime }$$E21

In Eq. (21), the bulk modulus ($K$) is defined. ${\sigma }_p^{\prime }$ is frequently used to define the hydrostatic effective stress. For nonhydrostatic triaxial tests, where ${\sigma }_z>{\sigma }_r$, and ${\sigma }_r^{\prime }={\sigma }_{\theta }^{\prime }$ and ${\epsilon }_r={\epsilon }_{\theta }$ Hooke’s law in Eq. (19) simplifies to

$$\begin{array}{c}E{\epsilon }_r=\left(1-\nu \right){\sigma }_r^{\prime }-\nu {\sigma }_z^{\prime }\\ E{\epsilon }_z={\sigma }_z^{\prime }-2\nu {\sigma }_r^{\prime }\end{array}$$E22

4.2. The effective stress changes that drive deformation

Within the elastic domain, any change in the effective stress drive deformation in the sample, from here on abbreviated with the $\delta$-symbol used to rewrite Hooke’s law at quasistatic changes. The $\delta$-symbol is used to identify the variables that are changing during for example, a loading sequence

$$\begin{array}{c}\mathit{E\delta }{\epsilon }_r=\left(1-\nu \right)\delta {\sigma }_r^{\prime }-\mathit{\nu \delta }{\sigma }_z^{\prime }\\ E{\mathit{\delta \epsilon }}_z={\mathit{\delta \sigma }}_z^{\prime }-2\nu {\mathit{\delta \sigma }}_r^{\prime }\end{array}$$E23

In Eq. (23), the underlying assumption is that Young’s modulus and Poisson ratio remain fixed. Furthermore, when pore pressure is included, the effective stress changes due to both axial and radial stress and pore pressure

$${\sigma }_r^{\prime }=\delta {\sigma }_r-\mathit{\alpha \delta }{P}_f\text{and}\delta {\sigma }_z^{\prime }=\delta {\sigma }_z-\mathit{\alpha \delta }{P}_f$$E24

It is assumed that the Biot stress coefficient remains fixed during loading. Using these definitions into Eq. (32) enables us to fully describe the relation between the stress, pore pressure, and strain in Hooke’s law

$$\begin{array}{c}\mathit{E\delta }{\epsilon }_r=\left(1-\nu \right)\delta {\sigma }_r-\nu {\mathit{\delta \sigma }}_z+\left(2\nu -1\right)\mathit{\alpha \delta }{P}_f\\ E{\mathit{\delta \epsilon }}_z={\mathit{\delta \sigma }}_z-2\mathit{\nu \delta }{\sigma }_r+\left(2\nu -1\right)\mathit{\alpha \delta }{P}_f\end{array}$$E25

4.3. Plasticity and irreversible deformation

For a highly porous chalk, a nonzero component of the observed strain is always irreversible when the load is released as exemplified by [8] and [9] where irreversible plasticity is seen also within the ‘elastic’ phase of the $\mathit{QP}$-plot (Figure 4). As such, loading and unloading may display history dependence in the elastic parameters. This may be caused by the way in which the porous material is being held together and the relative importance of the different forces responsible for determining the stiffness of the chalk. Now, at increasing stresses beyond ‘elasticity’, the type of irreversible deformation that develops, depend on the state of stress, as illustrated in Figure 4. Here, the mean effective stress is plotted on the horizontal axis and the deviatoric stress along the vertical axis. Considering cylindrical cases, the deviatoric stress equals $\mathrm{Q}={\mathrm{\sigma }}_{\mathrm{z}}-{\mathrm{\sigma }}_{\mathrm{r}}$ and the mean effective stress is $\mathrm{P}=\left({\mathrm{\sigma }}_{\mathrm{z}}+2{\mathrm{\sigma }}_{\mathrm{r}}\right)/3-\mathrm{\alpha }{\mathrm{P}}_{\mathrm{f}}.$ Since both $Q$ and $P$ are invariant, the results of core data can be used at any case in which the material is the same: (1) for hydrostatic systems, $\mathrm{Q}=0$, and pore collapse occurs when the mean effective stress exceeds a certain threshold; (2) tensile fractures develop at negative values of $P$ which can be found for high fluid pressures, or in Brazilian tests; and (3) shear failure occurs when the deviatoric stress exceeds a certain value. For frictional materials, it is typical to observe that the deviatoric stress required to induce shear failure is increasing with increasing mean effective stress. For Coulomb materials, this relation is proportional, and the slope is related to the frictional coefficient. Chalks have been found to be satisfactory described with such a frictional coefficient, while clays behave differently. The way in which the irreversible deformation affects the porosity evolution differs from case to case. Within shear zones, the porosity may both increase, because of dilation and de-compaction when tightly packed grains reorganize or reduce because of grain crushing when the imposed forces exceed a certain level.

Figure 4.

The failure envelope (solid line) shows that at which stresses plasticity and irreversible deformation occur (numbers are not applicable). Mean effective stress ($\mathrm{P}$) and $\mathrm{Q}$ is the deviatoric stress, hence the end cap depends solely upon the material and not the geometry.

5. Time-dependent pore volume reduction processes and compaction

To understand how the mechanical and chemical processes affect the porosity during pore collapse, it is important to take a closer look at how the observed bulk strain can be partitioned. In this section, we consider the simplest possible partition (below) in which the overall strain is partitioned additively into a solid volume and a pore volume component. The relative importance of these mechanisms may be found by the analysis of quantitative measurements of the bulk volume change and the change in the solid volume due to the dissolution/precipitation as the mineral mass and density change over time, while grains reorganize, crush, and solid contacts evolve. In this case, the observed volumetric change can be partitioned additively via

$${\epsilon }_{\mathit{vol}}={\epsilon }_{\mathit{\text{pore}}}+{\epsilon }_{\mathit{\text{solid}}}$$E26

This does not imply that imply that cross terms do not exist in which: (a) the rate of pore volume reduction is sensitive to the reduction in solid volume and (b) how the solid volume rate may depend on how grains reorganize to change the flow pattern and potentially expose new fresh mineral surfaces to the reactive brine. It is likely to assume that based on the accelerated strain presented in [10] minute changes to the solid volume increase the rate of pore collapse (also seen in [4]).

Given the simple partitioning above, a model can be developed to describe the observed creep curve with a few physical parameters (see Eqs. (20)–(23) in [10]). In this model, overall volumetric strain is additively partitioned into a pore and solid volume component in which the pore volume equals, $V_p=\varphi V_b$. Extending the rate of change in bulk volume is by the pore volume change (via using the product rule) and the solid volume change rate

$$\frac{d{V}_b}{\mathit{dt}}=\frac{d{V}_p}{\mathit{dt}}+\frac{d{V}_s}{\mathit{dt}}=\varphi \frac{d{V}_b}{\mathit{dt}}+V_b\frac{\mathit{d\varphi }}{\mathit{dt}}-\beta$$E27

The solid volume rate is assumed to be constant ($\beta$, in cm³/day determined from ion chromatography data). For Mons chalk at 130°C and 92°C at 1 PV/day of 0.219 MgCl₂ brine, the solid volume changes approximately by 0.01 and 0.005 cm³/day, respectively. The porosity reduction rate can be proportional to porosity

$$\frac{\mathit{d\varphi }}{\mathit{dt}}=-\xi {\left(\varphi -{\varphi }_c\right)}^n$$E28

where $\xi$ is the proportionality constant, ${\varphi }_c$ is a terminal porosity (grain reorganization cannot continue until zero porosity), and the power $n$ is used to model nonlinear behavior. For simplicity, if we assume $n=1$ and ${\varphi }_c=0$, the volumetric strain is explicitly given as

$${\epsilon }_{\mathit{vol}}\left(t\right)=\frac{\left(1-{\varphi }_0\right)e^{\mathit{\xi t}}}{e^{\mathit{\xi t}}-{\varphi }_0}-\frac{\mathit{\beta t}}{V_{b,0}\left(1-e^{-\mathit{\xi t}}\right)}-1$$E29

This model takes the initial porosity (${\varphi }_0$) and bulk volume ($V_{0,b}$), while the porosity rate constant $\xi$ is a free variable.

The mathematical models aimed to match observed creep data have a long history, and several, more or less physically based models were reported. Generally, these models do not consider the underlying solid and pore volume contribution, but may still satisfactorily match the observed strains. Three models that have been used are:

Power law with cut-off:

$${\epsilon }_P=A{t}^B{e}^{-\raisebox{1ex}{$t$}\!\left/ \!\raisebox{-1ex}{$t_0$}\right.}$$E30a

De Waal [11]:

$${\epsilon }_{\mathit{dW}}=A\text{log}\left(\mathit{Bt}+1\right)$$E30b

Griggs [12]:

$${\epsilon }_G=A\text{log}\left(t+1\right)+\mathit{Bt}$$E30c

The model parameters $\left(A,B\text{and}{t}_0\right)$ are found when the residual strain $\mathit{RES}=\frac{1}{N}{\sum }_n\mid {\epsilon }_{\mathit{exp}}-{\epsilon }_{\mathit{\text{model}}}\mid$ is minimized.

5.1. Pore collapse and grain reorganization: the constant solid volume case

The movement of grains relative to each other at high-mean effective pressures causes pore volumes to collapse. It has been experimentally verified that for chalks, the rate of compaction may sometimes accelerate when the fluid composition change, a process termed water weakening. Water weakening has been used to understand reservoir processes [13, 14, 15] and to interpret core experiments as exemplified in Figure 5, where an additional strain of $\sim$1% is seen at the first days of seawater (SSW) flow (approximately 2 pore volumes). This process cannot be attributed to chemical reactions leading to solid volume changes, since the ions in the produced effluent water are inadequate to cause any solid volume change from the mass loss and increased density often seen in these cases when chemical reactions occur. As such, the additional bulk strain is caused by pore collapse. This does not exclude how long-term chemical reactions can weaken rocks over longer time periods when porous chalk cores are continuously flooded.

Figure 5.

Axial creep strain over time at uniaxial strain condition performed on a Kansas chalk sample. The injection of seawater (SSW) leads to accelerated creep. The accelerated creep period is associated with the loss of sulfate ions in the effluent samples. Mixing CO₂ into the SSW from 120 days and onwards does not induce additional strain (from [20]).

To understand the immediate additional deformation (i.e., 1–2 days corresponding to a flow of 1-2 pore volumes as seen here), a grain-level approach is required. The grain-grain friction controls and cement bonds binding neighboring grains together control the relative movement of grains. Friction between grains is given by the frictional coefficient times the normal force, $F_{\mathit{\text{fric}}}=\mu F_N$. The normal force arises from the externally imposed load and the attractive Van der Waal forces that induce the cohesive forces grains. This has been shown to be reduced by negative disjoining pressures in the overlapping double layer between adjacent mineral grains when surface-active divalent ions adsorb onto the charged chalk surfaces [16, 17, 18, 19].

5.2. Pressure solution and other grain-reorganization mechanisms

Pressure and temperature are the state variables that control the Gibbs chemical potential [21]. During diagenesis and burial, the chemical stability of mineral phases is altered as the temperature, hydrostatic and lithostatic pressure increases. Pressure solution of stressed grain contacts, and precipitation in unstressed parts of the rock framework, have been used as one of the primary rock-forming mechanisms during diagenesis. Pressure solution can occur in closed systems, in which the overall mass and density remain fixed. For high Biot coefficients, the local stress at particle contacts may become significant (see Figure 7 in [10]), and thus, a stress-dependent production of Ca-ions is observed where more Ca-production for high stress than low stress (see Figure 6 acquired from [10]).

Figure 6.

(a) Outlet cation concentration of 0.219 M MgCl₂ (dashed line, same ion strength as seawater) flooded through chalk from the Obourg saint vast formation (Mons, Belgium) at 130°C for 0.5, 3.5, and 12.3 MPa effective stress. (b) Calcium production at varying stresses at 130°C and 92°C. The amount of Ca ions in the produced effluent depends on time, temperature, and stress (acquired from Figure 6 in [10]).

It has been a long-standing discussion how to mathematically describe the relevant thermodynamic pressure for accurate determination of the chemical potential from the stress tensor. Several candidates coexist. The stress tensor in reservoir systems (and core scale experiments) depends on the weight of the overburden (lithostatic weight), side stress (tectonic forces), pore pressure, and the Biot coefficient. The simplest of determining the thermodynamic pressure is using the pore pressure. This way of thinking may seem reasonable at first glance since it is at the interface between the solid and the fluid where the chemical reactions occur. Simultaneously, at the rock-fluid interface, the stresses through the solid framework could also play a role in determining chemical solubility. In that case, the continuum mechanics provide a range of choices for calculating the thermodynamic pressure: (1) the average compressive stress (i.e., the first invariant of the solid framework stress tensor), (2) the principal stresses, thereby leading to different solubility in the different spatial directions. In sedimentary systems, this would often lead to enhanced solubility in the vertical direction as the first principal stress direction is vertical. This may explain the formation of the horizontal stylolites that are sometimes found in calcitic, carbonate, and limestone rocks [6]. (3) The relevant thermodynamic pressure could be related to the stress gradients that have been observed throughout porous materials, termed force chains. At grain-grain contacts, through which the externally imposed loads are being carried, the stresses can be significantly higher than the average. In these regions, the solid-solid stress is given by ${\sigma }_{\mathit{ss}}^{\prime }={\sigma }^{\prime }/\left(1-\alpha \right)$, where ${\sigma }^{\prime }$ is the effective stress and $\alpha$ is the Biot stress coefficient. As such, for unconsolidated sands and calcitic mudstone, in which $\alpha >0.9$, this fraction is significant and may be responsible for additional calcite dissolution [10]. As has been shown previously (in e.g., [18, 22] and also before that), the contact area ratio is linked to the Biot stress coefficient ($\alpha$).

Even though pressure solution is a process of chemical nature, it does not necessarily change the solid volume since the mass can be conserved (closed system, i.e., no larger scale mass flow) and the same mineral phase is precipitated as the one dissolved (i.e., same density). In that sense, pressure solution contributes to pore volume reduction rather than the solid volume in the strain partitioning presented here. Hence, pressure solution may fall under mechanical compaction even though the underlying mechanisms of pressure solution are chemically driven.

6. Time-dependent solid volume evolution mechanisms

In open, nonequilibrium systems with rock-fluid interactions, the solid volume is subject to change. It has been shown in a range of experiments how additional strain is accumulated during compaction at constant stress conditions when reactive brines are injected [4, 23]. The solid volume varies when solid mass ($M_s$) and mineralogical density (${\rho }_s$) change

$$V_s=\frac{M_s}{{\rho }_s}$$E31

The change in solid volume may be evaluated by

$${\mathrm{\Delta }V}_s\left(t\right)=\frac{M_s\left(t\right)}{{\rho }_s\left(t\right)}-\frac{M_{s,0}}{{\rho }_{s,0}}$$E32

Here, the solid volume change is given by the difference between the ratio of the mass and density at a given time and the values before chemo-mechanical processes initiated. The evolution of the solid mass over time is given by the difference between the chemical mass flux in and out of the system, and density changes as new minerals precipitate.

6.1. Mass transfer in open systems

When fluids continuously flow and react with the rock, the mass (and hence the solid volume) changes. The chemical flux can be monitored by evaluating the effluent concentration through difference between the ion concentrations in and out of the volume element (Figure 6a, b). This volume element may, in some cases, be between an injector and a producer in an oil field, or a core scale experiment in the laboratory [3]. The concentration of ions can be measured using ion chromatography, and over a time interval $\mathit{\delta t}$ the difference in mass is given by

$$\frac{{\mathrm{\delta }M}_s}{\mathit{\delta t}}={\sum }_j\left(c_{\mathit{\text{in}},j}-c_{\mathit{\text{in}},j}\right)q{m}_j$$E33

In Eq. (33), the factor $\left(c_{\mathit{\text{in}},j}-c_{\mathit{\text{in}},j}\right)$ is the difference in the ion concentration of chemical species $j$ (mole/L), $q$ is the flow rate (L/day), and $m_j$ is the molar mass of species $j$ (g/mole). Hence, the term ${\mathrm{\delta }M}_s/\mathit{\delta t}$ is given in g/day. The overall mass is estimated by summing over all measured ions in the chemical interaction, giving a unit (g/day), which can be used further. The total mass evolution of each species is determined by integration. Assessing rock-fluid interactions to real cases, for example during, seawater flooding of the Ekofisk field (North Sea, Norway), chemical reactions have been observed. Here, dissolution of 1–2 wt. % is anticipated from the analysis of the produced water [24, 25, 26].

6.2. Method to quantify the solid volume evolution

In Eq. (32), the change in solid volume depends on the change in both mass and in density as the minerals dissolve and precipitate. The overall mineral density, as $n$ minerals dissolve/precipitate in time is given by

$${\rho }_s=c_1{\rho }_1+c_2{\rho }_2+\dots +c_n{\rho }_n$$

$$1=c_1+c_2+\dots +c_n$$E34

As the concentrations of different minerals vary, the changes to ${\rho }_s$ can be estimated.

It is not always the case that a detailed kinetic chemical model exist tuned to take into account how different mineral mixtures react with fluids in each different case. If the overall density before (${\rho }_{s,0}$) and after (${\rho }_{s,f}$) the chemical experiment are known (from e.g., pycnometry) a reduced mass parameter ($\tilde{m}$) ranging from 0 to 1 can be defined from the initial and final mass using mineral $k$ ($M_{k,0}$ and $M_{k,f}$, respectively)

$$\tilde{m}\left(t\right)=\frac{M_k\left(t\right)-M_{k,0}}{M_{k,f}-M_{k,0}}$$E35

Then, the density at any given time may be estimated using

$${\rho }_s\left(t\right)=\tilde{m}\left(t\right){\mathrm{\rho }}_{\mathrm{s},\mathrm{f}}+\left(1-\tilde{m}\left(t\right)\right){\rho }_{s,0}$$E36

7. Predicting dynamic porosity evolution: an illustrative example

An example of a dynamic porosity development analysis is presented here based on material published in 2015 [3]. The results of an experiment performed over 1090 days where a Liegè (Belgium) chalk sample was exposed to hydrostatic stress (11.1 MPa, approximately 5 MPa above yield) and continuous flow of 0.219 M MgCl₂ (33 and 99 cm³/day, pore pressure 0.7 MPa and 130°C). Basic sample measurements were performed of dry/saturated mass, pore volume, solid and bulk volume and hence porosity before and after test, hence the mineral density estimated and confirmed using He-pycnometry (Table 1). The bulk core volume was reduced more than 15% and mineral mass is reduced by more than 18% while the density is increased from 2.7 to 2.9 g/cm³. The stresses are sufficient to induce pore collapse for these chalks, and the 0.219 M MgCl₂-brine (with equal ion strength as seawater) induced dissolution of the calcium carbonate and precipitation of denser Mg-bearing carbonates (e.g., magnesite and dolomite).

	Before test	After test (1090 days)	Change
Dry mass (on scale)	125.57 g	102.64 g	−22.93 g
Wet weight (saturated)	158.56 g	126.34 g	−32.22 g
Pore volume	32.99 cm^3	23.71 cm^3	−9.28 cm^3
Solid volume	46.85 cm^3	35.53 cm^3	−11.32 cm^3
Bulk volume	79.84 cm^3	59.23 cm^3	−20.61 cm^3
Mineral density (saturation and pycnometer)	2.68 and 2.70 g/cm^3	2.89 and 2.90 g/cm^3	0.21 and 0.20 g/cm^3
Porosity (saturation and pycnometer)	41.3 and 41.7%	40.0 and 40.1%	−1.3 and −1.6%

Table 1.

Basic measurements of the core before and after the 1090 days long-term test [3]. Figure 1 displays SEM images of the core material before and after the flow-through test.

During the test, the axial deformation and the ion concentration of the effluent fluids were measured. In Figure 7a, the Mg and Ca ion concentrations are measured through time. These measurements can be used to find the production rate in g/day using Eq. (33), seen in Figure 7b. A trebling of the inlet flow rate leads to more than a doubling in the calcite dissolution. The mass evolution is used to estimate the dynamic change in density using Eq. (36), and is then combined to estimate the solid volume as function of time as seen by the dotted line in Figure 7c plotted together with the bulk volume, estimated from the axial strain. Thus, the pore volume can be estimated (dashed line in Figure 7c). As can be seen, the pore volume is reduced when bulk compaction dominate the overall process until 200 days. A typical observation from primary creep experiments is that the overall creep rate decreases with time. After 200 days, when the compaction rate has reduced, the flow rate was increased thereby increasing the rate at which dissolution/precipitation occurs, see Eq. (33) where the flow-rate dependency is explicitly shown. At this point of time, the overall porosity dynamics change. In the initial compaction-dominated regime, the overall porosity reduced to a value of as low as 33%, and afterwards it starts increasing (solid line in Figure 7d). At approximately 400 days, the flow rate is then reduced again and the rate of change in porosity is changing accordingly.

Figure 7.

(a) Ion chromatography of the produced ion concentration of mg and Ca throughout the test. Mg is retained in the core while Ca is produced. (b) Calculated production rate of mg (solid) and Ca (dashed) using Eq. (32), (33) and (36), while (c) displays the total, pore, and solid volume evolution. (d) Observed volumetric creep (dashed line) and estimated porosity evolution as the relative importance of bulk compaction and dissolution/precipitation change.

From 900 days and onwards, the Ca that was initially found within the core had been produced, after the solid volume was interpreted to be constant and the bulk compaction is facilitated by pore volume reduction, and hence the porosity is decreased to 40.1%.

In the experiment presented here, both pore and solid volume are subject to change. Since only the bulk volume or the solid volume could be determined from axial strain IC data, respectively, the pore volume was determined. As is exemplified in the presented experiment, the porosity evolution dynamics display a complex behavior because of the reduction in pore volume and solid volume. Their rate depends upon stress, the way in which deformation is accumulated and the rate of dissolution/precipitation.

8. Summary

Porosity is an important parameter for understanding the diagenetic processes and petrophysical reservoir systems. Its importance to the mechanical stiffness and strength of porous rocks, and to the resource potential, and rate of hydrocarbons produced from reservoirs is evident. The porosity is a dynamic parameter from the strain and chemical reactions from injection of fluids out of equilibrium with the host rock (e.g., seawater brines at elevated temperature in chalks) that induce additional deformation over time. The adsorption of surface-active ions leads to alterations in the forces binding grains together, leading to instantaneous additional deformation.

To understand quantitatively how porosity changes dynamically through time, there are series of processes that needs to be incorporated. This chapter presents some of the ways in which the bulk strain can be partitioned into elastic/plastic components, time-dependent, and time-independent components, and solid volume and pore volume processes. For chalks, the dynamic porosity evolution depends on the relative importance of the different processes at play, that again are functions of the stress, strain, temperature, flow rate, and fluid chemistry.

The methods presented here do not cover all possibilities for porosity evolution determination depending upon measurements that are available. When the bulk volume strain and chemical composition of the effluent fluids are known, the following porosity evolution model applies

$$\varphi \left(t\right)=\frac{{\varphi }_0-{\epsilon }_{\mathit{vol}}\left(t\right)-\left(\frac{M_s\left(t\right)}{{\rho }_s\left(t\right)}-\frac{M_{s,0}}{{\rho }_{s,0}}\right)/V_{b,0}}{1-{\epsilon }_{\mathit{vol}}\left(t\right)}$$E37