# Sandstone compaction by intergranular pressure solution

F. K. Lehner, Y. M. Leroy

► To cite this version:

F. K. Lehner, Y. M. Leroy. Sandstone compaction by intergranular pressure solution. Mechanics of fluid saturated rocks, pp.115-168, 2004, 10.1016/S0074-6142(03)80019-0. hal-00112505

HAL Id: hal-00112505

https://hal.science/hal-00112505v1

Submitted on 17 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Sandstone Compaction by Intergranular Pressure Solution

Florian Lehner<fnref n="1" />

Yves Leroy<fnref n="2" />

## 3.1 Introduction

It has long been known that the mineral matrix of porous rocks tends to dissolve in pore water preferentially along highly stressed, fluid-permeated grain-to-grain contacts. The dissolved material is then transported by molecular diffusion and advective fluxes to nearby or possibly distant precipitation sites, typically the free faces of mineral grains where they meet the pore fluid. This process can dramatically change the rock fabric, as is illustrated by the example shown in Figure 3.1. It is considered one of the principle processes of the early diagenesis of sedimentary rocks (Bathurst, 1958; Wanless, 1979; Houseknecht, 1988; Houseknecht and Hathon, 1987; Tada and Siever, 1989), and it is also a process by which substantial macroscopic creep strains may accumulate. The implied low-temperature ductile deformation mechanism is most often termed “pressure solution” by geologists (Rutter, 1983), but is also known as solution-precipitation creep or solution transfer creep (Durney 1972).

That rocks deform by such a chemomechanical process was observed by geologists on a somewhat larger scale as early as 1863, when Sorby first explained the phenomenon of pitted pebbles (see Figure. 3.2) in terms of a selective removal of material from the more soluble of two impinging pebbles (Mosher 1976, 1981; for references to early work on pressure solution and related theories, see especially the review article by Durney 1978). Sorby also hypothesized that the (limestone)

<footnote n="1">¹Department of Geology and Palaeontology, University of Salzburg, Hellbrunnerstrasse 34, 5020 Salzburg, Austria. lehner.fk@sbg.at</footnote>

<footnote n="2">²LMS, Ecole Polytechnique, 91128 Palaiseau Cedex, France. leroyy@lms.polytechnique.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 3.1** ► Photomicrographs showing intergranular pressure solution (IPS) in a sandstone. (a) Pore-filling calcite cement has preserved original quartz grains by inhibiting IPS. (b) Sample from same depth as (a) without cement (dark areas are porosity), but some quartz overgrowth, showing significant IPS. (Courtesy Houseknecht and Hathon 1987.)</figcaption></figure>

material he studied had dissolved preferentially “where the pressure is greatest, and crystallized, where it is least.”

Characteristic features of pressure solution are indeed the microscopic (grain scale) and macroscopic solution seams or sutures, which tend to assume a preferred orientation perpendicular to the largest compressive stress. These appear most conspicuous in the form of stylolites, the digitated solution seams most frequently formed in porous, water-saturated carbonate rocks (Dunnington 1954, 1967; Glover 1968; Carrio-Schaffhauser et al. 1990). Stylolites owe their name to the stylo-like shapes that are seen in cases of extreme column-and-socket interdigitations, when the solution surface is exposed in three dimensions. These solution seams appear in cross section as serrated veins and are easy to spot by their dark color, which they owe to less soluble, usually clay-type minerals that are left behind
by the dissolved rock matrix. The stylos, columns, or teeth of a stylolite can reach
lengths of the order of meters and the seam thickness itself can reach similar
dimensions in rare cases, suggesting an enormous reduction in layer thickness by
dissolution. The stylolitic solution seams that can be spotted in everyday life in
marble tables or floors have typical seam thicknesses of the order of a millimeter
and typical overall lengths (in cross section) of the order of a meter.

It is not altogether clear what causes the solution process to localize macro-
scopically along discrete and often quite evenly spaced, preferentially oriented
solution seams, but a process of “anticrack growth” (Fletscher and Pollard, 1981)
and propagation (Carrio-Schaffhauser et al., 1990) appears to be involved that
remains to be fully explained. Stylolytic solution interfaces are found also on the
scale of individual grains, for example, in sandstones (Houseknecht and Hathon,
1987), where variations in solubility across individual grains or the presence of
insoluble minerals might explain the sutured morphology of the grain-to-grain
contacts (see Figure 3.1). On a grain scale, intergranular pressure solution (IPS)
therefore occurs by necessity in a highly localized fashion. On a coarser macro-
scopic scale the phenomenon may, however, affect substantial layer thicknesses in
a pervasive and seemingly continuous fashion, apparently depending on composi-
tional factors (Wanless, 1979). The distinction between discontinuous and contin-
uous modes ultimately becomes a question of the scale at which the phenomenon
is viewed.

To a theorist, these observations would suggest a two-scale approach toward a
final macroscopic theory, and the principle goal of this chapter is indeed to obtain a
pressure solution “creep law” by such an approach. Our points of departure are the
phenomena seen in Figure 3.1 and Figure 3.2. These suggest to us the concepts of an
*aggregate* of grains and that of a fluid-infiltrated, *grain-to-grain contact*, concepts
on which we shall base the *grain-scale model* of Section 3.2. Next, the transition
from the microscopic grain scale to the macroscopic continuum scale must be
accomplished. In general, this will involve a difficult *averaging step*, to which we
shall give some attention, without actually carrying out the necessary computations
for the general case of a realistic grain fabric. In Section 3.3 we circumvent the
difficulty by assuming a simple periodic arrangement of grains, as is often done.
This can be justified in view of the fundamental nature of certain questions that
a model should answer before addressing the purely geometrical complexity of a
realistic grain fabric. A further goal of future theoretical developments will be the
identification of parameters in a more complete theory—parameters critical to the
appearance of macroscopically localized modes of pressure solution.

Figure 3.1 also illustrates the following important point. If a porous and per-
meable rock mass is undergoing pervasive IPS, then any volume element of it con-
stitutes an open thermodynamic system, allowing soluble matter to be exported
from dissolution sites or imported for deposition at precipitation sites. In such a
context, pressure solution will be of interest to geologists primarily as a process
that can explain large-scale mass movements, accounting for a substantial loss in
layer thickness at a location and providing a source of cement at a different location
<figure><img src="image_3.png" /><figcaption>**Figure 3.2** ► Pitted quartzite pebbles. (Courtesy Mosher 1976)</figcaption></figure>

(Tada and Siever, 1989). But because the average path length from dissolution to reprecipitation sites may shorten, depending on the balance of relevant rate parameters, down to the size of the smallest representative elementary volume (REV) of the rock, the solid skeleton within an REV can eventually also behave like a closed system on a macroscopic scale. Material dissolved at grain-to-grain contacts will then be transported by grain boundary diffusion and diffusion through an open pore space to nearby precipitation sites at low-stress pore walls, giving rise to a characteristic microtexture with overgrowth of the original grains in so-called pressure shadows. This mass transfer from high-energy dissolution sites to low-energy precipitation sites enables an aggregate of grains to accumulate macroscopic creep strains in response to an applied load. Under closed-system conditions, the process closely resembles a type of grain boundary diffusion creep known as Coble creep in the materials science literature (Coble, 1963). Under open-system conditions, the coupling of deformation to long-range solute transport will complicate the theoretical description of creep, since now the solute concentration will enter as an additional dependent variable into the problem. Moreover, the local pore-scale kinetic processes of intergranular dissolution, grain boundary diffusion, and free-face precipitation along the pore walls will each be governed by distinct characteristic times that may differ substantially from the relevant time scales of macroscopic advective–diffusive transport.

This chapter is intended to provide the theoretical means for exploring prob-
lems of stress-sensitive dissolution and coupled macroscopic deformation and
transport. Choosing a simple model system for a grain-scale description of IPS,
we discuss the thermodynamic foundations of the subject of pressure solution in
Section 3.2. The resulting pore-scale model is shown to furnish pertinent source terms in the macroscopic mass balance equations that are introduced in Section 3.3 along with other elements of a macroscopic theory of compaction creep by IPS. This macroscopic description is partially developed to provide interested readers with some background. In Section 3.4, the theory is further specialized to yield descriptions of essentially one-dimensional boundary value problems for a closed and an open system. As a problem of the latter kind, we then study the compaction of a thick sedimentary layer, embedded in fixed hydrostatic pressure and geothermal temperature fields, in which IPS is suddenly “switched on.”

## 3.2 Grain-Scale Model of IPS

### 3.2.1 Model Assumptions

As shown in Figure 3.3, we consider a nominally flat contact between two identical spherical grains composed of a homogeneous solid substance. The solid phase is in contact with its own aqueous solution, the pore fluid, which is at a uniform

<figure><img src="image_4.png" /><figcaption>**Figure 3.3** ► IPS affecting two identical spherical grains of initial diameter *d*; nominally flat, fluid-permeated contact zone shown in magnification. Current grain size normal to contact is *l*<sub>n</sub>*d*.</figcaption></figure>
pressure $p$. Both solid and fluid phase are at a uniform temperature $T$. The two grains shown in Figure 3.3 belong to an assembly of grains, left unspecified for the moment, that together with the pore fluid sustains compressive macroscopic normal stresses that exceed the pore-fluid pressure in magnitude in all directions. There exists therefore a compressive intergranular normal stress that tends to press the grains together and is assumed to enhance the solubility of the grains in any aqueous pore fluid present at their contact.

On a finer scale, the nominally flat grain-to-grain contact is visualized here as a “contact zone” of some finite thickness $\delta$. This thickness may vary between a few nanometers for “clean” contacts between single crystals of halite or quartz and tens of micrometers for intergranular solution seams that contain impurities shed by the dissolving grains and other less soluble minerals. We further stipulate that this contact zone remain, in a roughened state during active IPS, such that adjacent grains maintain solid-solid contact across insular ridges that are separated by fluid-filled grooves (Raj, 1982; Lehner and Bataille, 1984/85; Spiers and Schutjens, 1990; Cox and Paterson, 1991; Schutjens and Spiers, 1999). In plan view such a grain boundary will then exhibit a characteristic island-and-channel structure, whereas in cross section the contact zone forms a thin, porous, and permeable interlayer that is wetted and invaded by the pore fluid. Any dissolved solid material can thus be carried by diffusive transport through this grain boundary fluid phase from interior points of a contact zone toward its rim where it meets the interconnected pore space.

Assuming now that the grain centers converge normal to the contact and ignoring small elastic shape changes, it follows that the rate of grain convergence must equal the (uniform) jump in the normal components of the grain velocities across the contact interface. Accordingly, $\llbracket v_n \rrbracket = v_n^+ - v_n^- \le 0$ (see Figure 3.3) and the uniform solution rate (in kg/m²/s) along the contact is given by $\rho^s \llbracket v_n \rrbracket = \rho^s \dot{l}_n d$, where $d$ is the initial grain diameter and $l_n d$ measures the current distance between two parallel contacts of the same grain. If $J_r^{gb}$ denotes the radial component of the diffusive mass flux of dissolved material (in kg/m²/s), averaged over an effective grain boundary thickness $\delta$, then the total rate of outflow across a cylindrical section of area $2\pi r\delta$ and the total rate of dissolution within this circumference must balance. Consequently,

$$2\delta J_r^{gb}(r) + \rho^s \llbracket v_n \rrbracket r = 0, \quad (3.1)$$

where $\rho^s$ denotes the density of the solid.

We shall further assume that the diffusive mass flux of the solute species, $J_r^{gb}$, is governed by Fick's law as applicable to bulk diffusion in a liquid phase. We introduce this relationship in the form $J_r^{gb} = -\rho^f L^{gb} d\mu/dr$, where $L^{gb}$ is a phenomenological mobility coefficient, $\mu$ is the mass-specific chemical potential of the solute component in the grain boundary solution phase, and $\rho^f$ is the bulk density of the latter. For a sufficiently dilute solution, Fick's law can be written
in terms of the gradient in the solute mass fraction C as (see, e.g., De Groot and Mazur 1962, Chapter 11)

$$J_r^{gb} = -\rho^f D_{gb} DC/dr. \qquad (3.2)$$

Here $D_{gb} = L^{gb}(d\mu/dC)_{p,T} \approx L^{gb}kT/\rho^s\Omega^s C$ is an approximately constant mass diffusivity (in m²/s) that controls solute diffusion through the intergranular fluid phase, where k is the Boltzmann constant and $\Omega^s$ denotes the molecular volume of the solid. This relationship follows from the standard expression

$$\mu = \mu^*(p, T) + (kT/\rho^s\Omega^s) \ln \gamma C \qquad (3.3)$$

for the chemical potential (in J/kg) of the solute component of a binary solute-diluent solution, when the concentration dependence of the activity coefficient $\gamma$ can be disregarded.

We may now combine equations (3.1) and (3.2) and integrate the resulting differential equation in $C(r)$ to obtain the radial distribution of the solute mass fraction along the intergranular contact:

$$C(r) = C(a) - \frac{\rho^s \llbracket v_n \rrbracket}{4\rho^f \delta D_{gb}} (a^2 - r^2). \qquad (3.4)$$

As expected for grain convergence ($\llbracket v_n \rrbracket < 0$), the solute concentration and therefore also its chemical potential are higher inside the grain contact than at its rim, in the interconnected pore space. By taking a surface average of this distribution over the grain contact, one obtains the following expression for the rate of grain convergence:

$$\llbracket v_n \rrbracket = - \frac{8\rho^f \delta D_{gb}}{\rho^s a^2} [\bar{C} - C(a)] \qquad (3.5)$$

in terms of the average concentration $\bar{C}$ along the grain boundary and the concentration $C(a)$ at the rim of the contact, i.e., in the pore space where $\mu$ and hence C are taken to be uniform.

To determine the unknown average concentration $\bar{C}$ in equation (3.5), we must now enter into a discussion of the chemical thermodynamics of stress-enhanced dissolution.

### 3.2.2 Chemical Equilibrium Between a Stressed Solid and Its Aqueous Solution

It was Sorby who in 1863 remarked that an effect similar in principle to the lowering of the freezing temperature of water with increasing pressure, first established experimentally by W. Thomson (Lord Kelvin) in 1850, must hold true with respect to the solubility of salts in water and might thus explain pressure solution phenomena such as pitted pebbles; and indeed in 1861 J. Thomson had already
demonstrated the stress-enhanced solubility of salt crystals in water experimen-
tally. But it remained for Gibbs (1961) to furnish a full theoretical explanation
of the effect of stress on the chemical equilibrium of an arbitrarily stressed solid
with its own solution phase. In particular, Gibbs showed that the condition for
local chemical equilibrium at a flat interface between a stressed pure solid and its
solution is given by (Gibbs 1961, Eq 387)<fnref n="3" />

$$
u^s - Ts^s + p/\rho^s = \mu, \tag{3.6}
$$

where $u^s$ and $s^s$ are the specific internal energy and specific entropy, respectively, of the pure solid phase, $p$ is the pressure in the solution phase, and $\mu$ is the chemical potential of the solute component in the solution phase (in the present context always a binary aqueous solution) as given by equation (3.3). Because $-p$ must equal the normal component of stress $\sigma_n$ at the solid side of the interface, the left side of condition (3.6) may be expressed solely in terms of quantities associated with the solid phase. Choosing to replace the first two terms by the specific Helmholtz free energy $f^s = u^s - Ts^s$, we may thus write equivalently

$$
f^s - \sigma_n / \rho^s = \mu. \tag{3.7}
$$

Condition (3.6) characterizes a strictly local state of equilibrium, pertaining
only to points along a solid-solution interface. For curved interfaces an interfacial
energy term of the form $(1/R_1 + 1/R_2)\gamma^{sf}/\rho^s$ must be added to its left side,
in which $R_1$ and $R_2$ denote finite principal radii of interfacial curvature (positive
when their centers lie on the side of the solid) and $\gamma^{sf}$ denotes the interfacial
energy (see Gibbs 1961, Eq 661; see also Heidug 1991).

One can picture a variety of situations involving a nonhydrostatically stressed
solid in *global* equilibrium with its own solution. In each case, however, the cri-
terion of global interfacial equilibrium will amount to nothing less than the re-
quirement that condition (3.6) be satisfied pointwise along the entire solid-fluid
interface in the system under consideration. It is only in the special case in which
a homogeneous solid is completely surrounded by a solution phase at uniform
chemical potential $\mu$ and pressure $p$ that the term $f^s + p/\rho^s$ becomes a constant
independent of position, thus appearing in the role of a specific Gibbs free energy
characterizing the bulk of the solid phase—on equal footing with the chemical po-
tential $\mu$ of the dissolved solid. In general, however, there exists no useful concept
of a scalar Gibbs free energy or chemical potential that could serve to define the
heterogeneous equilibrium of a solid composed of different phases or of a system
comprising a stressed solid in contact with a liquid solution. In fact, it has become
clear that the formal development of Gibbs’ thermodynamics of heterogeneous sys-
tems could be clarified in this point and in a sense made more complete through

<footnote n="3"><sup>3</sup>See also Lehner and Bataille (1984/85) for an alternative derivation of this result based on a defi-
nition of interfacial equilibrium as corresponding to a state of vanishing interfacial entropy production.</footnote>
the explicit recognition of the concept of a “chemical-potential tensor” (Bowen and Wiese, 1976; Grinfeld, 1982, 1991; Truskinovskiy, 1984; Heidug and Lehner, 1985). It will be noticed, indeed, that the condition of interfacial equilibrium shown in equation (3.7) can be written in the form

$$ \mu_n^s = \mu \qquad (3.8) $$

