Mass Balance

In a two-phase porous medium, two mass balance equations must be formulated - one for each phase. For a saturated soil, this means one balance of mass for the solid skeleton and one for the water phase. To derive these, we begin with the general integral form of balance laws.

General Form of Balance Equations

Let \(\Omega \subset \mathbb{R}^3\) denote an arbitrary material control volume, with boundary \(\partial \Omega\) and outward unit normal vector \(\mathbf{n}\). The general integral form of a balance equation for an arbitrary physical quantity \(\psi\) is given by:

\[ \dfrac{\mathrm{d}}{\mathrm{d}t} \int_{\Omega} \boldsymbol{\psi} \mathrm{d}v = \int_{\partial \Omega} \boldsymbol{\Psi} \cdot \mathbf{n} \mathrm{d}s + \int_{\Omega} \boldsymbol{\Pi} \mathrm{d}v + \int_{\Omega} \hat{\boldsymbol{\psi}} \mathrm{d}v. \]

where:

\(\boldsymbol{\psi}\) is the physical quantity stored in the domain (e.g. partial mass density),
\(\boldsymbol{\Psi}\) is the flux vector of \(\psi\) across the boundary \(\partial \Omega\) (e.g. mass flux),
\(\boldsymbol{\Pi}\) is an external source (e.g. injection),
\(\hat{\boldsymbol{\psi}}\) is an internal production or consumption term (e.g. phase change, typically zero for water).

This equation states that the rate of change of the quantity inside the control volume equals the net inflow across the boundary, plus external sources, plus internal generation.

In the finite element method (FEM), we discretise the domain into elements and evaluate field variables (displacement, pressure, etc.) locally - either at nodes or integration points. The local form of the balance equation tells us what must hold at every point in the material, not just for an entire control volume. To derive the local form of this balance - which must hold at every point in space - we proceed in two steps.

Application of Gauss's Theorem

