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 a control volume fixed in space, 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 is equal to the integral of the divergence of that field over the volume. In our context, it equates the net outflow through the surface to the local expansion or compression of the quantity 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. In physical terms, it represents the rate at which a quantity (e.g. mass, momentum, fluid) is expanding or contracting at a point. A positive divergence means the point behaves like a source (more is leaving than entering), while a negative divergence indicates a sink.

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 express the time derivative of the 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. \]

Reynolds' theorem accounts for the fact that the control volume may deform or move (e.g. due to solid deformation). It states that the rate of change of a quantity inside a moving/deforming volume equals the local time derivative plus the effect of the volume change (divergence of the velocity field).

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. \]

The total time derivative of any scalar field following the water phase must account for both the local and convective changes. It is expressed as:

\[ \dfrac{\mathrm{d}^\alpha \sqcup}{\mathrm{d} t} = \dfrac{\mathrm{d}^s \sqcup}{\mathrm{d}t} + \bar{\mathbf{w}}^w \cdot \nabla \sqcup. \]

Applying this to the pressure and 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 \varphi^w \cdot \bar{\mathbf{w}}^w + \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 \cdot \nabla \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\), the porosity, we arrive at the final form of the mass balance of pore water:

\[\label{eq:mb-water-final} n (\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 the pore volume due to volumetric deformation, \(n (\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, \(\nabla \cdot \mathbf{w}^w\).