in terms of the normal component $\mu_n^s = n \cdot \mu^s \cdot n$ at the solid-solution interface (with unit normal $n$) of the tensor

$$ \mu^s = f^s \mathbf{1} - \sigma / \rho^s, \qquad (3.9) $$

where $\mathbf{1}$ is an isotropic unit tensor and $\sigma$ is the Cauchy stress. Under hydrostatic conditions, when $\sigma = -p\mathbf{1}$, $\mu^s$ is also isotropic and its orientation-independent normal component becomes identical with the specific Gibbs free energy $f^s + p/\rho^s$ or scalar chemical potential of a single-component solid.

### 3.2.3 A Phenomenological Rate Law for IPS

In the following we wish to derive a model for IPS that is free of any a priori bias toward one or the other rate-limiting process. We must therefore avoid the frequently made assumption that the intergranular solution phase is in chemical equilibrium with the stressed solid phase. Similarly, we shall admit a state of disequilibrium between the two phases along the hydrostatically loaded pore walls and allow for precipitation from a locally supersaturated pore fluid at such sites, or possibly free-face dissolution in an undersaturated pore fluid. We shall follow earlier work by Lehner and Bataille (1984/85) and Lehner (1990, 1995, 1997) and employ thermodynamic arguments that allow us to set up a phenomenological description of IPS.

Consider first the solid-fluid interface comprising the pore walls or free faces of the grains. We have already seen that Gibbs' condition, shown in equation (3.8), provides a local condition of equilibrium along this boundary. In the interpretation of the thermodynamic theory of irreversible processes, violation of equation (3.8) implies the existence of a thermodynamic force $\mu_n^s - \mu = 0$ that will drive a flux of the soluble solid substance across the phase boundary, implying either dissolution or precipitation.

Let $w_n$ denote the local speed of displacement of the solution surface in the direction of an outwardly directed unit normal and let $v_n$ be the normal component of the solid-particle velocity at the interface. The difference $w_n - v_n$ therefore designates the velocity of the interface relative to the material, while the product $r_{sf} = \rho^s(w_n - v_n)$ yields the flux (in kg/m²/s) of solute mass per unit area across the pore wall ($r_{sf} < 0$ for solution). The product $r_{sf}(\mu_n^s - \mu)$ can be shown to equal the rate at which free energy is dissipated per unit area as the material passes through the phase boundary (see Lehner and Bataille 1984/85). This rate
is expressible as the product of the absolute temperature T (which is assumed to remain continuous across the phase boundary) and a nonnegative interfacial entropy production rate $\dot{\sigma}$. Accordingly, one has

$$-r_{sf}(\mu_n^s - \mu) = T\dot{\sigma} \ge 0, \qquad (3.10)$$

which must hold for any value of the flux. It follows that the work-conjugate thermodynamic force (or “affinity”) $\mu_n^s - \mu$ and flux $r_{sf}$ must be functionally dependent. The simplest phenomenological relation of this kind is obtained upon truncating a Taylor expansion of the function $r_{sf} = f(\mu_n^s - \mu)$ after the linear term, i.e., putting $r_{sf} \approx -L(\mu_n^s - \mu)$. This may be written in the form

$$r_{sf} = -\rho^s K^{sf} (\rho^s \Omega^s / kT) (\mu_n^s - \mu), \qquad (3.11)$$

where the temperature-dependent, nonnegative phenomenological coefficient $L$ is replaced by the more convenient rate “constant” $K^{sf}$ with the dimension of a velocity. These coefficients assume different values, $L_+$ and $L_-$ or $K_+^{sf}$ and $K_-^{sf}$, for solution and precipitation, respectively.

A phenomenological relation of the type of equation (3.11), which governs the propagation of a phase boundary, was apparently employed first by Machlin (1953) in an analysis of the growth of precipitates in metals. The advantage of such a thermodynamic description, apart from its simplicity, lies in the clarification of the nature of the driving force (i.e., one that involves the normal component of a chemical potential tensor). This can be of great conceptual help in the development of rational theories of deforming two-phase materials. The fact remains, however, that the linear dependence in equation (3.11) of the flux $r_{sf}$ on the driving force $\mu_n^s - \mu$ is indeed only a consequence of the linearization of an appropriate nonlinear kinetic model relation (see Dibble and Tiller 1981). One must therefore restrict the validity of (3.11) to small driving forces, implying $\rho^s \Omega^s |\mu_n^s - \mu| / kT \ll 1$ in the present case. Relation (3.11) was proposed by Lehner and Bataille (1984/85) as a means to account for the effects of nonhydrostatic stresses in a dissolving solid. Including the surface energy term (Heidug, 1991), it has since been employed in a substantial number of theoretical modeling studies by Heidug and Leroy (1994), Leroy and Heidug (1994), and Ghoussoub and Leroy (2001) concerned with the stability and evolution of solid–fluid interfaces under conditions of pressure solution and redeposition.

Let us also observe that the thermodynamic force-flux relation equation (3.11) may be brought into the more familiar form of a first-order kinetic relation by use of equations (3.3) and (3.6), which allow us to express the solid chemical potential $\mu_n^s$ in terms of an equilibrium concentration $C_\sigma$, i.e., the solute concentration required to maintain chemical equilibrium between the stressed solid and the fluid phase locally along the solid–fluid interface. For dilute concentrations this gives

$$r_{sf} = \rho^s K^{sf} (C - C_\sigma) / C_\sigma, \qquad (3.12)$$

$C$ being the actual solute concentration in the pore fluid. We shall consider this type of kinetic relation further later in the discussion (see equation [3.26]).
Turning now to the kinetics of intergranular dissolution and diffusive transport, let us denote by $r_{gb} = \rho^s \llbracket v_n \rrbracket = \rho^s \dot{l}_n d$ the uniform dissolution rate along a given intergranular contact, i.e., the flux of solid material (in kg/m²/s) from both grains into the contact zone and across the actual solid-solution phase boundary lying within that zone (see Figure 3.3). Also, let $\mu_n^s(r)$ denote the normal component of the solid-phase chemical-potential tensor along the smooth outer margins of the contact zone of Figure 3.3—defined in terms of the values of $\sigma_n$, $\rho^s$, and $f^s$ along these margins—and $\mu(r)$ denote the chemical potential of the dissolved solid component in the grain boundary fluid. The product $-r_{gb}[\mu_n^s(r) - \mu(r)]$ then equals the rate at which work is dissipated per unit nominal contact area as the material passes through the solid contact zone and the phase boundary.

Furthermore, the work dissipated per unit area in the isothermal diffusive transfer of material, driven by the potential drop $\mu(r) - \mu(a)$ at the rate $r_{gb}$ from its solution site at $r$ to the contact periphery at $r = a$, is equal to $-r_{gb}[\mu(r) - \mu(a)]$. It follows that the total rate of dissipation associated with irreversible processes in an intergranular contact zone, a quantity that we denote by $\dot{\Delta}_{gb}$, must equal the sum of the integrated dissipative work rates for this zone, i.e.,

$$ \dot{\Delta}_{gb} = -r_{gb} \int_{S_{gb}} [\mu_n^s(r) - \mu(r)]dA - r_{gb} \int_{S_{gb}} [\mu(r) - \mu(a)]dA, \quad (3.13) $$

where the integrals are taken over the nominally flat grain-to-grain contact $S_{gb}$ of Figure 3.3. Note that on carrying out a formal integration over $S_{gb}$ and canceling terms, this becomes

$$ \dot{\Delta}_{gb} = -\pi a^2 r_{gb} [\bar{\mu}_n^s - \mu(a)], \quad (3.14) $$

from which it is apparent that the potential difference $\bar{\mu}_n^s - \mu(a)$ represents the average total driving force for the combined (serial) processes of intergranular dissolution and grain boundary diffusion.

If the first integral in equation (3.13) can be neglected against the second, grain boundary diffusion acts as the rate-limiting process and the driving chemical potential difference associated with the fast dissolution step may equivalently be set equal to zero, i.e., $\mu_n^s \approx \mu$. It is this Stefan-type approximation that allows equation (3.7) to be treated as an equilibrium condition as in the work of Paterson (1973).

In general, however, one cannot disregard a priori the dissipation associated with the dissolution of the solid phase within the contact zone, so that the first term in equation (3.13) must be retained. As in the above, we may now invoke the second law to assert that this product must remain nonnegative at any value of the driving force and to conclude that the force and flux must therefore be functionally related. In the linear approximation, this yields a force-flux relation of the form

$$ r_{gb} = -\rho^s K^{gb} (\rho^s \Omega^s / kT) [\mu_n^s(r) - \mu(r)], \quad (3.15) $$

where $K^{gb} > 0$ is phenomenological rate coefficient.
The use of relation (3.15) can be justified, as long as the solid material may be assumed to deform elastically while advancing through the contact zone toward the solution surface, except for the immediate vicinity of load-bearing islands. These islands are imagined to be undercut by dissolution but to deform by dissipative inelastic micromechanisms only in the final phase of this undercutting process, when the load is about to be transferred from such a deforming contact point onto neighboring stiffer contacts. Dissipation in the solid portion of the contact zone is therefore neglected against the dissipation associated with dissolution. While the local undercutting and removal of islands by dissolution is essential to achieving grain convergence by IPS, inelastic solid deformation can be viewed as an associated phenomenon, made possible only by the former and progressing therefore at the same rate. In the following we shall apply equation (3.15) without any correction for work that may be dissipated in deforming the solid material within the contact zone, and we shall treat $K^{gb}$ as a phenomenological coefficient governing the dissolution rate at a stressed but fluid-infiltrated grain boundary that possesses a dynamically stable island-channel structure (Raj, 1982; Lehner and Bataille, 1984/85; Lehner, 1990; Spiers and Schutjens, 1990; Cox and Paterson, 1991).

The second integral in equation (3.13) represents the rate of dissipation associated with solute diffusion along the grain boundary. The relevant kinetic relation is Fick's law (equation [3.2]), which has already been used in equation (3.5) to obtain the rate of grain convergence in terms of the average concentration $\bar{C}$ along the nominally flat circular contact $S_{gb}$. Let us first rewrite $\bar{C}$ in terms of the averaged chemical potential along $S_{gb}$. Under the assumption of dilute solute concentrations, it follows from equation (3.3) that $[C(r) - C(a)]/C(a) \approx \ln C(r)/C(a) \approx (\rho^s \Omega^s/kT)[\mu(r) - \mu(a)]$, and $[\bar{C} - C(a)]/C(a) \approx (\rho^s \Omega^s/kT)[\bar{\mu} - \mu(a)]$, so that after averaging over $S_{gb}$, equation (3.5) may be written

$$r_{gb} = -(8\rho^f C_{eq} \delta D_{gb} / a^2) (\rho^s \Omega^s / kT) [\bar{\mu} - \mu(a)], \quad (3.16)$$

where $C_{eq}$, the equilibrium concentration (mass fraction), has here been introduced as an approximation of the actual concentration $C(a)$ in the open pore space. This approximation can be justified for the contemplated geological applications by the uncertainties in the numerical values of other parameters, such as the grain boundary diffusion factor $\delta D_{gb}$.

Similarly, one has

$$r_{gb} = -\rho^s K^{gb} (\rho^s \Omega^s / kT) [\bar{\mu}_n^s - \bar{\mu}] \quad (3.17)$$

for the surface-averaged form of equation (3.15). Elimination of the chemical potential $\bar{\mu}$ from equations (3.16) and (3.17) now produces the relation

$$[v_n] = r_{gb} / \rho^s = -d(\rho^s \Omega^s / kT) [\bar{\mu}_n^s - \mu(a)] / [\tau_S + (4a^2/d^2)\tau_D] \quad (3.18)$$

for the uniform rate of grain convergence at the circular contact between two identical spherical grains. Here it has been found appropriate to introduce the
initial grain size *d* as a reference length, yielding the following definition of the characteristic times:

$$ \tau_s = d/K^{gb} \quad \text{and} \quad \tau_D = \rho^s d^3/(32\rho^f C_{\text{eq}}\delta D_{gb}), \qquad (3.19) $$

associated with the rates of intergranular solution and grain boundary diffusion, respectively.

In essence, the above result already represents a “creep law” for our simple model system. We note, in particular, that the difference $\bar{\mu}_n^s - \mu(a)$ between the average solid-phase chemical potential in the grain-to-grain contact and the (uniform) solute chemical potential in the pore space appears in the role of the thermodynamic force that drives IPS—a result that is consistent with expression (3.14) for dissipation associated with this process. We must now seek an expression for this driving force, by use of equation (3.9), in terms of the intergranular contact stress. For this purpose, we shall find it convenient to introduce as a reference chemical potential the potential $\mu_0 = \mu^*(p, T) + (kT/\rho^s\Omega^s) \ln \gamma C_{\text{eq}}$ of a solution that is in equilibrium with the solid phase, when both phases are subjected to the currently prevailing fluid pressure and temperature, *p* and *T*, respectively. Along a flat solution surface, the above value of $\mu_0$ therefore defines a reference equilibrium value for the normal component of the solid-phase chemical potential tensor, in agreement with Gibbs’ condition (3.8), equal to $\mu_n^s = f_0^s + p/\rho_0^s = \mu_0$. In terms of this reference value, the local potential difference $\mu_n^s(r) - \mu(a)$ along an intergranular contact may be written as the following sum:

$$ \mu_n^s(r) - \mu(a) = [\mu_n^s(r) - \mu_0] - [\mu(a) - \mu_0] = (f^s - f_0^s) + (1/\rho^s - 1/\rho_0^s)p - (\sigma_n + p)/\rho^s - [\mu(a) - \mu_0]. \quad (3.20) $$

Before making use of this expression in equation (3.18), we introduce the approximation

$$ \mu_n^s - \mu_0 \approx -(\sigma_n + p)/\rho^s, \qquad (3.21) $$

implying the neglect of the first two terms on the right side of equation (3.20) against the third term. This approximation is frequently used in combination with the Kamb-Paterson equilibrium assumption $\mu_n^s \approx \mu$ (Paterson, 1973; Rutter, 1976; Lehner and Bataille, 1984/85). Its validity rests on the observation that the ratio of the third over the first two terms in equation (3.20) is of the order of $|\sigma_n + p|/E$ for elastic grains with a Young modulus *E*, and that this ratio remains typically much smaller than 1 in the situations of interest. Accepting this approximation yields the following expression for the local driving force in equation (3.15):

$$ \mu_n^s(r) - \mu(r) \approx -(\sigma_n + p)/\rho^s - [\mu(r) - \mu_0], \qquad (3.22) $$

whereas the overall average driving force in equation (3.18) becomes

$$ \bar{\mu}_n^s - \mu(a) \approx -(\bar{\sigma}_n + p)/\rho^s - [\mu(a) - \mu_0]. \qquad (3.23) $$
We therefore have an expression for the rate of grain convergence in which the thermodynamic driving force is expressed in terms of an effective intergranular normal stress and the change in solute chemical potential from the hydrostatic equilibrium level $\mu_0$:

$$ \llbracket v_n \rrbracket = d(\rho^s \Omega^s/kT)[(\bar{\sigma}_n + p)/\rho^s + \mu(a) - \mu_0]/[\tau_s + (4a^2/d^2)\tau_D], \quad (3.24) $$

or alternatively, in terms of concentrations (mass fractions),

$$ \llbracket v_n \rrbracket = d[(\rho^s \Omega^s/kT)(\bar{\sigma}_n + p)/\rho^s + (C - C_{eq})/C_{eq}]/[\tau_s + (4a^2/d^2)\tau_D]. \quad (3.25) $$

To complete this discussion of the kinetics of pressure solution, we now return to relation (3.11) to note that the driving potential difference may again be replaced by $(\rho^s \Omega^s/kT)(\mu_n^s - \mu_0) - (C - C_{eq})/C_{eq}$, where the first term could be set equal to $(f^s - f_0^s) + (1/\rho^s - 1/\rho_0^s)p$ as in equation (3.20) for entirely elastic changes in the free energy of the solid. Although this would seem appropriate for a description of free-face pressure solution, precipitation may well occur by incoherent addition of a hydrostatically stressed solid phase, as has already been argued by Gibbs from considerations of stability. Within the framework of our approximate model, it therefore seems reasonable to neglect the contribution from nonhydrostatic stresses to the free-face solution or precipitation rate $r_{sf}$ and, instead of equation (3.12), use the first-order kinetic relation

$$ r_{sf} = \rho^s(w_n - v_n) = \rho^s K^{sf} (C - C_{eq})/C_{eq}. \quad (3.26) $$

## 3.3 Elements of a Continuum Theory of IPS

In this section we shall turn to the macroscopic continuum theory of porous materials that are assumed to compact by a process of IPS in the manner that was described in Section 3.2. This will require some attention to the construction of the pertinent macrovariables from the corresponding grain-scale fields, for example, when setting up mass balance equations with a source term that expresses an average rate of loss or gain of material by solution or precipitation processes. The same interpretation problem exists of course for the constitutive macrovariables of stress and deformation that we wish to employ. To expose and elucidate the structure of a certain type of two-scale theory forms the principal aim of this chapter, but for this to be accomplished we find it necessary to simplify not only the chemistry and kinetics, but also the geometry at the grain scale. Our goal has been the development of a theory of creep coupled to macroscopic solute transport in an *open system* up to a point at which it would yield predictions for the behavior of simple macroscopic systems, for example, the one-dimensional compaction of a layer of sandstone by IPS or of a sandstone sample in a laboratory experiment. We think it important to have a relatively simple description of this type available
at this stage as a testing ground for further theoretical developments as well as for
guiding the design of future laboratory experiments.

**3.3.1 Volume Averaging and Kinematics of IPS**