First, we convert the surface integral into a volume integral using the divergence theorem (Gauss's theorem):

\[ \int_{\partial \Omega} \boldsymbol{\Psi} \cdot \mathbf{n} \mathrm{d}s = \int_{\Omega} (\nabla\cdot \boldsymbol{\Psi}) \mathrm{d}v. \]

Gauss's theorem states that the total flux of a vector field through the boundary of a control volume equals the integral of its divergence over that volume. In our context, it equates the net outflow of \(\boldsymbol{\Psi}\) across the surface to the integrated local source/sink intensity (i.e. the divergence) inside the volume.

\(\nabla\) operator and divergence

The nabla operator \(\nabla\) is a vector differential operator that acts on scalar, vector, or tensor fields. It provides a compact notation to express spatial derivatives.

One of its key applications is the divergence, which is defined as follows for a vector field \(\mathbf{a} = a_i \mathbf{e}_i\):

\[ \nabla \cdot \mathbf{a} = \dfrac{\partial a_i}{\partial x_i}. \]

The divergence of a vector field measures the net outflow of that field from an infinitesimal control volume, per unit volume. A positive divergence means the point behaves like a source (more is leaving than entering); a negative divergence indicates a sink.

Two specific physical interpretations are useful here:

When the vector field is a flux \(\boldsymbol{\Psi}\) (e.g. mass flux), \(\nabla\cdot\boldsymbol{\Psi}\) is the rate of accumulation or loss of the transported quantity per unit volume.
When the vector field is a velocity \(\mathbf{v}\), \(\nabla\cdot\mathbf{v}\) is the volumetric expansion rate of the material — i.e. the relative rate at which an infinitesimal material volume changes its size.

Both interpretations appear in the derivations below: the first in Gauss's theorem applied to fluxes, the second in Reynolds' transport theorem applied to a material volume.

Similarly, the divergence of a second-order tensor field \(\boldsymbol{\sigma}\) (e.g. the stress tensor) is a vector field:

\[ (\nabla \cdot \boldsymbol{\sigma}) = \dfrac{\partial \sigma_{ij}}{\partial x_j}. \]

Here, the divergence gives the net force per unit volume resulting from spatial variations in the tensor field. In continuum mechanics, this term appears in the momentum balance equation and represents internal forces due to stress gradients.

Substituting into the balance equation:

\[ \dfrac{\mathrm{d}}{\mathrm{d}t} \int_{\Omega} \boldsymbol{\psi} \mathrm{d}v = \int_{\Omega} (\nabla\cdot \boldsymbol{\Psi}) \mathrm{d}v + \int_{\Omega} \boldsymbol{\Pi} \mathrm{d}v + \int_{\Omega} \hat{\boldsymbol{\psi}} \mathrm{d}v. \]

Application of Reynolds Transport Theorem

Next, we apply the Reynolds transport theorem to a material control volume — one that moves with the material at velocity \(\mathbf{v}\). The theorem expresses the time derivative of a volume integral as a pointwise quantity:

\[ \dfrac{\mathrm{d}}{\mathrm{d}t} \int_{\Omega} \boldsymbol{\psi}\, \mathrm{d}v = \int_{\Omega} \left( \dfrac{\mathrm{d} \psi}{\mathrm{d}t} + \psi\, (\nabla\cdot \mathbf{v}) \right) \mathrm{d}v, \]

where \(\dfrac{\mathrm{d} \psi}{\mathrm{d}t} = \dfrac{\partial \psi}{\partial t} + \mathbf{v} \cdot \nabla \psi\) is the material time derivative (the rate of change of \(\psi\) following a material point), and \(\nabla \cdot \mathbf{v}\) is the volumetric expansion rate of the material.

The theorem therefore states that the rate of change of a quantity inside a material volume equals the material time derivative of \(\psi\) plus the contribution arising from the expansion or compression of the volume itself as it follows the material.

In a multi-phase setting, the theorem is applied separately to each phase \(\alpha\), with \(\mathbf{v}\) replaced by the intrinsic velocity \(\bar{\mathbf{v}}^\alpha\) of that phase and \(\dfrac{\mathrm{d}}{\mathrm{d}t}\) understood as the corresponding phase-specific material derivative \(\dfrac{\mathrm{d}^\alpha}{\mathrm{d}t}\).

Substituting this into the balance equation:

\[ \int_{\Omega} \left( \dfrac{\mathrm{d} \psi}{\mathrm{d}t} + \boldsymbol{\psi} (\nabla\cdot \mathbf{v}) \right) \mathrm{d}v = \int_{\Omega} (\nabla\cdot \boldsymbol{\Psi}) \mathrm{d}v + \int_{\Omega} \boldsymbol{\Pi} \mathrm{d}v + \int_{\Omega} \hat{\boldsymbol{\psi}} \mathrm{d}v. \]

Since this equality must hold for any arbitrary control volume \(\Omega\), we obtain the local form of the balance equation:

\[ \dfrac{\mathrm{d} \boldsymbol{\psi}}{\mathrm{d}t} + \boldsymbol{\psi} (\nabla\cdot \mathbf{v}) = (\nabla\cdot \boldsymbol{\Psi}) + \boldsymbol{\Pi} + \hat{\boldsymbol{\psi}}. \]

Mass Balance of the Solid

The mass balance of the solid phase is obtained by replacing the physical field \(\boldsymbol{\psi}\) with the scalar field of the partial solid density \(\rho^s = \varphi^s \bar{\rho}^s\). This quantity describes how the amount of solid material per unit volume of the mixture changes over time. In other words, we track how the solid content within the control volume evolves.

\[\label{eq:mb-solid} \dfrac{\mathrm{d} \rho^s}{\mathrm{d}t} + \rho^s (\nabla \cdot \mathbf{\bar{v}^s}) = 0 \]

From the definition of partial density:

\[\label{eq:partial-density-solid} \rho^s = \varphi^s \bar{\rho}^s = (1 - n) \bar{\rho}^s, \]

we see that changes in \(\rho^s\) may arise from two sources:

A change in the intrinsic density \(\bar{\rho}^s\) (e.g. due to compressibility of the solid grains),
A change in the volume fraction \(\varphi^s\) at constant grain density, i.e. a rearrangement of grains.

Expanding the total derivative of \(\rho^s\) gives:

\[ \dfrac{\mathrm{d} \rho^s}{\mathrm{d}t} = \dfrac{\mathrm{d} (\varphi^s \bar{\rho}^s)}{\mathrm{d}t} = \dfrac{\mathrm{d} ((1 - n)\bar{\rho}^s)}{\mathrm{d}t} \]

For most practical applications, the solid grains are effectively incompressible compared to the bulk skeleton. That is, the intrinsic bulk modulus \(\bar{K}^s\) of the solid grains is much larger than the effective bulk modulus \(K^s\) of the skeleton, so that:

\[\label{eq:mb-solid-1} \dfrac{\mathrm{d} ((1 - n)\bar{\rho}^s)}{\mathrm{d}t} = \dfrac{\mathrm{d} (1 - n)}{\mathrm{d}t} \bar{\rho}^s + (1 - n) \underbrace{\dfrac{\mathrm{d} \bar{\rho}^s}{\mathrm{d}t}}_{= 0} \]

The second term in Equation (\(\ref{eq:mb-solid}\)) represents the effect of volume change due to the deformation of the skeleton. Substituting Equations (\(\ref{eq:partial-density-solid}\)) and (\(\ref{eq:mb-solid-1}\)) into (\(\ref{eq:mb-solid}\)), and dividing through by \(\bar{\rho}^s\), yields:

\[\label{eq:mb-solid-2} \dfrac{\mathrm{d} (1 - n)}{\mathrm{d}t} + (1 - n)(\nabla \cdot \mathbf{\bar{v}^s}) = 0 \]

Rewriting Equation (\(\ref{eq:mb-solid-2}\)) gives:

\[\label{eq:mb-solid-3} \dfrac{\mathrm{d} n}{\mathrm{d}t} = (1 - n)(\nabla \cdot \mathbf{\bar{v}^s}) \]

Recognising that the divergence of the solid velocity corresponds to the volumetric strain rate \(\dfrac{\mathrm{d} \varepsilon^v}{\mathrm{d}t}\), the final form of the mass balance for the solid phase can be written as an evolution equation for the porosity:

\[\label{eq:mb-solid-final} \dfrac{\mathrm{d} n}{\mathrm{d}t} = (1 - n)\dfrac{\mathrm{d} \varepsilon^v}{\mathrm{d}t} \]

This result shows that the porosity increases when the skeleton undergoes volumetric expansion and decreases during compression. It provides an essential link between kinematics and the storage properties of the porous medium.

Mass Balance of Pore Water

Neglecting mass exchange between the solid and fluid phases, and substituting \(\rho^w\) for \(\psi\) in the general balance equation, the mass balance of the pore water phase reads:

\[\label{eq:mb-water} \dfrac{\mathrm{d} \rho^w}{\mathrm{d}t} + \rho^w (\nabla \cdot \mathbf{\bar{v}^w}) = 0. \]

Analogous to the derivation for the solid phase, we substitute \(\rho^w = \varphi^w \bar{\rho}^w\) into Equation \eqref{eq:mb-water}:

\[\label{eq:mb-water-1} \dfrac{\mathrm{d} (\varphi^w \bar{\rho}^w)}{\mathrm{d}t} + \varphi^w \bar{\rho}^w (\nabla \cdot \mathbf{\bar{v}^w}) = 0. \]

Expanding the material derivative yields:

\[\label{eq:mb-water-2} \varphi^w \dfrac{\mathrm{d} \bar{\rho}^w}{\mathrm{d}t} + \bar{\rho}^w \dfrac{\mathrm{d} \varphi^w}{\mathrm{d}t} + \varphi^w \bar{\rho}^w (\nabla \cdot \mathbf{\bar{v}^w}) = 0. \]

In contrast to the solid phase, the pore water cannot always be assumed incompressible. Even small amounts of entrapped air significantly reduce the effective bulk modulus. For example, for an ideally saturated soil (\(S^e = 1\)), the intrinsic bulk modulus of water is \(\bar{K}^w(S^e = 1) \approx 2.2\) GPa, which may be treated as incompressible in practical applications. However, a small air content of only 0.1% reduces the bulk modulus to approximately \(\bar{K}^w(S^e = 0.999) \approx 100\) MPa.

Assuming isothermal conditions, the density of pore water can be expressed as pressure-dependent:

\[\label{eq:compressibility-water} \dfrac{\mathrm{d} \bar{\rho}^w}{\mathrm{d}t} = \dfrac{\bar{\rho}^w}{\bar{K}^w} \dfrac{\mathrm{d} \bar{p}^w}{\mathrm{d}t}, \]

where \(\bar{K}^w\) is the bulk modulus of water and \(\bar{p}^w\) is the intrinsic pore pressure.

Substituting Equation \(\eqref{eq:compressibility-water}\) into \(\eqref{eq:mb-water-2}\) and dividing by \(\bar{\rho}^w\) gives:

\[\label{eq:mb-water-3} \varphi^w \left( \dfrac{1}{\bar{K}^w} \dfrac{\mathrm{d} \bar{p}^w}{\mathrm{d}t} \right) + \dfrac{\mathrm{d} \varphi^w}{\mathrm{d}t} + \varphi^w (\nabla \cdot \mathbf{\bar{v}^w}) = 0. \]

We now recall that the motion of the pore water is described in a modified Eulerian framework, i.e. relative to the motion of the solid. According to the relations introduced in the section Kinematic Relations, the water velocity is given by:

\[ \bar{\mathbf{v}}^w = \bar{\mathbf{v}}^s + \bar{\mathbf{w}}^w, \]

where \(\bar{\mathbf{w}}^w\) is the seepage velocity. Substituting into Equation \eqref{eq:mb-water-3} gives:

\[\label{eq:mb-water-4} \varphi^w \left( \dfrac{1}{\bar{K}^w} \dfrac{\mathrm{d} \bar{p}^w}{\mathrm{d}t} \right) + \dfrac{\mathrm{d} \varphi^w}{\mathrm{d}t} + \varphi^w \left( \nabla \cdot (\bar{\mathbf{v}}^s + \bar{\mathbf{w}}^w) \right) = 0. \]

Applying the relation between the material derivatives following the water and the solid phase introduced in Kinematic Relations — i.e. \(\dfrac{\mathrm{d}^w \sqcup}{\mathrm{d}t} = \dfrac{\mathrm{d}^s \sqcup}{\mathrm{d}t} + \bar{\mathbf{w}}^w \cdot \nabla \sqcup\) — to the pressure and to the volume fraction, Equation \eqref{eq:mb-water-4} becomes:

\[\label{eq:mb-water-5} \varphi^w \left( \dfrac{1}{\bar{K}^w} \left( \dfrac{\mathrm{d}^s \bar{p}^w}{\mathrm{d}t} + \bar{\mathbf{w}}^w \cdot \nabla \bar{p}^w \right) \right) + \left( \dfrac{\mathrm{d}^s \varphi^w}{\mathrm{d}t} + \bar{\mathbf{w}}^w \cdot \nabla \varphi^w \right) + \varphi^w \left( \nabla \cdot \bar{\mathbf{v}}^s + \nabla \cdot \bar{\mathbf{w}}^w \right) = 0. \]

Although this form is suitable for numerical implementation, it can be further simplified. First, we combine terms as follows:

\[\label{eq:mb-water-6} \dfrac{\varphi^w}{\bar{K}^w} \dfrac{\mathrm{d}^s \bar{p}^w}{\mathrm{d}t} + \dfrac{1}{\bar{K}^w} \nabla \bar{p}^w \cdot \mathbf{w}^w + \dfrac{\mathrm{d}^s \varphi^w}{\mathrm{d}t} + \nabla \cdot \mathbf{w}^w + \varphi^w (\nabla \cdot \bar{\mathbf{v}}^s) = 0. \]

To express everything in terms of the Darcy velocity, we note that \(\mathbf{w}^w = \varphi^w \bar{\mathbf{w}}^w\), so that:

\[ \nabla \cdot \mathbf{w}^w = \nabla \cdot (\varphi^w \bar{\mathbf{w}}^w) = \bar{\mathbf{w}}^w (\nabla \cdot \varphi^w) + \varphi^w (\nabla \cdot \bar{\mathbf{w}}^w), \]

and:

\[ \bar{\mathbf{w}}^w \cdot \nabla \bar{p}^w = \dfrac{1}{\varphi^w} \mathbf{w}^w \cdot \nabla \bar{p}^w. \]

Substituting these into Equation \eqref{eq:mb-water-6}, and acknowledging that for a saturated soil \(\varphi^w = n\), gives:

\[ \dfrac{n}{\bar{K}^w} \dfrac{\mathrm{d}^s \bar{p}^w}{\mathrm{d}t} + \dfrac{1}{\bar{K}^w} \nabla \bar{p}^w \cdot \mathbf{w}^w + \dfrac{\mathrm{d}^s n}{\mathrm{d}t} + n (\nabla \cdot \bar{\mathbf{v}}^s) + \nabla \cdot \mathbf{w}^w = 0. \]

The last two terms involving the solid skeleton can be combined using the previously derived solid mass balance, Equation \eqref{eq:mb-solid-final}. With \(\varphi^s = 1-n\) and \(\dfrac{\mathrm{d}^s n}{\mathrm{d}t} = (1-n)(\nabla \cdot \bar{\mathbf{v}}^s)\):

\[ \dfrac{\mathrm{d}^s n}{\mathrm{d}t} + n (\nabla \cdot \bar{\mathbf{v}}^s) = (1-n)(\nabla \cdot \bar{\mathbf{v}}^s) + n (\nabla \cdot \bar{\mathbf{v}}^s) = \nabla \cdot \bar{\mathbf{v}}^s. \]

We arrive at the final form of the mass balance of pore water:

\[\label{eq:mb-water-final} (\nabla \cdot \bar{\mathbf{v}}^s) + \dfrac{n}{\bar{K}^w} \dfrac{\mathrm{d}^s \bar{p}^w}{\mathrm{d}t} + \dfrac{1}{\bar{K}^w} \nabla \bar{p}^w \cdot \mathbf{w}^w + \nabla \cdot \mathbf{w}^w = 0. \]

Equation (\(\ref{eq:mb-water-final}\)) states that a change in pore water mass content can result from:

a change in pore volume due to skeleton deformation, \(\nabla \cdot \bar{\mathbf{v}}^s\),
a change in pressure leading to compression or expansion of the pore water, \(\dfrac{n}{\bar{K}^w} \dfrac{\mathrm{d}^s \bar{p}^w}{\mathrm{d}t}\),
pressure-driven redistribution of pore water, \(\dfrac{1}{\bar{K}^w} \nabla \bar{p}^w \cdot \mathbf{w}^w\),
and the net inflow or outflow of water across the boundary, \(\nabla \cdot \mathbf{w}^w\).