The approximate theory of solution and diffusion processes at a wet grain boundary that has been discussed in the previous section is essentially a grain-scale continuum description; it also defines the microscale considered in this chapter, whereas the macroscale is the scale at which we wish to formulate and solve boundary value problems for porous sedimentary layers that compact by a pressure solution creep process. It would be difficult now to give a macroscale description of the aggregate deformation by such processes without some attention to an averaging procedure through which the variables of the macroscopic theory must be related to the microscale picture of the previous section. We shall therefore begin by recollecting a few facts about spatial averaging, from which we shall be able to proceed with a discussion of the requisite kinematic variables of a macroscopic theory of creep by IPS.

A microscale description usually deals with processes that take place within an REV of some porous medium (in the present case typically a porous sandstone). An REV occupies a closed region $V$ in space and comprises a large enough mass of solid skeleton material such that the relevant field variables defined on $V$ will possess meaningful macroscopic (spatial) averages. With reference to Figure 3.4, let $V_f$ denote the interconnected, fluid-saturated pore space lying within the boundary $\partial V$ of $V$, and $V_s$ the portion occupied by the solid skeleton. The latter would, in fact, include any isolated cavities or fluid inclusions that do not form part of the interconnected porosity. For convenience, we use the same symbols to denote a region and its volume. As shown in Figure 3.4, the boundary $\partial V$ comprises the fluid–fluid and solid–solid intersections $S_{ff}$ and $S_{ss}$, respectively, of the imagined cut. Within $V$, one distinguishes the solid–fluid phase boundary (pore walls), $S_{sf}$ and the grain boundaries $S_{gb}$.

The volume fractions $\phi_v$ ($v = s, f$) of the solid or void-filling fluid phase are defined in the usual manner by $\phi_v = V_v/V$ ($v = s, f$), it being understood that the volumes $V$ and $V_v$ are determined by suitable techniques on a rock sample. Because $V_s$ and $V_f$ add up to $V$, it is clear that the volume fractions of the two phases satisfy the relation $\phi_s + \phi_f = 1$, and we shall therefore eventually revert to the standard notation $\phi = \phi_f$ for the porosity and $1 - \phi$ for the volume fraction of the solid phase.

We note here that the existence of a useful macroscopic point function with the significance of a volume fraction $\phi_v$ depends on the existence of a size range of the REV such that the value of $\phi_v$ remains stationary within that range. This property of stationarity of any field variable is typically taken for granted in continuum theories of porous media. Certain macroscopic variables of these theories, such as a partial density of an extensive property or an average particle velocity of the solid or fluid phase may thus be interpreted as one or the other of the
<figure><img src="image_5.png" /><figcaption>**Figure 3.4** ► REV of a porous rock.</figcaption></figure>

following simple volume averages of some microscale field variable $\psi$, i.e., its phase average:

$$ \psi_v = \frac{1}{V} \int_{V_v} \psi \, dv, \quad v = s, f, \qquad (3.27) $$

or its intrinsic phase average:

$$ \psi^v = \frac{1}{V_v} \int_{V_v} \psi \, dv, \quad v = s, f. \qquad (3.28) $$

From these definitions, it is apparent that $\psi_v$ represents a *partial* quantity that is related to $\psi^v$ by

$$ \psi_v = \phi_v \psi^v, \quad v = s, f. \qquad (3.29) $$

When a phase average exists for each of the two phases $s$ and $f$, then the sum $\psi = \psi_f + \psi_s$ defines a bulk property of the porous medium, for example the bulk mass density $\rho = \rho_f + \rho_s = \phi_f \rho^f + \phi_s \rho^s$, which is usually defined in terms of the actual intrinsic mass densities $\rho^f$ and $\rho^s$ and the respective volume fraction as in the second expression. Note that this notation does not distinguish between a local variable and its average, because the scale to which a variable applies will
usually be clear from the context<fnref n="4" /> so that no ambiguity should arise. Only where the distinction is to be emphasized shall we employ the notation $\langle\psi\rangle$ for a spatial average.

Suppose now that the quantity $\psi$ represents a component, in a given Cartesian coordinate system, of the particle velocity of the solid phase. Its intrinsic phase average is then given by

$$v_i^s = \frac{1}{V_s} \int_{V_s} v_i dv. \qquad (3.30)$$

In terms of this average velocity of the solid skeleton, one may evidently define as usual the macroscopic rate of deformation as the symmetric part of the velocity gradient $\partial v_i^s / \partial x_j$ by

$$D_{ij} = \frac{1}{2} \left( \frac{\partial v_i^s}{\partial x_j} + \frac{\partial v_j^s}{\partial x_i} \right), \qquad (3.31)$$

where we omit a superscript on $D_{ij}$. But to see how $D_{ij}$ may actually be computed in terms of the microscale motions that occur within an REV, we shall now make use of the following result<fnref n="5" /> for the derivative of a volume integral of physical variable $\psi$ (scalar, vector, or tensor component) over the region $V_s$ (the same result holds of course for $V_f$):

$$\frac{\partial}{\partial x_j} \int_{V_s} \psi dv = \int_{S_{ss}} \psi n_j da. \qquad (3.32)$$

Here $n_j$ denotes a component of the outer unit normal to the surface $S_{ss}$. Division by $V_s$ and an application of the divergence theorem to the term on the right yields a result that has become known as the averaging rule or even averaging theorem in the literature on volume-averaged transport equations for two-phase media (Slattery, 1972):

$$\frac{\partial \psi^s}{\partial x_j} = \frac{1}{V_s} \int_{V_s} \frac{\partial \psi}{\partial x_j} dv - \frac{1}{V_s} \int_{S_{sf}} \psi n_j da + \frac{1}{V_s} \int_{S_{gb}} [[\psi]] n_j da, \qquad (3.33)$$

where the unit normal vector along the phase boundary $S_{sf}$ points into the fluid phase. Along the grain boundary $S_{gb}$ the variable $\psi$ is allowed to experience a

<footnote n="4">⁴As in the case of the sum $\rho = \rho_f + \rho_s$, for example (see p. 130).</footnote>

<footnote n="5">⁵The result may be deduced without difficulty by inverting the steps taken in the proof of the divergence theorem (see O.D. Kellogg's *Foundations of Potential Theory*, Dover Publs., 1953, Sec. VIII.7), generalizing Leibniz's rule for the differentiation of an integral with respect to a parameter. It may also be deduced from Reynolds' transport theorem. In the form of equation (3.33) the result is familiar from the theory of distributions but may again be obtained directly, if not from equation (3.32), then through elementary arguments (Slattery, 1972).</footnote>
jump discontinuity $\llbracket \psi \rrbracket(\mathbf{x}, t) = \psi^+(\mathbf{x}, t) - \psi^-(\mathbf{x}, t)$ from its value in the grain on the negative side to its value in the neighboring grain on the positive side of $S_{gb}$, the surface normal defining the (arbitrarily fixed) positive side (see Figure 3.4).

Application of equations (3.32) and (3.33) to equation (3.31) now yields

$$ 
\begin{align} 
D_{ij} &= \frac{1}{2V_s} \int_{S_{ss}} (v_i n_j + v_j n_i) da \\ 
&= \frac{1}{2V_s} \int_{V_s} \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right) dv - \frac{1}{2V_s} \int_{S_{sf}} (v_i n_j + v_j n_i) da \\ 
&\quad + \frac{1}{2V_s} \int_{S_{gb}} (\llbracket v_i \rrbracket n_j + \llbracket v_j \rrbracket n_i) da. \tag{3.34}
\end{align} 
$$

This shows that during active IPS there may be three distinct contributions to the macroscopic rate of deformation. In the applications considered in this chapter, the first contribution, which results from the internal deformation of the grains, will typically be neglected against the last one, which arises from IPS. To appreciate the significance of the second term in equation (3.34), one may think of two extreme situations. Suppose, first, that the grains are approximately undeformable, so that the bulk aggregate deformation results almost entirely from the velocity discontinuity along grain-to-grain contacts; the second term is then determined by the third and the grain or pore shape. For a spatially periodic arrangement of rigid grains that are truncated by intergranular dissolution, one may in fact choose a single grain to define an REV (as we have done in the previous section) and consider the grain as being fixed in space. In this case the third term in equation (3.34) must therefore account entirely for the bulk deformation rate. If, on the other hand, one were to study an insoluble and incompressible, perhaps rubber-like solid skeleton material, then the bulk volume deformation would obviously result from the displacement of the pore walls and would thus be accounted for by the second term in equation (3.34).

Let us now recall a few basic definitions and results from continuum mechanics, which we shall require subsequently (see, e.g., Chadwick 1999). It is assumed here that the motion of a given *body* of porous rock can be described, in the usual manner, in terms of that of its *particles* or *material points*, by specifying the position $\mathbf{x}$ of any given particle (relative to a chosen origin) in the current configuration of the body as a function of its position $X$ in some reference configuration and of time $t$, writing<fnref n="6" />

$$ \mathbf{x} = \mathbf{x}(X, t). \tag{3.35} $$

One may think here of a *material point* as corresponding to the center of mass of a small portion of the solid skeleton, which remains in the solid state throughout

<footnote n="6">⁶For simplicity, a function and its value are often denoted by the same symbol in this chapter.</footnote>
the entire motion; for example, a suitable subset of the solid part of an REV in the reference configuration, the existence of which we shall stipulate. An alternative description of the motion of a body is furnished by the inverse function $X = X(x, t)$ of equation (3.35). If the scalar, vector, and tensor fields that define the relevant properties of a deforming body are defined on the reference configuration as functions of $X$ and $t$, they are said to be given in the reference description, while if expressed as functions of $x$ and $t$, the fields are given in the spatial description. For example, the velocity of a particle is defined by the partial derivative of its motion $x(X, t)$ with respect to time, holding $X$ fixed. Through use of the inverse motion, it may, however, be represented by the fields $\dot{x}(X, t)$ and $v(x, t)$ according to

$$ \dot{\mathbf{x}} = \frac{\partial \mathbf{x}}{\partial t}(X, t) = \frac{\partial x}{\partial t}[X(x, t), t] \equiv v(x, t) \quad (3.36) $$

in the referential and the spatial description, respectively. In this chapter, a superimposed dot always denotes the so-called *material derivative* with respect to time, i.e., the rate of change of a quantity associated with a fixed particle. When a field is given in the spatial description, its material derivative is therefore obtained through substitution of the motion and application of the chain rule. For the scalar field $\psi(\mathbf{x}, t)$, for example,

$$ \begin{aligned} \dot{\psi}(\mathbf{x}, t) &= \frac{\partial}{\partial t} \psi[x(X, t), t]_{X} \\ &= \frac{\partial \psi}{\partial t}(\mathbf{x}, t) + \frac{\partial \psi}{\partial x_i} \frac{\partial x_i}{\partial t} [X(\mathbf{x}, t), t] \\ &= \frac{\partial \psi}{\partial t}(\mathbf{x}, t) + \{(\mathrm{grad}\psi)(\mathbf{x}, t)\} \cdot v(\mathbf{x}, t). \end{aligned} \quad (3.37) $$

When studying the motion of particles in a neighborhood of a given particle $X$, the local deformation of a body is brought into view through the gradient of the motion (3.35), that is, the deformation gradient $F = \mathrm{Grad}x(X, t)$ (or $F_{i\alpha} = \partial x_i / \partial X_{\alpha}$ in Cartesian component notation). The motion $x = x(X, t)$ is usually assumed to meet the requirement $J \equiv \det F > 0$ for its Jacobian determinant at all $X$. Here the sign ensures the transformation of a positive volume element $dV$ with position $X$ in the reference configuration into a positive volume element $dv = JdV$ at $\mathbf{x}$ in the current configuration. Although it is natural as a rule to stipulate the condition of noninterpenetrability of matter, which is implied by the fixed sign of $J$, such a requirement may well be unwarranted for soluble porous solids. For now it can happen in principle that a portion of the original solid skeleton is completely removed by dissolution with or without replacement by newly precipitated material. Although the kinematics of soluble porous materials could be developed in sufficient generality to allow for such situations, it is assumed here that they will not occur, the implication being that the solid skeleton of the rock under study will always preserve enough of its initial solid substance to allow a meaningful definition of the motion (3.35) in terms of a single reference
configuration. It is clear then that the velocity defined by equation (3.36) in terms of this motion can be identified, in a spatial description, with the average velocity $v^s$ of the solid phase as given by equation (3.30). The trace $D_{ii} = \text{div}v^s$ of the rate of deformation is therefore also linked to the relative rate of change in bulk volume $\dot{J} = (dv) \cdot /dV$ through Euler's relation

$$ \dot{J} = J \text{div}v^s. \qquad (3.38) $$

When dealing with large deformations, a referential description will often be found more concise and transparent than its spatial counterpart, and this is true also for the present problem, where we shall find it convenient to discuss the change of pore volume in terms of the pore volume per unit referential volume, or *pore volume fraction* $v_f = \phi J$, and in terms of $v_s = (1-\phi)J$, the solid volume per unit referential volume. Because the rates of change of these quantities are related by

$$ \dot{v}_f + \dot{v}_s = \dot{J}, \qquad (3.39) $$

$\dot{v}_f$ may be determined from $\dot{J}$, via Euler's relation (3.38) and the trace of the deformation rate (3.34) and from $\dot{v}_s$. An expression for the latter is obtained, by taking the material derivative of its defining expression, as follows:

$$ \begin{aligned} \dot{v}_s &= \frac{d}{dt} \left\{ \frac{1}{V_0} \int_{V_s} dv \right\} \\ &= \frac{1}{V_0} \int_{\partial V_s} w_n da = \frac{1}{V_0} \int_{S_{sf}} w_n da + \frac{1}{V_0} \int_{S_{ss}} v_n da \\ &= \frac{1}{V_0} \int_{S_{sf}} (w_n - v_n) da + \frac{1}{V_0} \int_{S_{gb}} \llbracket v_n \rrbracket da + \frac{1}{V_0} \int_{V_s} \text{div}v dv. \end{aligned} \qquad (3.40) $$

Here we have applied a form of Reynolds' transport theorem (see, e.g., Chadwick 1999, Section 3.6) and the divergence theorem, subtracting and adding a term $V_0^{-1} \int_{S_{sf}} v_n da$. In this way the speed of propagation, $w_n - v_n$, of the grain surface $S_{sf}$ relative to the material, rather than its speed of displacement, $w_n$, in a given coordinate system, enters into the surface integral over $S_{sf}$ as the kinetically relevant rate.

We note further that, in accordance with equations (3.37) and (3.38), the material derivative of any quantity $\psi J$ can be written

$$ \begin{aligned} (\psi J) \cdot &= J \dot{\psi} + \psi \dot{J} = J(\partial_t \psi + v^s \cdot \text{grad}\psi) + \psi J \text{div}v^s \\ &= J\{\partial_t \psi + \text{div}(\psi v^s)\}. \end{aligned} \qquad (3.41) $$
### 3.3.2 Balance Equations for Mass

To write a mass balance equation for each of the two phases $v = f, s$, let us assume now that the mass densities $\rho^v$ and velocities $v^v$ are well-defined macroscopic point functions of $x$ that may be interpreted as intrinsic phase averages in the above sense. The local mass balances for the two phases then are

$$ \partial_t (\phi_s \rho^s) + \mathrm{div}(\phi_s \rho^s v^s) = r_s \quad (\text{for the solid phase}) \quad (3.42) $$

and

$$ \partial_t (\phi_f \rho^f) + \mathrm{div}(\phi_f \rho^f v^f) = -r_s \quad (\text{for the fluid phase}), \quad (3.43) $$

where $\partial_t$ denotes partial differentiation with respect to time at a fixed spatial position $x$. The source term $r_s$ represents the rate of change of the solid skeleton mass per unit bulk volume by precipitation ($r_s > 0$) or solution ($r_s < 0$) processes. It has the following interpretation as an average over grain-scale interfacial mass transfer rates:

$$ r_s = \frac{1}{V} \int_{S_{sf}} \{\rho(w_n - v_n)\} da - \frac{1}{V} \int_{S_{gb}} [[\rho(w_n - v_n)]] da. \quad (3.44) $$

Mass conservation requires the specific mass flux to remain continuous across $S_{sf}$, so that it could in principle be determined in terms of the appropriate densities and particle velocities on either side of this phase boundary. However, unless required otherwise, it is best to think of the integration over $S_{sf}$ as being performed on its negative side, i.e., in the solid phase. Along the grain boundary, on the other hand, the mass flux across individual grain faces will experience a jump across $S_{gb}$ during active IPS, as has been discussed already for the grain-scale model of the previous section.

Using the notation $\phi \equiv \phi_f = 1 - \phi_s$ for the porosity, one can rewrite equation (3.43) in terms of the Darcy filter velocity

$$ q = \phi(v^f - v^s) \quad (3.45) $$

to obtain

$$ \partial_t (\phi \rho^f) + \mathrm{div}(\phi \rho^f v^s) + \mathrm{div}(\rho^f q) = -r_s. \quad (3.46) $$

The fluid phase forms a binary mixture of aqueous silica (the solute) and water (the solvent). For these we define the partial mass densities $\rho_1$ (for silica) and $\rho_2$ (for water), which add up to the (intrinsic) fluid density according to

$$ \rho^f = \sum_{k=1}^{2} \rho_k. \quad (3.47) $$
As concentration variables we shall employ the mass fractions

$$C_k = \frac{\rho_k}{\rho^f}, \quad k = 1, 2, \qquad (3.48)$$

which satisfy the condition $\sum C_k = 1$.

To each of the two components of the fluid mixture one may ascribe an average particle velocity $v^k$ such that the partial mass flow rates $\rho_k v^k$ associated with these particle velocities add up to the total mass flow rate of the fluid and

$$v^f = \sum_{k=1}^{2} C_k v^k \qquad (3.49)$$

therefore becomes the velocity of the local center of mass, or barycentric velocity of the fluid mixture.

A balance equation for the solute mass that accounts for the addition or loss by solution or precipitation processes must be of the form

$$\partial_t(\phi\rho_1) + \text{div}(\phi\rho_1 v^1) = -r_s, \qquad (3.50)$$

featuring the same source term as equation (3.46). The mass flux appearing in the second term of this equation may now be written in a physically meaningful way as a sum of three fluxes. Clearly, because we wish to account for molecular diffusion as a mode of solute transport, the motion of one species relative to the mixture as a whole must be described in terms of an appropriate, kinetically relevant diffusive flux. Here we shall consider the velocity of a particle species relative to the barycentric velocity and define the diffusive mass flux vector for this species by

$$\dot{\boldsymbol{j}}_k = \rho_k \phi (\boldsymbol{v}^k - \boldsymbol{v}^f) = C_k \rho^f \phi (\boldsymbol{v}^k - \boldsymbol{v}^f), \quad k = 1, 2, \qquad (3.51)$$

where $\sum \dot{\boldsymbol{j}}_k = \mathbf{0}$ follows from equation (3.49). There exists therefore only one independent diffusive flux $\dot{\boldsymbol{j}} \equiv \dot{\boldsymbol{j}}_1$, say, just as there is only one independent concentration $C \equiv C_1$. Making use of definition (3.51) one may thus express the solute mass flow rate in equation (3.50) in terms of $\dot{\boldsymbol{j}}$, $\boldsymbol{v}^s$, and the filter velocity $\boldsymbol{q}$, putting $\phi\rho^1 v^1 = \dot{\boldsymbol{j}} + C\rho^f\phi\boldsymbol{v}^f = C\rho^f\phi\boldsymbol{v}^s + C\rho^f\boldsymbol{q} + \dot{\boldsymbol{j}}$. Substitution of this sum for the mass flux then yields the balance equation for the solute mass

$$\partial_t(C\phi\rho^f) + \text{div}(C\phi\rho^f\boldsymbol{v}^s) + \text{div}(C\rho^f\boldsymbol{q}) + \text{div}\dot{\boldsymbol{j}} = -r_s. \qquad (3.52)$$

To cast the balance equations (3.42), (3.46), and (3.52) in a referential description, we now introduce the referential filter velocity $\boldsymbol{Q}$, defined by

$$\boldsymbol{Q} = JF^{-1}\boldsymbol{q} \qquad \boldsymbol{q} = J^{-1}F\boldsymbol{Q} \qquad (3.53)$$

and the referential diffusive mass flux vector $\boldsymbol{J}$, defined by

$$\boldsymbol{J} = JF^{-1}\boldsymbol{j} \qquad \boldsymbol{j} = J^{-1}F\boldsymbol{J}. \qquad (3.54)$$
Because an oriented material surface element nda in the current configuration becomes $NdA = J^{-1}F^T nda$ in the reference configuration ($F^T$ denotes the transpose of $F$), the definitions (3.53) and (3.54) are seen to satisfy the invariance requirement $q \cdot nda = Q \cdot NdA$ and $j \cdot nda = J \cdot NdA$, respectively, for the scalar fluxes through a material surface element (see Chadwick, 1999, Chapter 2, Equation 21). Thus, upon integrating these fluxes over the closed boundary of a material region in the two configurations, one must have $\int_{\partial\Omega_0} Q \cdot NdA = \int_{\partial\Omega} q \cdot nda$ or, by application of the divergence theorem and a subsequent change from the spatial to the referential description, $\int_{\Omega_0} \text{Div} Q dV = \int_{\Omega} \text{div} q dv = \int_{\Omega_0} J \text{div} q dV$. From this follows $\text{Div} Q = J \text{div} q$ or $\partial Q_\alpha / \partial X_\alpha = J \partial q_i / \partial x_i$ (in referential and spatial Cartesian coordinates $X_\alpha$ and $x_i$, respectively) and the analogous result for the vector fields $J$ and $j$.

In addition to the referential flux vectors, we define the referential partial mass densities (masses per unit of bulk volume in the reference configuration)

$$m_s = \rho^s (1 - \phi) J = \rho^s v_s \qquad m_f = \rho^f \phi J = \rho^f v_f \qquad (3.55)$$

for the solid and the fluid phase. Multiplication of equations (3.50) and (3.46) by $J$ and use of application of relation (3.41), together with the above definitions and properties, now yields the referential mass balances for the two phases:

$$\dot{m}_s = J r_s \qquad (3.56)$$

and

$$\dot{m}_f + \text{Div}(\rho^f Q) = -J r_s, \qquad (3.57)$$

while for the solute mass per unit referential volume $m_f C$, one has the referential balance equation

$$(m_f C)^\bullet + \text{Div}(C \rho^f Q) + \text{Div} J = -J r_s$$

or, by use of equation (3.57),

$$m_f \dot{C} + \rho^f Q \cdot \text{Grad} C + \text{Div} J = -(1-C) J r_s. \qquad (3.58)$$

The reader may have noted that hydrodynamic dispersion phenomena have been disregarded in setting up the last transport equation. This neglect is in fact fully justified in view of the small filter velocities that are to be expected for anticipated applications, where even at the highest recorded sedimentation rates molecular diffusion would by far outpace hydrodynamic dispersion, the relevant criterion for this to happen requiring Péclet numbers $Vd/D \le 1$ (see, e.g., Bear 1988, Chapter 10.4), where $V$ represents an average fluid particle velocity (say, 5 m/yr), $d$ an average grain diameter (say, 1 mm), and $D$ is the molecular diffusivity of the solute species in water (say, $10^{-8} \text{m}^2/\text{s}$) (see Section 3.4.2).
### 3.3.3 Force-Flux Relations

It is assumed that the filter velocity is governed by *Darcy's law*

$$ q = -\eta^{-1} K (\text{grad} p - \rho^f g). \qquad (3.59) $$

In the *referential description* this must take the form

$$ Q = -\eta^{-1} K^* (\text{Grad} p - \rho^f F^T g), \qquad (3.60) $$

where

$$ K^* = J F^{-1} K (F^{-1})^T \qquad (3.61) $$

is a *referential permeability tensor* that was first introduced by Biot (1972).

The diffusive mass flux vector of the solute species is assumed to obey a macroscopic form of *Fick's law* appropriate for the diffusion at dilute concentrations through a porous medium, which is

$$ j = -\rho^f D \text{grad} C, \qquad (3.62) $$

or, in the *referential description*,

$$ J = -\rho^f D^* \text{Grad} C, \qquad (3.63) $$

with

$$ D^* = J F^{-1} D (F^{-1})^T \qquad (3.64) $$

as an appropriate *referential effective diffusivity tensor*. The diffusivity tensor in *Fick's law*, as shown in equation (3.62), can in fact be decomposed into the three factors: porosity $\phi$, molecular diffusivity of the solute $D$ (in m²/s), and geometric conductivity $C$, according to

$$ D = \phi D C. \qquad (3.65) $$

The significance of the dimensionless symmetric tensor $C$ follows from the fact that the quantity $1-C$ can be interpreted as a *tortuosity tensor*.$^7$

<footnote>$^7$In the literature on porous media, the tensor $C$ itself is sometimes named “tortuosity” (Bear, 1988). As may easily be seen for a porous medium consisting of a bundle of straight rods or capillaries, the component of $C$ in the direction parallel to these will attain its maximum value of 1, implying zero tortuosity, whereas the component of $D$ in the same direction is simply determined by the molecular diffusivity and by the porosity, i.e., by the available cross section. But if the rod or capillaries were bent into a sinuous shape, for example, the actual diffusion path length would increase with decreasing wave length in relation to the straight distance in the overall direction, yielding the limit $l_{straight}/l_{sinuous} \to 0$. Thus, while the component of $C$ in the overall direction must vanish in this limit, the corresponding component of the tensor $1-C$ will approach its maximum value 1 so that it is in fact the latter quantity that deserves to be called tortuosity (see also Lehner 1979, for an analytical justification of this definition of tortuosity).</footnote>
### 3.3.4 Equations of Equilibrium

In this chapter we consider only slow, quasi-static motions of the two phases of a fluid-saturated porous medium in a given inertial frame. The equations of balance of linear and angular momentum for the “bulk” medium are therefore reduced to the equation of equilibrium

$$ \mathrm{div} \, \boldsymbol{\sigma} + \rho \boldsymbol{g} = \mathbf{0}, \qquad (3.66) $$

where $g$ denotes the acceleration of gravity, $\rho = (1 - \phi)\rho^s + \phi\rho^f$ the bulk density, and $\boldsymbol{\sigma}$ the symmetric Cauchy stress. This stress may be called *total stress* in the present context, where averaging considerations suggest the decomposition $\boldsymbol{\sigma} = \phi\boldsymbol{\sigma}^f + (1 - \phi)\boldsymbol{\sigma}^s$ in terms of the partial stresses of each phase. The partial stress $\boldsymbol{\sigma}^f$ is usually equated to the pore fluid pressure $-p\mathbf{1}$ in macroscopic theories that invoke a condition of local equilibrium. More relevant to our subject is the decomposition $\boldsymbol{\sigma} = \boldsymbol{\sigma}' - p\mathbf{1}$ into an *effective stress* $\boldsymbol{\sigma}'$ and pore fluid pressure, in terms of which equation (3.66) may be written

$$ \mathrm{div} \, \boldsymbol{\sigma}' + (\rho - \rho^f)\boldsymbol{g} - (\mathrm{grad} \, p - \rho^f \boldsymbol{g}) = \mathbf{0} \qquad (3.67) $$

or, upon recognizing that the last term in parentheses represents the driving force in Darcy's law (3.59),

$$ \mathrm{div} \, \boldsymbol{\sigma}' + (\rho - \rho^f)\boldsymbol{g} + \eta \mathbf{K}^{-1} \boldsymbol{q} = \mathbf{0}. \qquad (3.68) $$

Equation (3.66) may be cast in a referential description in terms of a nonsymmetric *nominal stress* $s$, related to the Cauchy stress $\boldsymbol{\sigma}$ by

$$ s = JF^{-1}\boldsymbol{\sigma}, \qquad \boldsymbol{\sigma} = J^{-1}Fs. \qquad (3.69) $$

This definition is consistent with the invariance requirement $\boldsymbol{\sigma}^T n da = s^T N dA$ for the force exerted upon a material surface element in the two configurations. By use of the divergence theorem one can therefore show that $J \mathrm{div} \, \boldsymbol{\sigma} = \mathrm{Div} \, s$, so that on multiplying equation (3.66) by $J$ one obtains the referential form of the equilibrium equation

$$ \mathrm{Div} \, s + m\boldsymbol{g} = \mathbf{0}, \qquad (3.70) $$

in which $m = m_f + m_s$ is a bulk mass density per unit referential volume.

Similarly, one can rewrite equation (3.68) in terms of a *nominal effective stress* $s' = JF^{-1}\boldsymbol{\sigma}'$ and the quantities entering Darcy's law (3.60):

$$ \mathrm{Div} \, s' + m'\boldsymbol{g} + \eta J(K^*F^T)^{-1}\boldsymbol{Q} = \mathbf{0}, \qquad (3.71) $$

where $m' = J(\rho - \rho^f) = v_s(\rho^s - \rho^f)$ is the so-called buoyant mass of the solid phase per unit referential volume.
### 3.3.5 Equations of State

The rate of change of the fluid mass content per unit referential volume can be written as

$$ \dot{m}_f = (\rho^f v_f) \cdot = \rho^f v_f \frac{\dot{\rho}^f}{\rho^f} + \rho^f \dot{v}_f. \qquad (3.72) $$

The density of the fluid phase satisfies an equation of state

$$ \rho^f = \rho^f(T, p, C) \qquad (3.73) $$

so that

$$ \frac{\dot{\rho}^f}{\rho^f} = -\alpha_f \dot{T} + \beta_f \dot{p} + \gamma_f \dot{C}, \qquad (3.74) $$

with

$$ \alpha_f(T, p, C) = -\frac{1}{\rho^f} \frac{\partial \rho^f}{\partial T} \Big|_{p,C}, \qquad (3.75) $$

$$ \beta_f(T, p, C) = \frac{1}{\rho^f} \frac{\partial \rho^f}{\partial p} \Big|_{T,C}, \qquad (3.76) $$

$$ \gamma_f(T, p, C) = \frac{1}{\rho^f} \frac{\partial \rho^f}{\partial C} \Big|_{p,T}. \qquad (3.77) $$

For solute concentrations close to equilibrium, so that the fluid density is well approximated by

$$ \rho^f = \rho^f[T, p, C_{eq}(T, p)], \qquad (3.78) $$

equation (3.74) may be simplified to

$$ \frac{\dot{\rho}^f}{\rho^f} = -\alpha_f^* \dot{T} + \beta_f^* \dot{p}, \qquad (3.79) $$

with

$$ \alpha_f^*(T, p) = \alpha_f - \gamma_f^* \partial C_{eq}/\partial T, \qquad (3.80) $$

$$ \beta_f^*(T, p) = \beta_f + \gamma_f^* \partial C_{eq}/\partial T, \qquad (3.81) $$

$$ \gamma_f^*(T, p) = \gamma_f[p, T, C_{eq}(T, p)]. \qquad (3.82) $$

The pore volume fraction $v_f$, which enters into equation (3.72), will undergo elastic as well as inelastic changes. These will be brought about by elastic and inelastic changes in bulk volume and the volume of the solid phase, the latter being due to IPS and solution or precipitation processes along the pore walls. We shall thus determine the pore volume fraction from the relation $v_f = J - v_s$ in terms of these quantities and, on adopting equation (3.79), write

$$ \dot{m}_f = \rho^f(J - v_s)(-\alpha_f^* \dot{T} + \beta_f^* \dot{p}) + \rho^f(\dot{J} - \dot{v}_s). \qquad (3.83) $$
The changes in the density of the solid phase with temperature and stress will contribute negligibly to the total deformation of a granular aggregate that is compacting by IPS. We therefore disregard any changes in $\rho^s$ and approximate the rate of change of the solid mass $m_s = \rho^s v_s$ in the balance equation (3.42) by $\rho^s \dot{v}_s$, so that

$$ \dot{v}_s = Jr_s / \rho^s = \frac{1}{V_0} \int_{S_{sf}} (w_n - v_n) da - \frac{1}{V_0} \int_{S_{gb}} [[w_n - v_n]] da \quad (3.84) $$

is obtained from equation (3.44) when the density of the solid phase is uniformly equal to $\rho^s$. This result is seen to be consistent with equation (3.40) when the last term (which equals $-v_s \dot{\rho}^s$) is omitted from that equation and upon noticing that $w_n$ is continuous across $S_{gb}$. It allows $\dot{v}_s$ to be determined through the use of the kinetic relations (3.26) and (3.25).

In this section we have discussed certain key elements of a continuum theory of pressure solution creep. The treatment of the constitutive relations remains clearly incomplete (i.e., focused exclusively on pressure solution creep) and the same is true for certain important aspects of a general field theory. Although some of these, such as the coupled transport of heat and mass, would represent standard developments, others would point in the direction of substantial future efforts, e.g., the major task of modeling the evolution of a general grain fabric numerically. Here we have deliberately kept the theoretical discussion within the scope of the applications that are within reach at this point. In addressing these, Section 3.4 will begin with a description of the simplest possible microstructure—that of a simple cubic packing of identical spheres. Although strongly idealized, this structure offers the advantage of allowing a complete formulation of the theory in closed form, from a grain-scale description to a macroscopic continuum theory.

## 3.4 Compaction Creep for a Simple Microstructure

### 3.4.1 Unit Cell Behavior

As a first application of the theory developed in the foregoing section, let us now consider the compaction by IPS for a simple possible microstructure (i.e., that of a simple cubic packing of spherical grains with uniform initial grain diameter $d$). We may associate a cube of volume $d^3$ as a unit cell with this initial grain diameter, such that the centers of its faces mark the contact points of neighboring spheres. This will be the largest possible initial unit cell before development of grain truncations by IPS. We shall assume here that the directions of the surface normals to the faces of this cell coincide with the principal directions of the macroscopic stress acting
on the granular aggregate throughout its deformation history. We therefore refer to these directions simply as principal directions.

By the process of compaction by IPS, as described for the single grain-to-grain contact of Figure 3.3, the edges of the unit cell will have shrunk at a certain moment to the lengths $l_\nu d$, $\nu = 1, 2, 3$, the now prismatic cell truncating the grain along six pairwise orthogonal grain-to-grain contacts, the cell volume having been reduced to $V = l_1 l_2 l_3 d^3$ as shown in Figure 3.5.

To introduce a certain additional degree of freedom into the description of the initial grain shape, we shall view the actual initial state of the chosen grain assembly as one in which the grain contact areas have already acquired a certain finite size, precisely as if the grains had experienced some fictitious episode of initial IPS, producing truncations along the principal planes and a corresponding reduction in length of the edges of the unit cell from the maximum length $d$ to an actual initial length $l_\nu^0 d$, $\nu = 1, 2, 3$. We further assume the grain diameter to be equal to $d$ in this initial state, which we shall henceforth identify with the reference state of our grain assembly. The volume of the unit cell in the reference state is thus given by $V_0 = l_1^0 l_2^0 l_3^0 d^3$.

Our aim is now to obtain expressions for the creep rate, i.e., the deformation rate $D_{ij}$ of the unit cell and for rate of change of the solid volume fraction $\dot{v}_s$ in terms of the relevant macroscopic driving force. The latter will be found to differ for closed- and open-system conditions. Clearly, the unit cell response will be equivalent to the homogeneous deformation response of a macroscopic continuum that possesses the simple microstructure of the assumed cubic packing of spheres.

<figure><img src="image_6.png" /><figcaption>**Figure 3.5** ► Single grain from a simple cubic packing of identical grains showing six pairwise orthogonal grain truncations.</figcaption></figure>
The volume of the unit cell at the current time $t$ equals $V = l_1l_2l_3d^3 = \lambda_1\lambda_2\lambda_3V_0$, where the true stretches $\lambda_\nu$ from the reference configuration are related to the fictitious initial and current stretches $l_\nu^0$ and $l_\nu$ by

$$l_\nu = \lambda_\nu l_\nu^0. \qquad (3.85)$$

From this we have $J = V/V_0 = \lambda_1\lambda_2\lambda_3$ and consequently

$$\dot{J} = J \sum_{\nu=1}^{3} \frac{\dot{\lambda}_{\nu}}{\lambda_{\nu}} = J \sum_{\nu=1}^{3} D_{\nu}, \qquad (3.86)$$

the $D_\nu$ denoting the eigenvalues of the deformation rate $D_{ij}$ for the known principal directions. Focusing once more on Figure 3.3, we note that (since the grain is assumed to be rigid) we must have

$$\dot{\lambda}_{\nu} = \llbracket v_{\nu} \rrbracket / (l_{\nu}^{0} d) \qquad (3.87)$$

for the rate of shortening in each of the three principal directions, $\nu = 1, 2, 3$, normal to a grain contact; this will suffice to determine the deformation rate in the present case.

For reasons of simplicity we make the fundamental assumption that free-face solution or precipitation will preserve the spherical shape of the free grain surface, leading to a simple shrinkage or growth of the grain radius $R(t)$ at the rate $\dot{R}(t)$, with $R_0 = R(t = 0) = d/2$ as initial radius. As long as the grain truncations remain circular (we shall assume that the deformation stays within this range), the solid volume fraction is given by the expression

$$v_s = V_s/V_0 = \frac{\pi}{6l_1^0 l_2^0 l_3^0} \left\{ r^3 - \frac{1}{2} \sum_{\nu=1}^{3} (r - l_\nu)^2 (2r + l_\nu) \right\}, \qquad (3.88)$$

where $r = R/R_0$ stands for the normalized grain radius. Note that $l_\nu$ must ultimately be replaced in equation (3.88), as well as in all subsequent appearances, by its expression (3.85) in terms of the true stretch $\lambda_\nu$.

Because the simplicity of the assumed grain assembly gave us the explicit expression (3.88) for $v_s$, we can circumvent the determination of the rate $\dot{v}_s$ from equation (3.84) (a step that cannot be avoided in general for an incrementally evolving, irregular grain fabric) and obtain $\dot{v}_s$ directly by differentiation of equation (3.88). The result is

$$\dot{v}_s = \frac{\pi}{4l_1^0 l_2^0 l_3^0} \left\{ \sum_{\nu=1}^{3} (r^2 - l_\nu^2) l_\nu^0 \dot{\lambda}_\nu + \left[ 2r^2 - \sum_{\nu=1}^{3} \{(r - l_\nu)^2 + r^2 - l_\nu^2\} \right] \dot{r} \right\}. \qquad (3.89)$$
The rate of change of the pore volume fraction may be determined immediately in terms of these results from the relation $\dot{v}_f = J - \dot{v}_s$.

To obtain $\dot{\lambda}_v$ we only have to substitute equation (3.18)—with $a_v$ instead of $a$ for the contact radius—in equation (3.87) to arrive at the following general kinetic relation for the rate of grain convergence

$$l_v^0 \dot{\lambda}_v = -(\rho^s \Omega^s / kT)(\bar{\mu}_v^s - \mu)/[\tau_s + (4a_v^2/d^2)\tau_D]. \quad (3.90)$$

Here and subsequently we denote by $\mu \equiv \mu(a)$ the locally (on the macroscale) uniform chemical potential of the solute in the pore space.

For the description of open systems, we shall employ the approximate result of equation (3.25) in equation (3.87), giving

$$l_v^0 \dot{\lambda}_v = [(\Omega^s / kT)(\bar{\sigma}_v + p) + (C - C_{eq})/C_{eq}]/[\tau_s + (4a_v^2/d^2)\tau_D]. \quad (3.91)$$

To obtain the rate of change $\dot{r}$ of the grain radius, we fix our frame of reference on the grain shown in Figure 3.5 to note that $\dot{r} = 2w_n/d = 2r_{sf}/(\rho^s d)$ if we make the assumption of uniform growth or shrinkage of the grain radius. To remain consistent with this assumption, we also replace the local normal component of the chemical potential along $S_{sf}$ by its average $\tilde{\mu}_n^s$, taken over $S_{sf}$, so that equation (3.11) yields

$$\dot{r} = -\tau_1^{-1}(\rho^s \Omega^s / kT)(\tilde{\mu}_n^s - \mu), \quad (3.92)$$

whereas relation (3.26) produces the approximate result

$$\dot{r} = \tau_1^{-1}(C - C_{eq})/C_{eq}. \quad (3.93)$$

Here we have introduced a third characteristic time that is associated with the interface reaction along the free faces of a grain and is defined by

$$\tau_1^{\pm} = d/(2K_{\pm}^{sf}). \quad (3.94)$$

We must now seek an expression for the local (grain boundary) average stress $\bar{\sigma}_v$ in terms of the macroscopic average stress, which we shall denote by $\langle \sigma_{ij} \rangle$ for clarity, departing from the above usage. Similarly, let $\langle \sigma_{ij} \rangle_s = V^{-1} \int_{V_s} \sigma_{ij} dV$ denote the macroscopic phase-averaged stress in the solid phase. Fixing the origin of a Cartesian coordinate system at the center of the unit cell, the phase-averaged stress may also be expressed by the surface integral $V^{-1} \int_{\partial V_s} \sigma_{ik} n_k x_j da$ taken over the entire (closed) boundary $\partial V_s = S_{sf} + S_{gb}$ of the grain, where $n_k$ and $x_j$ denote components of the outward unit normal to $\partial V_s$ and of the position vector, $S_{gb}$ comprises the six intergranular contacts for a unit cell, and $S_{sf}$ stands for
the pore walls. Using the fact that the only traction along $S_{sf}$ is supplied by the uniform pore fluid pressure $p$, we split the integral over $\partial V_s$ and write

$$ \begin{aligned} \langle \sigma_{ij} \rangle_s &= -p \frac{1}{V} \int_{S_{sf}} n_i x_j da + \frac{1}{V} \int_{S_{gb}} \sigma_{ik} n_k x_j da \\ &= -p \frac{1}{V} \int_{\partial V_s} n_i x_j da + \frac{1}{V} \int_{S_{gb}} (\sigma_{ik} + p\delta_{ik}) n_k x_j da, \end{aligned} \quad (3.95) $$

and since the first integral on the right equals $(1 - \phi)p\delta_{ij}$, where $\phi = V_f/V$ denotes the aggregate porosity, the relationship

$$ \langle \sigma_{ij} \rangle + p\delta_{ij} = \langle \sigma_{ij} \rangle_s + (1 - \phi) p\delta_{ij} = \frac{1}{V} \int_{S_{gb}} (\sigma_{ik} + p\delta_{ik}) n_k x_j da \quad (3.96) $$

exists between the integrated effective grain boundary traction and the macroscopic effective stress. Note that the total stress is defined as the sum of the partial or phase-averaged stresses of the solid and fluid phases according to $\langle \sigma_{ij} \rangle = \langle \sigma_{ij} \rangle_s + \langle \sigma_{ij} \rangle_f = \langle \sigma_{ij} \rangle_s - \phi p \delta_{ij}$. In fact, because we assume that the macroscopic principal stress directions remain perpendicular to the three orthogonal pairs of grain-to-grain contacts, we need to consider only a relation in the principal stresses of the form (in which no summation is implied over a repeated index $v$)

$$ \langle \sigma_v \rangle + p = \frac{1}{V} \int_{S_{gb}} (\sigma_v + p) x_v n_v da, \quad v = 1, 2, 3. \quad (3.97) $$

From the geometry of the unit cell it follows that $x_v n_v = l_v d/2$ as well as $V = l_1 l_2 l_3 d^3$ and $a_v^2 = (r^2 - l_v^2)d^2/4$ for the current volume of the unit cell and the radius of the circular grain contacts with surface normal in the $x_v$-direction. The integral in equation (3.97) therefore yields the following expression for the average effective normal stress on a contact in terms of the macroscopic principal effective stress:

$$ \bar{\sigma}_v + p = \frac{4}{\pi} [l_1 l_2 l_3 / (r^2 - l_v^2) l_v] (\langle \sigma_v \rangle + p), \quad v = 1, 2, 3. \quad (3.98) $$

This is the desired relationship that, upon substitution in equation (3.91), furnishes the result (with $l_v = l_v^0 \lambda_v$, $v = 1, 2, 3$)

$$ l_v^0 \dot{\lambda}_v = \frac{1}{\tau_s + (r^2 - l_v^2)\tau_D} \left[ \frac{4l_1 l_2 l_3}{\pi(r^2 - l_v^2)l_v} \left( \frac{\langle \sigma_v \rangle + p}{kT/\Omega^s} \right) + \frac{C - C_{eq}}{C_{eq}} \right]. \quad (3.99) $$

In the special case of a simple cubic packing undergoing isotropic compaction under an effective all-round compressive stress $\langle \sigma \rangle + p$, one has $l \equiv l_1 = l_2 = l_3$ and $l_0 \equiv l_1^0 = l_2^0 = l_3^0$, so that $a_v^2 = a^2 = (r^2 - l^2)d^2/4$ and relation (3.99) specializes to (with $l = l_0 \lambda$)

$$ l_0 \dot{\lambda} = \frac{1}{\tau_s + (r^2 - l^2)\tau_D} \left[ \frac{4l^2}{\pi(r^2 - l^2)} \left( \frac{\langle \sigma \rangle + p}{kT/\Omega^s} \right) + \frac{C - C_{eq}}{C_{eq}} \right]. \quad (3.100) $$
The description of open-system compaction creep by IPS in a simple cubic packing of spheres is completed with the determination of $\dot{v}_s$ by substitution of the expressions (3.93) and (3.99) for $\lambda_v$ and $\dot{r}$ in equation (3.89).

The pressure solution creep law (3.99) was first obtained by Lehner (1990, 1995, 1997), who also drew attention to the following particular features. First of all, the fact that equation (3.99) applies to open systems, that is, systems with significant solute import or export into or from an REV (here the unit cell), so that the concentration $C$ of the solute is coupled to its transport. This makes for the appearance in the thermodynamic driving force (the sum of terms within brackets) of a chemical as well as a mechanical term. We note here that the absence of the chemical term from most published pressure solution creep relations suggests that these should apply either under closed-system conditions, as will be seen, or else under conditions of sufficiently fast reprecipitation of the dissolved material such that significant supersaturations could not exist. We note that in the present circumstances the chemical term ensures a physically consistent behavior at vanishing effective pressure, when $\lambda_v$ need not vanish (as is usually assumed) and compaction creep can continue as long as the solution phase remains undersaturated ($C < C_{eq}$). Relation (3.99) also predicts that IPS should cease to operate at a certain supersaturation. Here we must caution, however, that even if that supersaturation level were attained, the prediction itself may well turn out too inaccurate, principally because the complex process of grain boundary healing is not addressed by the present model.

Secondly, although linear in the stress and concentration, relation (3.99) displays a strong geometrical nonlinearity. In other words, at fixed values of these variables the creep rate will depend strongly on the accumulated strain through a strong dependence on the current size of the intergranular contacts. This geometrical nonlinearity has been taken into account for some time in experimental studies of compaction creep behavior caused by IPS. For example, care has been taken in such studies to determine the dependence of creep rates on stress at a fixed strain (see, e.g., the work of Spiers and Schutjens (1990), which also contains a theoretical analysis of this nonlinearity).

A further observation relates to the appearance of the characteristic times $\tau_s$ and $\tau_D$ in the leading factor of equation (3.99). Clearly, the larger one of the two will tend to govern the rate of the overall process, as it should in the case of serial processes. As may be expected, however, $\tau_D$ is weighted by the grain contact size, becoming less effective at smaller relative contact sizes. A further important feature of these characteristic times, apparent from their definitions, is the different dependence on grain size that they exhibit. It follows that grain boundary diffusion-controlled creep rates should display a dependence on the third power of the grain size, whereas dissolution-controlled creep rates should vary linearly with $d$. This theoretically expected result has been exploited for some time by experimentalists so as to determine the rate-controlling process in experimentally compacted grain aggregates of various salts (e.g., NaCl, KCl, $NaNO_3$) that served as rock analogues (Spiers and Schutjens, 1990; Spiers et al., 1990; Spiers and Brzesowsky, 1993).
Let us now consider briefly the compaction of the unit cell of Figure 3.5 under closed-system conditions. These we shall define by the simple requirement of conservation of solid mass or—because the density of the solid is assumed to remain constant—the requirement $\dot{v}_s = 0$, which from equation (3.89) is seen to demand the relationship

$$ \dot{r} = \frac{-\sum_{\nu=1}^{3} (r^2 - l_{\nu}^2) l_{\nu}^0 \dot{\lambda}_{\nu}}{2r^2 - \sum_{\nu=1}^{3} [(r^2 - l_{\nu}^2) + (r - l_{\nu})^2]} \quad (3.101) $$

between $\dot{r}$ and the $\dot{\lambda}_{\nu}$. A second relationship between these rates, obtained from equations (3.91) and (3.93) by elimination of the concentration term $(C - C_{eq})/C_{eq}$ and by use of equation (3.98) is

$$ l_{\nu}^{0} \dot{\lambda}_{\nu} = \left\{ \frac{\frac{4}{\pi} [l_1 l_2 l_3 / (r^2 - l_{\nu}^2) l_{\nu}] (\Omega^s / kT) (\langle \sigma_{\nu} \rangle + p) + \tau_1 \dot{r}}{\tau_s + (r^2 - l_{\nu}^2) \tau_D} \right\} / [\tau_s + (r^2 - l_{\nu}^2) \tau_D]. \quad (3.102) $$

Next, we substitute this result in relation (3.101), to obtain

$$ \dot{r} = \frac{-(4/\pi) \sum_{\nu=1}^{3} (l_1 l_2 l_3 / l_{\nu}) (\Omega^s / kT) (\langle \sigma_{\nu} \rangle + p) / [\tau_s + (r^2 - l_{\nu}^2) \tau_D]}{2r^2 - \sum_{\nu=1}^{3} [r^2 - l_{\nu}^2 + (r - l_{\nu})^2] + \tau_1 \sum_{\nu=1}^{3} (r^2 - l_{\nu}^2) / [\tau_s + (r^2 - l_{\nu}^2) \tau_D]} \quad (3.103) $$

for the rate of change of the (nondimensional) grain radius under closed-system conditions. This expression for $\dot{r}$ may now be used in equation (3.102) to obtain $\dot{\lambda}_{\nu}$.

The previous system of equations for $\dot{r}$ and $\dot{\lambda}_{\nu}$ takes on a simpler form in the case of isotropic compaction, with $l_{\nu} = l$, $\lambda_{\nu} = \lambda$, and $\langle \sigma_{\nu} \rangle = -P$ ($\nu = 1, 2, 3$), when

$$ \dot{r} = \frac{(4/\pi)l^2(\Omega^s/kT)P_{\text{eff}}}{(2rl - 4r^2/3)[\tau_s + (r^2 - l^2)\tau_D] + (r^2 - l^2)\tau_1} \quad (3.104) $$

and, with $l = l_0 \lambda$,

$$ l_0 \dot{\lambda} = \frac{-(4/\pi)l^2(2rl - 4r^2/3)(\Omega^s/kT)P_{\text{eff}}}{(r^2 - l^2)\{(2rl - 4r^2/3)[\tau_s + (r^2 - l^2)\tau_D] + (r^2 - l^2)\tau_1\}} \quad (3.105) $$

in terms of a macroscopic effective pressure $P_{\text{eff}} = P - p$. Note that the assumptions of nonexistent intergranular precipitation and of circular grain truncations impose the constraints $\lambda \le 1$ and $1 \le r/\lambda l_0 \le \sqrt{2}$, respectively, on admissible values of $\lambda$ and $r$ in the last two expressions.

### 3.4.2 One-Dimensional Compaction of a Thick Layer

The previous description of the unit cell compaction behavior will now be put to use for the solution of a first boundary value problem, that of the one-dimensional compaction of a thick sedimentary layer. Our analysis of this problem can be seen as an attempt to explore the theoretically predicted process of compaction of a thick
sedimentary sequence during a relatively short interval of time. We prefer to think of it as a first test case, allowing us to explore the consequences of certain model assumptions and to gain a perspective on necessary future developments. While our problem may be placed in the context of earlier work on rock densification by pressure solution (Angevine and Turcotte, 1983; Lemee and Gueguen, 1996; Fowler and Yang, 1999), we shall not attempt a direct comparison in this article. Rather, we shall focus on a simple problem of the deformation of a thick layer of sediment (sandstone) in which the process of IPS is suddenly switched on (i.e., allowed to occur) only after its deposition. We shall also make a number of simplifying assumptions, specific to the problem at hand, that can later be abandoned in studies of a broader kind that are based on the general theoretical framework of this chapter.

To proceed to the formulation of the governing equations of our problem, let us now introduce the spatial and referential Cartesian coordinate systems {$x_1, x_2, x_3$} and {$X_1, X_2, X_3$}, where $x_3$ and $X_3$ are pointing vertically upward. The particle motion that occurs during one-dimensional compaction in the negative $x_3$-direction is described by the component equations $x_1 = X_1$, $x_2 = X_2$, and $x_3 = x_3(X_3, t)$, from which the matrix of the deformation gradient $F_{i\alpha} = \partial x_i / \partial X_\alpha$ and that of its inverse are readily determined as

$$ (F_{i\alpha}) = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \lambda \end{pmatrix}, \quad (F_{i\alpha}^{-1}) = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \lambda^{-1} \end{pmatrix}, \quad (3.106) $$

where $\lambda = \partial x_3 / \partial X_3$ is the stretch (here in fact shortening) in the vertical direction. We also note that $F = F^T$ and $J = \det F = \lambda$.

Consider next the equation of equilibrium (3.70). Because we assume all derivatives with respect to $X_1$ and $X_2$ to vanish, we must have

$$ \mathrm{Div} s = \frac{\partial}{\partial X_3} \left( J F^{-1} \boldsymbol{\sigma} \right)_{33} e_3 = \frac{\partial}{\partial X_3} \left( \lambda F^{-1}_{3j} \sigma_{j3} \right) e_3 = \frac{\partial \sigma_{33}}{\partial X_3} e_3, \quad (3.107) $$

and after scalar multiplication of equation (3.70) by the unit vector $e_3$ in the positive (upward oriented) $X_3$-direction, there follows the only nontrivial equilibrium equation for the vertical direction

$$ \frac{\partial \sigma_{33}}{\partial X_3} - mg = 0, \quad (3.108) $$

in terms of the (Cauchy) normal stress $\sigma_{33}$ which, by the above assumptions, must also be a principal stress.

The equivalent form of equation (3.71) in the effective stress is

$$ \frac{\partial \sigma'_{33}}{\partial X_3} - (\rho^s - \rho^f) v_s g + \frac{\eta}{K_3^*} Q_3 = 0. \quad (3.109) $$
Here we assume that the principal axes of the symmetric permeability **K** and effective diffusivity **D** coincide with the coordinate directions, so that by virtue of equations (3.61) and (3.64) the matrices of the corresponding referential quantities are of the form

$$ (K_{ij}^*) = \lambda \begin{pmatrix} K_1 & 0 & 0 \\ 0 & K_2 & 0 \\ 0 & 0 & \lambda^{-2} K_3 \end{pmatrix}, \quad (D_{ij}^*) = \lambda \begin{pmatrix} D_1 & 0 & 0 \\ 0 & D_2 & 0 \\ 0 & 0 & \lambda^{-2} D_3 \end{pmatrix}. \qquad (3.110) $$

The only nonvanishing components of the referential filter velocity and diffusive mass flux therefore satisfy the force-flux relations

$$ Q_3 = -\eta^{-1} K_3^* \left( \frac{\partial p}{\partial X_3} + \lambda \rho^f g \right) \qquad (3.111) $$

and

$$ J_3 = -\rho^f D_3^* \frac{\partial C}{\partial X_3}. \qquad (3.112) $$

These enter into the one-dimensional forms of the balance equations (3.57) and (3.58).

## Problem Statement for Hydrostatic Pore Pressures

Let us now focus on an interesting situation that is likely to occur in sandstones of relatively high permeability when compaction rates by IPS are too slow to produce significant nonhydrostatic excess fluid pressures. The theoretical description of one-dimensional compaction may then be further simplified by forsaking the use of Darcy's law (3.111) and treating the flux $Q_3$ as an unknown dependent variable of the problem, while assuming an initial hydrostatic distribution of the pore fluid pressure to remain unaltered during compaction. In consequence, $Q_3$ is carried along in the mass balance equations, but deleted from the equilibrium equation (3.109), which is thereby reduced to its static form

$$ \frac{\partial \sigma'_{33}}{\partial X_3} - \Delta \rho v_s g = 0, \qquad (3.113) $$

where $\Delta \rho = \rho^s - \rho^f$.

We further simplify the above general formulation by assuming the temperature field to be given in terms of a fixed geothermal gradient, by taking the fluid and solid mass densities to be constant, and by stipulating dilute solute concentrations ($C \ll 1$). Using an appropriately simplified equation of state (3.83), together with
relation (3.84), we now introduce equation (3.112) into equation (3.58) to obtain
the following system of equations in $Q_3$ and $C$ in one dimension:

$$
-\frac{\partial Q_3}{\partial X_3} = \left(\frac{\rho^s}{\rho^f} - 1\right) \dot{v}_s + \dot{\lambda} \quad (3.114)
$$

and

$$
-Q_3 \frac{\partial C}{\partial X_3} + \frac{\partial}{\partial X_3} \left( D_3^* \frac{\partial C}{\partial X_3} \right) = \frac{\rho^s}{\rho^f} \dot{v}_s + (\lambda - v_s) \dot{C}. \quad (3.115)
$$

Here $\lambda$, $\dot{\lambda}$, $v_s$, and $\dot{v}_s$ remain to be determined.

The effective diffusivity of the porous material is allowed to depend solely
on temperature (through the molecular diffusivity $\mathcal{D}$), i.e., we shall disregard any
changes in the relevant component $D_3^*$ that result from the changes in the pore
geometry during compaction.<fnref n="8" /> However, we shall disassociate the numerical value
used for the geometrical part $\phi C_3 / \lambda$ of $D_3^*$ from the particular microstructure of
a simple cubic packing and assume an empirical Archie-type relation of the form
$D_3^* = \mathcal{D} \phi^2$. In fact, we shall simply take

$$
D_3^* = \mathcal{D} \phi_0^2, \tag{3.116}
$$

the porosity being set equal here to an assumed initial, depth-dependent porosity
$\phi_0$ specified presently by equation (3.132).

To complete the formulation of our problem, we shall now adopt the above
model of a simple cubic packing of spheres for a grain-scale description of the
sedimentary rock. Although this implies a drastic simplification of reality, the
possibility to test a model of IPS that rests on a clear physical foundation matters
most to us at this stage. Imposing the constraint of one-dimensional compaction,
we assume $l_1^0 = l_2^0 = l_3^0 = l_0$ so that $l_1 = l_2 = l_0$, while $l_3 = l = l_0\lambda$; putting
these into equations (3.88), (3.89), and (3.99), we now have

$$
v_s = \frac{\pi}{6l_0^3} [r^3 - (r-l_0)^2(2r+l_0) - \frac{1}{2}(r-l_0\lambda)^2(2r+l_0\lambda)] \quad (3.117)
$$

and

$$
\dot{v}_s = \frac{\pi}{4l_0^3} \{(r^2 - l_0^2\lambda^2)l_0\dot{\lambda} + 2[l_0\lambda - 2(r-l_0)]r\dot{r}\}, \quad (3.118)
$$

<footnote n="8">⁸In the present case equation (3.110) implies $D_3^* = \mathcal{D}\phi C_3/\lambda = \mathcal{D}(1 - v_s/\lambda)C_3/\lambda$ according to definition (3.65), which shows that the use of the referential diffusivity will in fact compensate to some extent for a reduction of $C_3$ through the simultaneous decrease in $\lambda$ during deformation.</footnote>
where the rates $\dot{r}$ and $\dot{\lambda}$ must be found from the fundamental expressions (3.93) and (3.99), the latter assuming the form

$$l_0 \dot{\lambda} = \frac{1}{\tau_s + (r^2 - l_0^2 \lambda^2) \tau_d} \left[ \frac{4l_0^2}{\pi(r^2 - l_0^2 \lambda^2)} \left( \frac{\sigma_{33} + p}{kT/\Omega^s} \right) + \frac{C - C_{eq}}{C_{eq}} \right]. \quad (3.119)$$

It is through this last relation for $\dot{\lambda}$ that the coupling to the (effective) stress occurs. Note that the notation has now been simplified again, letting $\sigma_{33}$ stand for the component of the macroscopic total stress that appears in the equilibrium equation (3.108), in fact the principal stress in the vertical direction.

To return to the formulation of our boundary value problem, let us now consider a layer of sandstone with initial thickness $h(0) = H$, as measured vertically upward from a fixed (impermeable) base at $x_3 = X_3 = 0$. At time $t$ the top of the layer will have subsided to a level $x_3 = h(t)$, the sea level remaining at $x_3 = H$ at all times (Figure 3.6). No sediment is added at the top of the layer during the entire time interval studied.

The initial state of the layer is one in which it may be assumed to have undergone a substantial amount of compaction by mechanical and possibly also chemical processes of early diagenesis. We shall allow for this in a simple manner by assuming a certain initial distribution of the porosity, corresponding to an empirical porosity-depth profile, which we shall translate into a corresponding fictitious initial strain $l_0$, as discussed previously. We do not wish to attach an inordinate significance to certain quantitative aspects of our model, such as the absolute porosity values.

<figure><img src="image_7.png" /><figcaption>**Figure 3.6** ► A layer of sandstone undergoing vertical gravitational compaction by IPS; its initial thickness $h(0) = H$ is reduced to $h(t) < H$ at times $t > 0$ by the motion $x_3 = x_3(X_3, t)$ of horizons marked by their initial position $X_3$. Coordinates $X_3$ and $x_3$ are used in the referential description (a) and the spatial description (b).</figcaption></figure>
Rather, we shall focus on the changes that are predicted to occur once pressure solution has been “switched on.”

Figure 3.6 illustrates the difference between a *spatial description* (b) in terms of the vertical coordinate $x_3$ and a *referential description* (a) in terms of the coordinate $X_3$, corresponding to the initial position of a compacting horizon currently at $x_3$. When employing the second, any given function of $x_3$ must be translated into a function of $X_3$ by use of the motion $x_3 = x_3(X_3, t)$. In particular, we shall have to apply this mapping to the assumed temperature distribution

$$T(x_3, t) = T_0 + T'[h(t) - x_3], \qquad (3.120)$$

where $T_0 = 277\text{K}$ is the temperature at the sea floor and $T' = 30 \text{ K/1000 m}$ is a fixed geothermal gradient. The known temperature distribution allows the determination of a number of temperature-sensitive quantities that are specified here by the following functions (T in Kelvin):

* The solubility $C_{eq}$ (expressed as a mass fraction) of quartz in water at neutral pH:

$$\log_{10} C_{eq} = \log_{10} M^s - 0.254 - \frac{1107.12}{T}, \qquad (3.121)$$

from Rimstidt 1997 solubility data (given in mol/kg) upon multiplication by the molar mass of quartz, $M^s = 0.055 \text{ kg/mol}$.

* The rate “constants” $K_-^{sf}$ and $K_+^{sf}$ (in m/s) for the precipitation and solution reactions along the quartz-water interface (taken from Rimstidt and Barnes 1980):

$$\log_{10} K_-^{sf} = \log_{10} \bar{V}^s - 0.707 - \frac{2598}{T}, \qquad (3.122)$$

$$\log_{10} K_+^{sf} = \log_{10} \bar{V}^s + 1.174 - 0.002028 \times T - \frac{4158}{T}, \qquad (3.123)$$

obtained by multiplying the rate constants $k_-$ and $k_+$, as given by Rimstidt and Barnes (1980) in units of mol/m²/s, by the molar volume of quartz, $\bar{V}^s = 2.27 \times 10^{-5} \text{ m}^3/\text{mol}$.

No information is available on the “rate constant” $K^{gb}$ that was introduced in relation (3.15) by analogy with (3.11). In the following, we shall therefore equate $K^{gb}$ to $K_+^{sf}$.

* The molecular diffusivity $D$ (in m²/s) of the solute in water:

$$D = D_0 \exp\{-15000/RT\}. \qquad (3.124)$$

Here we follow Nakashima (1995) in assuming an activation energy of 15 KJ/mol for this process, while choosing the value $9 \times 10^{-7} \text{ m}^2/\text{s}$ for $D$, to match Applin (1987) value of $2.2 \times 10^{-9} \text{ m}^2/\text{s}$ for the diffusivity at ambient conditions.
► The grain boundary diffusion factor $\delta D_{gb}$ (in m$^3$/s) is estimated by adopting a range of values for the effective grain boundary width $\delta$ (which includes a grain boundary porosity factor), while taking $D_{gb} = 0.1D$, i.e., essentially assuming equal molecular diffusivities for the solute in the grain boundaries and the free pore space, but reducing the former to $D_{gb}$ by a tortuosity factor of approximately 1/10. This assumption is consistent with our view that the grain boundaries in sandstones of the type shown in Figure 3.1 possess a coarse island-channel structure that is infiltrated by the pore fluid. The grain boundary width $\delta$ is therefore taken to be approximately equal to the height of the channels and should not be confused with the thickness of a fluid film that may actually exist temporarily between the islands or asperities of two grains, where these define an actual grain-to-grain contact on the microscopic scale. Although the rate of undercutting of such asperity contacts is likely to be controlled by a local dissolution rate constant, the overall diffusive transport along the grain boundary of the solute species produced by this undercutting process must be governed by its molecular diffusivity in the network of channels as well as by the channel geometry. Mindful of the fact that the physical structure of the grain boundaries remains a controversial topic among investigators of rather diverse empirical backgrounds of individual investigators, we do not wish to assume a dogmatic position on this point. We therefore propose to examine the consequences for our model predictions of a variation in the effective grain boundary width $\delta$ between 2 nm, as appropriate for a thin fluid-film model, and 200 nm, as would apply to a coarse island-channel structure.

► From definitions (3.19), and (3.94) it is apparent that the above temperature functions will affect the characteristic times $\tau_D$, $\tau_S$, and $\tau_I$ through the appropriate rate constants and—in the case of $\tau_D$—through the solubility $C_{eq}$ as well.

The governing differential equations for our problem are equations (3.113), (3.114), and (3.115), together with equations (3.93) and (3.117) through (3.119), and with given functions (3.121) through (3.124) of the temperature distribution (3.120). To these are adjoined the following initial and boundary conditions:

$$C(X_3, 0) = C_{eq}(X_3, 0), \qquad (3.125)$$

$$r(X_3, 0) = 1, \qquad (3.126)$$

$$\lambda(X_3, 0) = 1, \quad \text{for } 0 < X_3 < H, \qquad (3.127)$$

and

$$C(X_3 = H, t) = C_{eq}(H, 0), \qquad (3.128)$$

$$\frac{\partial C}{\partial X_3}(X_3 = 0, t) = 0, \qquad (3.129)$$

$$Q_3(X_3 = 0, t) = 0. \qquad (3.130)$$
The initial solute concentration (3.125) is inconsistent with condition (3.129) at the basis of the domain if condition (3.121) is adopted together with condition (3.120). To satisfy the boundary condition (3.129), the distribution (3.121) is therefore amended in the following ad hoc manner in the subsequent numerical calculations: Below a fraction of the current column height (typically 10%), the equilibrium concentration (3.125) is changed to a parabolic distribution with a zero slope at the basis. The two parameters of the parabola are selected such as to ensure the continuity of the concentration and of its first spatial derivative.

We note further that by virtue of equation (3.117) the conditions (3.126) and (3.127) imply the relationship

$$ \phi_0(X_3) = 1 - v_s(X_3, 0) = 1 - (\pi/12)(9l_0^{-2} - 4l_0^{-3} - 3) \quad (3.131) $$

between the porosity in the initial (reference) state and $l_0$. This enables us to mimic an exponential initial porosity-depth distribution of the exponential Athy-type

$$ \phi_0(X_3) = \phi^* \exp\{-\kappa(H - X_3)\}, \quad (3.132) $$

via a suitable depth dependence of $l_0(X_3)$, which is obtained as the solution of the cubic equation that results from a substitution of equation (3.132) in equation (3.131). Here $\kappa$ is a decay constant and $\phi^*$ denotes the constant surface porosity at the top of the layer. Following Sclater and Christie (1980), we take $\kappa = 2.7 \times 10^{-4} \text{m}^{-1}$ for a sandstone, and $\phi^* = 1 - \pi/6$, corresponding to $l_0 = 1$ at the sedimentation boundary $X_3 = H$. For these parameter values it turns out that relation (3.131) is accurately approximated by the linear relation $\phi_0 = al_0 + b$, with $a = 1.632$ and $b = -1.1556$ within the relevant porosity and depth range. The following initial distribution of $l_0$ may therefore be used in the example

$$ \begin{aligned} l_0 &= (\phi^*/a) \exp\{-\kappa(H - X_3)\} - b/a \\ &= 0.2919 \exp\{-0.00027(5000 - X_3)\} + 0.7081. \end{aligned} \quad (3.133) $$

For the reader's convenience we include in Table 3.1 a summary of the essential parameters of our problem and of the numerical values assigned to them in the subsequent computations.

## Weak Formulation and Numerical Solution

The numerical solution of the compaction problem is based on a weak formulation of the differential equations (3.114) and (3.115), use of the evolution equations (3.93) and (3.119), together with equations (3.117) and (3.118), use of equations (3.113) and (3.120) through (3.132), and the selection of a time and a space
<table><caption>TABLE 3.1 ► Parameters of One-Dimensional Compaction Problem</caption><thead><tr><th>Par.</th><th>Definition</th><th>Value (Unit)</th></tr></thead><tbody><tr><td>C</td><td>Solute concentration (mass fraction)</td><td>variable</td></tr><tr><td>C<sub>eq</sub></td><td>Solute equilibrium concentration (equation 3.121)</td><td>function</td></tr><tr><td>d</td><td>Initial grain diameter; d = 2R<sub>0</sub></td><td>10<sup>-4</sup> (m)</td></tr><tr><td>D<sub>3</sub><sup>*</sup></td><td>Referential effective diffusivity (equation 3.116)</td><td>function (m<sup>2</sup>/s)</td></tr><tr><td>D<sub>gb</sub></td><td>Grain boundary diffusivity; D<sub>gb</sub> = 0.1D</td><td>function (m<sup>2</sup>/s)</td></tr><tr><td>D</td><td>Molecular diffusivity (equation 3.124)</td><td>function (m<sup>2</sup>/s)</td></tr><tr><td>D<sub>0</sub></td><td>Reference molecular diffusivity</td><td>9 × 10<sup>-7</sup> (m<sup>2</sup>/s)</td></tr><tr><td>g</td><td>Gravitational acceleration</td><td>9.81 (m/s<sup>2</sup>)</td></tr><tr><td>h(t)</td><td>Current column height (Figure 3.6)</td><td>variable (m)</td></tr><tr><td>H</td><td>Initial column height (Figure 3.6)</td><td>5000 (m)</td></tr><tr><td>k</td><td>Boltzmann's constant</td><td>1.38 × 10<sup>-23</sup> (J/K)</td></tr><tr><td>K<sub>+</sub><sup>sf</sup></td><td>Rate constant for quartz dissolution</td><td>equation 3.122 (m/s)</td></tr><tr><td>K<sub>-</sub><sup>sf</sup></td><td>Rate constant for quartz precipitation</td><td>equation 3.123 (m/s)</td></tr><tr><td>K<sub>gb</sub><sup>g</sup></td><td>Rate constant for intergranular solution</td><td>equation 3.122 (m/s)</td></tr><tr><td>l</td><td>Current fictitious stretch</td><td>variable</td></tr><tr><td>l<sub>0</sub></td><td>Initial fictitious stretch (Figure 3.5)</td><td>variable</td></tr><tr><td>M<sup>s</sup></td><td>Molar mass of quartz; M<sup>s</sup> = ρ<sup>s</sup>V<sup>s</sup></td><td>55 × 10<sup>-3</sup> (kg/mol)</td></tr><tr><td>N<sub>A</sub></td><td>Avogadro's number</td><td>6.022 × 10<sup>23</sup> (1/mol)</td></tr><tr><td>p</td><td>Pore fluid pressure</td><td>variable (Pa)</td></tr><tr><td>Q<sub>3</sub></td><td>X<sub>3</sub>-component of filter velocity</td><td>variable (m/s)</td></tr><tr><td>r</td><td>Normalized grain radius; r(t) = R(t)/R<sub>0</sub></td><td>variable</td></tr><tr><td>T</td><td>Absolute temperature (equation 3.120)</td><td>function (K)</td></tr><tr><td>T<sub>0</sub></td><td>Seafloor temperature</td><td>277.15 (K)</td></tr><tr><td>T'</td><td>Temperature gradient (equation 3.120)</td><td>3 × 10<sup>-2</sup> K/m</td></tr><tr><td>v<sub>s</sub></td><td>Solid volume fraction; v<sub>s</sub> = (1-φ)λ<sub>1</sub>λ<sub>2</sub>λ<sub>3</sub></td><td>variable</td></tr><tr><td>V&#773;<sup>s</sup></td><td>Molar volume of quartz</td><td>2.27 × 10<sup>-5</sup> (m<sup>3</sup>/mol)</td></tr><tr><td>α</td><td>Time integration constant (equation 3.137)</td><td>1/2</td></tr><tr><td>δ</td><td>Grain boundary width</td><td>2 × 10<sup>-9</sup> (m)</td></tr><tr><td>κ</td><td>Depth-decay constant for φ<sub>0</sub> (equation 3.132)</td><td>2.7 × 10<sup>-4</sup> (1/m)</td></tr><tr><td>λ</td><td>Current actual stretch; λ = l/l<sub>0</sub></td><td>variable</td></tr><tr><td>ρ<sup>s</sup></td><td>Constant mass density of solid phase</td><td>2650 (kg/m<sup>3</sup>)</td></tr><tr><td>ρ<sup>f</sup></td><td>Constant mass density of fluid phase</td><td>1000 (kg/m<sup>3</sup>)</td></tr><tr><td>σ'<sub>33</sub></td><td>Effective vertical stress (equation 3.113)</td><td>variable (Pa)</td></tr><tr><td>φ</td><td>Current porosity</td><td>variable</td></tr><tr><td>φ<sub>0</sub></td><td>Initial porosity (equation 3.132)</td><td>function</td></tr><tr><td>φ*</td><td>Initial porosity at layer top</td><td>1-π/6</td></tr><tr><td>Ω<sup>s</sup></td><td>Molecular volume of quartz; Ω<sup>s</sup> = V&#773;<sup>s</sup>/N<sub>A</sub></td><td>3.7 × 10<sup>-29</sup> (m<sup>3</sup>)</td></tr></tbody></table>
discretization. The various parts of this numerical scheme are now discussed. We begin by rewriting equations (3.114) and (3.115) as

$$Q_{3,3} + S_Q(r, \lambda, \dot{r}, \dot{\lambda}) = 0$$

and

$$Q_3 C_{,3} - (D_3^* C_{,3})_{,3} + S_c(r, \lambda, \dot{r}, \dot{\lambda}) + (\lambda - v_s) \dot{C} = 0, \quad (3.134)$$

with source terms $S_Q$ and $S_c$, the definitions and functional dependence of which on $r$, $\lambda$, $\dot{r}$, and $\dot{\lambda}$ will be clear from equations (3.114) and (3.115).

A weak formulation of the problem posed by these two equations is obtained by multiplying the first equation by a filter velocity $\tilde{Q}_3$ and the second by a solute concentration $\tilde{C}$, both of these “virtual” quantities satisfying essentially homogeneous boundary conditions. Next, we integrate the resulting equations over the reference configuration of the compacting layer. The second equation is further integrated by parts so as to eliminate a second-order spatial derivative, and the final results then are

$$\int_0^H \tilde{Q}_3 Q_{3,3} + \tilde{Q}_3 S_Q dX_3 = 0,$$

$$\int_0^H \tilde{C} Q_3 C_{,3} + \tilde{C}_{,3} D_3^* C_{,3} + \tilde{C} (S_c + (\lambda - v_s) \dot{C}) dX_3 = 0. \quad (3.135)$$

A spatial discretization of these equations is now performed, applying the finite-element method. Selecting two-noded Lagrange elements, this leads to a linear interpolation of the filter velocity, the solute concentration, and its rate. After discretization, the system of governing equations assumes the form

$$\int_0^H N_a N_{b,3} dX_3 Q_{3b} = - \int_0^H N_a S_Q dX_3,$$

$$\int_0^H (N_a Q_3 + N_{a,3} D_3^*) N_{b,3} dX_3 C_b + \int_0^H N_a (\lambda - v_s) N_b dX_3 \dot{C}_b = - \int_0^H N_a S_c dX_3, \quad (3.136)$$

in terms of the shape functions $N_a$ and the nodal unknowns $Q_{3b}$, $C_b$, and $\dot{C}_b$. In equation (3.136) and hereafter the summation over nodes with dummy indices $b$ is suppressed for the sake of simplicity.

The following generalized trapezoidal rule is now employed for integrating the solute concentration over time:

$$C_{n+1} = C_n + \Delta t [\alpha \dot{C}_{n+1} + (1 - \alpha) \dot{C}_n]. \quad (3.137)$$
Here the subscripts n and n + 1 designate values at the beginning and at the end, respectively, of the time step $\Delta t$. The parameter $\alpha$ is set to 0, 1/2, or 1, resulting in the selection of the classic forward-Euler, trapezoidal-rule, or backward-Euler time integration scheme, respectively. The system of governing equations finally becomes

$$ \begin{aligned} \dot{\bar{K}}_{ab} Q_{3b(n+1)} &= -\bar{F}_a(r_n, \lambda_n, \dot{r}_n, \dot{\lambda}_n) \\ \bar{K}_{ab}(Q_{3(n+1)}) C_{b(n+1)} + \dot{\bar{K}}_{ab} \dot{C}_{b(n+1)} &= -\bar{F}_a(r_n, \lambda_n, \dot{r}_n, \dot{\lambda}_n). \end{aligned} \quad (3.138) $$

The system must satisfy the initial conditions

$$ \begin{aligned} \dot{\bar{K}}_{ab} Q_{3b0} &= -\bar{F}_a(1, 1, 0, \dot{\lambda}_0) \\ \dot{\bar{C}}_{ab} \dot{C}_{b0} &= -\bar{F}_a(1, 1, 0, \dot{\lambda}_0) - \bar{K}_{ab}(Q_{30}) C_{b0}, \end{aligned} \quad (3.139) $$

where $\dot{r}_0 = 0$ follows directly from equation (3.93) for equation (3.125), while

$$ \dot{\lambda}_0 = \frac{4l_0}{\pi(1 - l_0^2)[\tau_s + (1 - l_0^2)\tau_D]} \left( \frac{\sigma'_{33}(X_3, 0)}{kT/\Omega^s} \right) \quad (3.140) $$

is deduced from equation (3.119) for equations (3.125) through (3.127). The initial effective stress distribution is obtained by integration of equation (3.113) as $\sigma'_{33}(X_3, 0) = -\Delta\rho g(H - X_3)\{1 - \phi^*[1 - \exp\{-\kappa(H - X_3)\}]/\kappa(H - X_3)\}$, noticing that $\sigma'_{33}[h(t)] = 0$ at all times. The characteristic times $\tau_s$ and $\tau_D$ are evaluated in equation (3.140) for the initial temperature profile $T(X_3, 0) = T_0 - T'(H - X_3)$.

The matrix notation used in the preceding should be clear from a comparison of equivalent terms in equations (3.136) and (3.138). The governing equation (3.138) is written for source terms that are evaluated at the beginning of each time step. This approximation was found convenient for obtaining an explicit solution by means of the following staggered five-step scheme:

1. Initialize $Q_{3b}$ and $\dot{C}_b$ with equation (3.139).

2. Solve for $Q_{3b(n+1)}$ from equation (3.138a).

3. Solve for $\dot{C}_{b(n+1)}$ and $C_{b(n+1)}$ from equations (3.138b) and (3.137).

4. Update $r_{n+1} = r_n + \Delta t \dot{r}_n$ and $\lambda_{n+1} = \lambda_n + \Delta t \dot{\lambda}_n$ and recompute $v_s$ from equation (3.117); determine the current position of every nodal point $x_3(X_3, t) = \int_0^{X_3} \lambda(X'_3, t) dX'_3$, deduce the current column height $h(t)$, and transfer this information to the quadrature points by interpolation. Then find the effective stress from $\sigma'_{33}(X_3, t) = -\Delta\rho g \int_{X_3}^H v_s(X'_3, t) dX'_3$, observing that $\sigma'_{33}(H, t) = 0$, and finally determine the temperature from $T(X_3, t) = T_0 + T'[h(t) - x_3(X_3, t)]$ and the equilibrium concentration $C_{eq}$ from equation (3.121).

5. Increment time by $\Delta t$ and proceed to step 2.
Note that the stiffness array in equation (3.138a) remains constant in time and therefore needs to be assembled only once. The solution procedure is thus shortened, because the factorization of the system of equations is not repeated.

Finally, throughout this computation we impose the constraints $\dot{\lambda} \le 0$ and $1 \le r/l_0 \le \sqrt{1+\lambda^2}$ on the variables $r$ and $\lambda$, so as to remain consistent with the assumption of nonexistent *intergranular* precipitation and of circular grain truncations.

The presentation of computed results in the remainder of this section is divided into two parts. The first of these is concerned with the evolution of the sedimentary column in the absence of free-face dissolution or precipitation. In the second part, both intergranular dissolution and free-face dissolution and precipitation processes are allowed to operate. This distinction was found useful in explaining the effectiveness of macroscopic advective and diffusive solute transport as well as the importance of a sink at short distances from the grain boundary to achieving significant compaction rates. The same finite-element mesh is used in both cases. It consists of a set of 150 two-noded elements of nonuniform size, distributed so as to permit decimeter-scale spatial boundary layers at the top and bottom of the layer to be captured. The number of elements was determined by a sensitivity analysis (which is not reported here), to make sure that the numerical results would not depend significantly on the spatial discretization. The time step was given values between $10^{-6}$ years during the early phase of compaction and $10^{-4}$ years in approaching a quasi-stationary state. The small value of the time step is a consequence of our choice of an explicit time integration scheme and also of the large differences between the pertinent characteristic times of the kinetic processes considered (Figure 3.7).

## Compaction in the Absence of Free-Face Precipitation

The results for compaction in the absence of free-face precipitation are presented in Figure 3.8. They were obtained by letting $\tau_I^+$ and $\tau_I^-$ approach infinity. Physically, the suppression of free-face dissolution or precipitation may, for example correspond to the presence of thin clay coatings on the pore walls. The filter velocity profile in Figure 3.8a shows a rapid variation of the velocity over a height of approximately 1500 m, starting from the value zero at the base of the layer. At shallower depths, the filter velocity is approximately constant, signaling the rigid translation of a noncompactant upper portion of the layer.

Clearly, compaction by pressure solution is seen to occur only in the lower part of the layer at the higher temperatures prevailing there. Initially, the rate of compaction is rather high, generating some 2.5 m/yr of displacement at the layer's top; but it then decreases rapidly by a factor of 1000 after a lapse of only one year. This boundary layer in time is depicted more clearly by Figure 3.8b, which shows the rapid decline of the filter velocity at the top of the layer. The vanishing of the filter flux is a clear indication of compaction ceasing to operate, even
<figure><img src="image_8.png" /><figcaption>**Figure 3.7** ► Characteristic times as functions of temperature.</figcaption></figure>

though one has an early burst in dissolution and associated rapid change in the solute concentration profile (Figure. 3.8c). We should note here that the initial solute concentration profile, which according to condition (3.125) equals the equilibrium distribution (3.121), has been altered in the lowermost 10% of the layer by reducing its slope to zero as $X_3 \to 0$, thereby enforcing a boundary condition of zero diffusive flux in agreement with Fick's law (3.112) and condition (3.129). Figure 3.8c shows that the solute concentration then evolves in time toward a second profile that may be interpreted as corresponding to a quasi-stationary state, with a concentration maximum approximately equal to $1.7C_{eq}$ at $X_3 = 0$. This second, quasi-stationary profile, which is well approximated by the dotted curve after $0.19 \times 10^6$ years, establishes itself first at times as short as 0.003 years near the base of the layer, where rapid intergranular dissolution produces a sharp rise in solute concentration. At shallower depths, this evolution toward a quasi-stationary state proceeds at a slower pace and dissolution becomes negligible at about 2000 m, the advective flux $Q_3$ consequently attaining a uniform value. At that point the concentration profile could be further modified only by diffusive solute transport. However, the relevant characteristic time for diffusion over the column height $H$ is approximately $H^2/D_3^*$, corresponding to some $1.5 \times 10^9$ years at a temperature of
<figure><img src="image_9.png" /><figcaption>**Figure 3.8** ► Computed compaction behavior in the absence of free-face solution or precipitation ($\dot{r} = 0$). The filter velocity, solute concentration, and rate of intergranular dissolution profiles are shown in (a), (c), and (d), respectively, for various times. The dotted curve in (c) pertains to $t = 0.2 \times 10^6$ years approximately. Plot (b) and its insert depict the history of the filter velocity at the top of the layer.</figcaption></figure>

$300K$, too large to justify an idea of diffusion-controlled compaction of the entire sedimentary column and thus validating the notion of a quasi-stationary state for the solute concentration found at time $0.19 \times 10^6$ years. We note that the corresponding quasi-stationary concentration profile could have been computed with sufficient accuracy directly from equation (3.119) and the known initial stress and concentration distributions by determining the concentration for which $\lambda$ vanishes, putting $\lambda=1$ as appropriate in the present case. The difference between the initial concentration profile and this quasi-stationary distribution is thus essentially controlled by the effective stress state in the layer.

The computed decline of the rate of intergranular dissolution with time is shown Figure 3.8d. Note that $\lambda$ evolves fastest at the base of the layer, a behavior that is obviously consistent with the first fast evolution of the solute concentration in that
region. Figure 3.8d also provides an indication of the total amount of compaction. The rate of compaction is initially as large as 0.25% per year. However the burst in dissolution only lasts for at most 0.01 years, extending over 1/20 of the column height and leading to an integrated compaction that cannot exceed 2 cm of surface subsidence. The main conclusion to be drawn from these runs must therefore be that in thick layers the absence of sites of deposition close to the source of dissolution will effectively slow down the compaction by IPS to the unrealistically large time scale that is associated with diffusive transport over the entire layer thickness.

## Compaction in the Presence of Free-Face Precipitation

When allowing for free-face solution and precipitation to occur, the computed compaction behavior became extremely unstable. The source of this problem lies in the large difference between the characteristic times $\tau_I^+$ and $\tau_I^-$ for free-face dissolution and precipitation (see Figure 3.7), a slight variation of the solute concentration around its equilibrium value activating either grain growth or dissolution on very different time scales.

From a number of test runs it became nevertheless clear that for the chosen set of parameters free-face dissolution would in fact never occur throughout the time interval covered by the computations. The only characteristic time relevant to the evolution of the grain diameter must therefore be $\tau_I^-$. To alleviate the problem of unstable behavior, the characteristic time for free-face dissolution was therefore set equal to that for precipitation, i.e., we used $\tau_I^+ = \tau_I^-$ as determined by equations (3.94) and (3.122).

In discussing the compaction behavior in the presence of free-face solution and precipitation ($\dot{r} = 0$) we now refer to the results shown in Figure 3.9. The structure of the velocity-depth profiles of Figure 3.9a is similar to the one found for $\dot{r} = 0$. The magnitude of $Q_3$ increases with increasing elevation, rapidly at first, but attaining uniform values at about half the layer thickness. These values are small enough to fully justify our assumption of approximately hydrostatic fluid pressures.

This upper half of the layer is thus subsiding rigidly into the lower half, where elevated temperatures give rise to significant compaction rates. The main difference between the two sets of profiles shown in Figure 3.8a and Figure 3.9a is that the filter velocity does not tend to zero after a short initial burst in dissolution but toward a value decreasing at a rate of approximately $10^{-2}$ m/yr/yr over a period of 30 years (Figure 3.9b). During this period the solute concentration profiles remain close to $C_{eq}$, a result very different from the one obtained in the absence of free-face dissolution and precipitation. Figure 3.9c shows that the ratio

$$ \Delta = \frac{(C - C_{eq})/C_{eq}}{4l_0^2(\Omega^s|\sigma'_{33}|/kT)/[\pi(r^2 - l_0^2\lambda^2)]} $$
<figure><img src="image_10.png" /><figcaption>**Figure 3.9** ► Computed compaction behavior in the presence of free-face precipitation. Depth dependence of filter velocity, relative contribution Δ of chemical term to driving force, grain growth rate, and compaction “stretch” are shown in (a), (c), (d), and (e); evolution of filter velocity at layer surface and of layer height, in (b) and (f). Note small values of Δ in (c) and agreement in (d) of computed $\dot{r}$ (fine solid line) with approximate value (3.144) corresponding to “closed-system” conditions (dashed curve).</figcaption></figure>
of the concentration term over the effective stress term in the driving force in equation (3.119) remains very small indeed. This suggests that the system behavior must be very close to a Stefan-type equilibrium along the grain surfaces, such that the kinetics of free-face precipitation can be disregarded altogether in relation to the rate-determining step of grain boundary dissolution. Moreover, in comparison with local precipitation rates, the rate of diffusive outflux of solute from a unit cell is likely to remain everywhere negligible, suggesting the closed-system approximation $\dot{v}_s \approx 0$ as in the previously discussed examples. From equation (3.118) it is directly apparent that the last approximation implies the following relationship between $\dot{r}$ and $\dot{\lambda}$:

$$ \dot{r} = \frac{-(r^2 - l_0^2 \lambda^2)l_0}{2r[l_0\lambda - 2(r - l_0)]} \dot{\lambda}, \quad (3.141) $$

and when this is combined with relation (3.93) and $(C - C_{eq})/C_{eq}$ is eliminated between the resulting expression and equation (3.119),

$$ \dot{\lambda} = \frac{2rl_0[l_0\lambda - 2(r - l_0)](\sigma'_{33}\Omega^s/kT)}{\frac{\pi}{2}r[l_0\lambda - 2(r - l_0)](r^2 - l_0^2\lambda^2)\tau_s + l_0^2(r^2 - l_0^2\lambda^2)\tau_I^{-} \quad (3.142) $$

Neglecting now the term containing $\tau_I^{-}$ against that containing $\tau_s$, in agreement with the assumption of interfacial equilibrium, we arrive at

$$ \dot{\lambda} = \frac{(4/\pi)l_0}{\tau_s(r^2 - l_0^2\lambda^2)} \left( \frac{\sigma'_{33}}{kT/\Omega^s} \right) \quad (3.143) $$

and by use of equation (3.141)

$$ \dot{r} = \frac{-(2/\pi)l_0^2}{\tau_s[l_0\lambda - 2(r - l_0)]r} \left( \frac{\sigma'_{33}}{kT/\Omega^s} \right). \quad (3.144) $$

Note that because $\dot{v}_s$ is assumed to vanish, the effective stress $\sigma'_{33}$ must remain independent of time and equal to its known initial distribution.

As can be seen from Figure 3.9d, relation (3.144) is indeed satisfied by the computed results. This suggests a very close approximation of the extreme situation, in which compaction by IPS is essentially dissolution-rate controlled and stress driven, behaving as under closed-system conditions (closed w.r.t. to the solid phase) on all time scales significantly below that set by diffusive transport across the entire layer thickness. Because the latter exceeds even geologically meaningful times, IPS would most likely have been predicted to proceed in this manner toward very low porosities if a sufficiently general microstructural evolution had been allowed for in the present model.

We note finally that intergranular dissolution below 3000 m depth is important enough to bring about an integrated compaction of some 40-m surface subsidence (see Figure 3.9f) during the extremely short initial period of only 31 years. At
that time, the constraint $1 \le r/l_0 \le \sqrt{1+\lambda^2}$ is violated at the base of the layer, signaling that the limit of applicability of our microstructural model has been reached.

## 3.5 Concluding Discussion

This chapter has covered the elements of a theory of pressure solution creep up to a level of detail that should facilitate some direct applications as well as further developments by individuals with an interest in this topic. Opportunities for future research in several directions will have become apparent at various places throughout the chapter. These range from an obvious need to further explore the one-dimensional compaction problem, in particular for layers that grow in thickness, to three-dimensional models that might include poroelastic response in addition to pressure solution creep. Further work on more realistic aggregate unit cells and grain truncation geometries will be needed for realistic modeling in the low-porosity range. Complications arising with polymineralic aggregates and more complex dissolution and precipitation kinetics in the presence of clay minerals and micas need to be addressed and incorporated in some tractable manner into a microscale continuum description of unit cell behavior. The kinetics of grain boundary diffusion and dissolution must be given further attention.

The simple one-dimensional compaction model discussed at the end of this chapter clearly represents only a first exploratory study of the behavior of the proposed model. Despite carrying out a numerical study for only one-layer thickness and a single set of initial and boundary conditions, the results we obtained for this case have yielded some basic insights into the behavior of our compaction model that can also point the way for future studies. When precipitation was inhibited, the exceedingly slow diffusive transport over 5000 m of layer thickness was found to lead to the build up of high (possibly unrealistic) supersaturations approximately 1.5 times the equilibrium concentration, which slowed down intergranular dissolution to extremely low rates. However, because this low-rate limit in compaction behavior is likely to be attained only in the absence of any significant advective influx of undersaturated fluid through the base of the layer, the role of the prevailing hydrological regime in compacting sedimentary layers is immediately apparent.

This observation suggests that in a sedimentary column that is open only at its top, macroscopic advection will be activated only if there exist local sinks—the free faces—on the grain scale, where the material at a nearby grain boundary can be deposited. This is confirmed by our numerical results for a sedimentary column in which precipitation is allowed to take place. In this case, the slow diffusive transport over 5000 m of layer thickness leads to another extreme behavior that is characterized by a Stefan-type equilibrium along the grain surfaces, by dissolution rate control, and by stress-driven compaction. Although we could not follow
the evolution toward low enough porosities for want of a sufficiently general microstructural description, IPS may be expected eventually to reach a limit in this case, at which the free-face surface area becomes too small to absorb the dissolved material. Such a development would entail a rise in solute concentration and slow down of the process. However, the corresponding porosity may well be close to zero.

The principal difficulty encountered in our numerical study originates from the fact that the rate processes involved in compaction by IPS operate on vastly different time scales. It was interesting to observe that, with the assumed kinetic data, grain boundary diffusion never played the role of a rate-determining step. However, the characteristic time $\tau_s$ for grain boundary dissolution is in fact not known and our assumption of equal-rate “constants” $K_{sf}^+$ for grain boundary and free-face dissolution remains only a convenient choice at this stage. With the characteristic time for precipitation $\tau_l^-$ being shorter than $\tau_s$, while the latter is still much smaller than that for macroscopic diffusion, our problem evolves on essentially three distinct time scales. More work on suitable numerical (stabilizing) schemes is therefore needed.

## Acknowledgment

The authors are grateful for the support received by F.K. Lehner as CNRS Associate Researcher at the LMS of the École Polytechnique and as visiting professor to the École Normale Supérieure.

## Bibliography

Angevine, C.L., and Turcotte, D.L. Porosity reduction by pressure solution: A theoretical model for quartz arenites. *Geol. Soc. Am. Bull.* 94: 1129–1134, 1983.

Applin, K.R. The diffusion of dissolved silica in dilute aqueous solution. *Geochim. Cosmochim. Acta* 50: 2147–2151, 1987.

Bathurst, R.C.G. Diagenetic fabrics in some British Dinantian limestones. *Liverpool Manchester Geol. J.* 2: 11–36, 1958.

Bear, J. *Dynamics of Fluids in Porous Media*. American Elsevier Publ. Co., Dover Edition. New York, 1988.

Biot, M.A. Theory of finite deformation of porous solids. *Indiana Univ. Math. J.* 21(7): 597–620, 1972.

Bowen, R.M., and Wiese, J.C. Diffusion in mixtures of elastic materials. *Int. J. Engng. Sci.* 7: 689–722, 1976.

Carrio-Schaffhauser, É., Raynaud, S., Latière, H.-J., and Mazerolle, F. Propagation and localization of stylolites in limestones. In *Deformation Mechanisms, Rheology and Tectonics*, Knipe, R.J., and Rutter, E.H., eds., Geol. Soc. Spec. Publ. No. 54, pp. 193–199, London, 1990.
Chadwick, P. *Continuum Mechanics*, 2nd ed. Dover Edition, Mineola, NY, 1999.

Coble, R.L. A model for grain boundary diffusion controlled creep in polycrystalline materials. *J. Appl. Physics* 34: 1679–1682, 1963.

Cox, S.F., and Paterson, M. Experimental dissolution-precipitation creep in quartz aggregates at high temperatures. *Geophys. Res. Lett.* 18: 1401–1404, 1991.

De Groot, S.R., and Mazur, P. *Non-equilibrium Thermodynamics*. North Holland, Amsterdam, 1962.

Dibble, W.E., and Tiller, W.A. Non-equilibrium water/rock interactions I. Model for interface-controlled reactions. *Geochim. Cosmochim. Acta* 45: 79–92, 1981.

Dunnington, H.V. Stylolite development postdates rock induration. *J. Sediment. Petrol.* 24: 27–49, 1954.

Dunnington, H.V. Aspects of diagenesis and shape change in stylolitic limestone reservoirs. In *Proc. Seventh World Petrol. Congress, Mexico*. Vol. 2, pp. 339–352. Elsevier, Amsterdam, 1967.

Durney, D.W. Solution transfer, an important geological deformation mechanism. *Nature* 237: 315–317, 1972.

Durney, D.W. Early theories and hypotheses on pressure-solution-redeposition. *Geology* 6: 369–372, 1978.

Fletscher, R.C. and Pollard, D.D. Anticrack model for pressure solution. *Geology* 9: 419–424, 1981.

Fowler, A.C., and Yang, X.-S. Pressure solution and viscous compaction in sedimentary basins. *J. Geophys. Res.* 104: 12989–12997, 1999.

Gibbs, J.W. On the equilibrium of heterogeneous substances. Reprinted in *The Scientific Papers of J. Willard Gibbs, Vol. I—Thermodynamics*. Dover Edition, New York, 1961.

Glover, G.E. Significance of stylolites in dolomitic limestones. *Nature* 217: 835–836, 1968.

Ghoussoub, J., and Leroy, Y.M. Morphology evolution of solid-fluid interfaces within grain boundaries of porous rock. *J. Mech. Phys. Solids*, 49: 2385–2430, 2001.

Grinfeld, M. Phase transitions of the first kind in nonlinear elastic materials. *Mechanics of Solids* 17: 92–101, 1982. (Translated from *Izv. AN SSSR Mekhanika Tverdogo Tela* 17 (1): 99–109.)

Grinfeld, M. *Thermodynamic Methods in the Theory of Heterogeneous Systems*. Longmans, Green, and Co., Toronto, 1991.

Heidug, W.K. A thermodynamic analysis of the conditions of equilibrium at nonhydrostatically stressed and curved solid-fluid interfaces. *J. Geophys. Res.* 96 (21): 909–21, 921, 1991.

Heidug, W.K., and Lehner, F.K. Thermodynamics of coherent phase transformations in nonhydrostatically stressed solids. *Pure and Applied Geophysics* 123: 91–98, 1985.
Heidug, W.K., and Leroy, Y.M. Geometrical evolution of stressed and curved solid-fluid phase boundaries. Transformation kinetics. *J. Geophys. Res.* 99: 505–515, 1994.

Houseknecht, D.W. Intergranular pressure solution in four quartzose sandstones. *J. Sediment. Petrol.* 58: 228–246, 1988.

Houseknecht, D.W., and Hathon, L.A. Petrographic constraints on models of intergranular pressure solution in quartzose sandstones. *Applied Geochemistry* 2: 507–521, 1987.

Lehner, F.K. On the validity of Fick's Law for transient diffusion through a porous medium. *Chem. Eng. Sci.* 34: 821–825, 1979.

Lehner, F.K. Thermodynamics of rock deformation by pressure solution. In *Deformation Processes in Minerals, Ceramics and Rocks*, Barber, D.J., and Meredith, P.G., eds., pp. 296–333. Unwin Hyman, London, 1990.

Lehner, F.K. A model for intergranular pressure solution in open systems. *Tectonophysics* 245: 153–170, 1995.

Lehner, F.K. Theory of pressure solution creep in wet compacting sediments. In *Mechanics of Solids with Phase Transformations*, Berveiller, M., and Fischer, F.D., eds., pp. 239–258. CISM Courses and Lectures—No. 368. Springer-Verlag, Wien, New York, 1997.

Lehner, F.K., and Bataille, J. Nonequilibrium thermodynamics of pressure solution. *Pure and Applied Geophysics* 122: 53–85, 1984-85.

Lemée, C., and Guéguen, Y. Modeling of porosity loss during compaction and cementation of sandstones. *Geology* 24: 875–878, 1996.

Leroy, Y.M., and Heidug, W.K. Geometrical evolution of stressed and curved solid-fluid phase boundaries: Stability of cylindrical pores. *J. Geophys. Res.* 99: 517–530, 1994.

Machlin, E. S. Some applications of the thermodynamic theory of irreversible processes to physical metallurgy. *J. Metals Trans. AIME*, March 1953, 437–445, 1953.

Mosher, S. Pressure solution as a deformation mechanism in Pennsylvanian conglomerates from Rhode Island. *J. Geology* 84: 355–363, 1976.

Mosher, S. Pressure solution deformation of the Purgatory Conglomerate from Rhode Island. *J. Geology* 89: 37–55, 1981.

Nakashima, S. Diffusivity of ions in pore water as a quantitative basis for rock deformation rate estimates. *Tectonophysics* 245: 185–203, 1995.

Paterson, M.S. Nonhydrostatic thermodynamics and its geologic applications. *Rev. Geophys. Space Phys.* 11: 355–389, 1973.

Raj, R. Creep in polycrystalline aggregates by matter transport through a liquid phase. *J. Geophys. Res.* 87: 4731–4739, 1982.

Rimstigt, J.D. Quartz solubility at low temperatures. *Geochim. Cosmochim. Acta* 61: 2553–2558, 1997.
Rimstidt, J.D., and Barnes, H.L. The kinetics of silica-water reactions. *Geochim Cosmochim. Acta* 44: 1683–1699, 1980.

Rutter, E.H. The kinetics of rock deformation by pressure solution. *Phil. Trans. R. Soc. London* A 283: 203–219, 1976.

Rutter, E.H. Pressure solution in nature, theory and experiment. *J. Geol. Soc. London* 140: 725–740, 1983.

Schutjens, P.M., and Spiers, C.J. Intergranular pressure solution in NaCl: Grain-to-grain contact experiments under the optical microscope. *Oil and Gas Science and Technology-Rev. IFP* 54: 729–750, 1999.

Sclater, J.G., and Christie, P.A.F. Continental stretching: An explanation of the post Mid-Cretaceous subsidence of the North Sea basin. *J. Geophys. Res.* 85: 3711–3739, 1980.

Slattery, J.C. *Momentum Energy and Mass Transfer in Continua*. McGraw-Hill, New York, 1972.

Spiers, C.J., and Brzesowsky, R.H. Deformation behaviour of wet granular salt: Theory versus experiment. In *Seventh Symposium on Salt*. Vol. I. pp. 83–92. Elsevier Sci. Publishers, Amsterdam, 1993.

Spiers, C.J., and Schutjens, P.M. Densification of crystalline aggregates by fluid-phase diffu-sional creep. In *Deformation Processes in Minerals, Ceramics and Rocks*, Barber, D.J., and Meredith, P.G., eds., pp. 334–353. Unwin Hyman, London, 1990.

Spiers, C.J., Schutjens, P.M., Brzesowsky, R.H., Peach, C.J., Liezenberg, J.L., and Zwart, H.J. Experimental determination of constitutive parameters governing creep of rock salt by pres-sure solution. In *Deformation Mechanisms, Rheology and Tectonics* Knipe, R.J., and Rutter, E.H., eds., Spec. Publ. No. 54, pp. 215–227. The Geological Society, London, 1990.

Tada, R., and Siever, R. Pressure solution during diagenesis. *Ann. Rev. Earth Planet. Sci.* 17: 89–118, 1989.

Truskinovskiy, L.M. The chemical-potential tensor. *Geochemistry International* 21: 22–36, 1984. (Translated from *Geokhimiya*, no. 12, pp. 1730–44, 1983).

Wanless, H.R. Limestone response to stress: Pressure solution and dolomitization. *J. Sediment. Petrol.* 49: 437–462, 1979.