Comprehensive Study Notes on Unsaturated Zone Processes and Hydrogeophysics

Introduction to Unsaturated Zone Processes

The vadose zone, representing the part of the subsurface located above the water table, is the site of critical processes controlling mass and energy exchanges between the subsurface and the atmosphere. This zone provides vital boundary conditions for atmospheric processes, including micro-meteorology and climate change, and governs subsurface water migration. The rates, patterns, and timing of aquifer recharge are primarily controlled by percolation through this zone. It acts as a filter where contaminants released near the ground surface may be physically, chemically, or biologically altered, retarded, or removed before reaching aquifers.

Unsaturated processes also dictate water availability for agriculture and serve as driving mechanisms for slope stability, flooding, and various engineering geology challenges. Despite its importance, the vadose zone is often ignored or simplified in hydrologic modeling due to limited data availability and the natural variability of soil properties across different scales.

Principles of Unsaturated Flow in Porous Media

Water flow in unsaturated porous media is driven by differences in hydraulic head $H$ [L], which is the sum of pressure head $\psi$ and elevation head $z$:

$H = \psi(\theta) + z$

Water Retention and the van Genuchten Model

The water retention curve describes the relationship between pressure head and volumetric water content $\theta$ [-]. A widely used empirical four-parameter model for this relationship is the van Genuchten (1980) model:

$\theta = \theta_r + \frac{\theta_s - \theta_r}{(1 + (\alpha |\psi|)^n)^m}$

In this equation:

• $\theta_r$ is the residual water content at asymptotically infinite suction.

• $\theta_s$ is the water content at full saturation (approximately equal to porosity).

• $\alpha$ [L$^{-1}$] is inversely proportional to the mean soil grain size.

• $n$ is an exponent; higher values indicate a wider grain and pore size distribution.

• The model is often simplified by setting $m = 1 - \frac{1}{n}$.

The relationship exhibits hysteresis, meaning the pressure-saturation curve differs during drainage versus imbibition. When a process is reversed mid-way, it follows an intermediate path known as a scanning curve.

The Richards' Equation

Water flow occurs according to the Darcy-Buckingham law:

$q = -K(\psi) \nabla H$

Where $q$ is the water flux [L/T] and $K(\psi)$ is the unsaturated hydraulic conductivity [L/T]. Combining this with the principle of mass conservation leads to the Richards' equation for flow in three dimensions:

$\frac{\partial}{\partial x_i} \left( K(\psi) \frac{\partial H}{\partial x_i} \right) = \frac{\partial \theta(\psi)}{\partial t}$

Where $i = 1, 2, 3$ represents spatial coordinates and $t$ is time.

Hydrogeophysical Petrophysical Relationships

Quantitative hydrogeophysics relies on relating measured electrical properties—electrical resistivity $\rho$ ($\Omega \cdot m$) or conductivity $\sigma$ (S/m) and relative dielectric permittivity $\kappa$—to hydrologic variables like volumetric water content.

Electrical Conductivity Models

For rock samples at full saturation where fluid conductivity $\sigma_w$ dominates, Archie’s Law (1942) applies:

$\sigma_b = \frac{\sigma_w}{F}$

Where the formation factor $F$ is defined by porosity $\phi$ and empirical constants $a$ and $m$ (the cementation exponent, typically $1.0 - 2.5$):

$F = a \phi^{-m}$

If the matrix conductivity $\sigma_s$ is non-negligible, the relationship is generalized as:

$\sigma_b = \frac{\sigma_w}{F} + \sigma_s$

For unsaturated conditions, Archie's Law is extended using water saturation $S_w$ and an empirical exponent $n$ (often close to $2$):

$\sigma_b = \sigma_w F^{-1} S_w^n$

The Waxman and Smits (1968) model for shaly sandstones further accounts for matrix effects:

$\sigma_b = S_w^n F^{-1} \left( \sigma_w + \frac{\sigma_s}{S_w} \right)$

Time-lapse conductivity measurements allow the calculation of saturation changes without knowing $F$:

$\frac{\sigma_b(t_2)}{\sigma_b(t_1)} = \left( \frac{S_w(t_2)}{S_w(t_1)} \right)^n$

Dielectric Constant Models

1. CRIM (Complex Refractive Index Model): A volume averaging model incorporating porosity $\phi$, volumetric water content $\theta$, and constants for solid matrix ($\kappa_s$), air ($\kappa_a$), and water ($\kappa_w$):

$\kappa^{\alpha} = (1 - \phi) \kappa_s^{\alpha} + \theta \kappa_w^{\alpha} + (\phi - \theta) \kappa_a^{\alpha}$

Using $\alpha=1/2$, changes in $\theta$ can be estimated as:

$\theta(t_2) - \theta(t_1) = \frac{\kappa(t_2)^{\alpha} - \kappa(t_1)^{\alpha}}{\kappa_w^{\alpha} - \kappa_a^{\alpha}}$

2. Topp Model (1980): An empirical relationship widely used for agricultural soils that relates permittivity directly to water content:

$\theta = -5.3 \times 10^{-2} + 2.92 \times 10^{-2}\kappa - 5.5 \times 10^{-4}\kappa^2 + 4.3 \times 10^{-6}\kappa^3$

Characterization of the Shallow Vadose Zone

Investigating the top few meters of the subsurface requires techniques with high spatial and temporal resolution. Direct methods like gravimetric sampling are invasive, destructive, and labor-intensive.

Time Domain Reflectometry (TDR)

TDR involves inserting parallel metal rods (waveguides) into the soil. A voltage step is applied, and the travel time of the pulse to the rod ends reveals the propagation velocity, which is used to calculate $\kappa$ and $\theta$.

Example: Borden Site, Ontario In a homogeneous sand site, experiments used pairs of waveguides (lengths 40, 60, 80, 100, 120, and 140 cm) to monitor infiltration from a $3\,m \times 3\,m$ dripline system ($15\,cm \times 15\,cm$ grid). While TDR is effective, interval-differenced water content profiles (calculating $\theta$ for specific depth intervals) showed errors, such as unreasonably low $\theta$ at 90 cm. This underscored the risk of mixing lateral and vertical variations when comparing different sample volumes.

Surface Ground-Penetrating Radar (GPR)

GPR uses electromagnetic wave velocity between transmitters and receivers at the surface to estimate $\theta$.

• WARR (Wide Angle Reflection and Refraction): One antenna is fixed; the other moves.

• CMP (Common Mid Point): Both antennas move from a central point.

Example 1: Grugliasco, Italy A 28 September 2004 irrigation experiment used 200 MHz antennas and 0.2 ns sampling. As the wetting front advanced, a low-velocity layer formed. Critically refracted waves allowed for the estimation of the wet layer thickness and velocity. The estimated $\theta$ values (5% dry, 38% wet) matched TDR data.

Example 2: Montemezzo, Italy In this mountain slope study, a thin soil layer over bedrock acted as a refractive waveguide. 100 MHz antennas recorded dispersive ground waves (lower frequencies traveled faster). Processing data in the frequency-wavenumber ($f-k$) domain produced dispersion curves used to invert for soil thickness ($h$) and velocity ($v$). Results identified soil water content ranging from 0.27 to 0.36 across different seasons.

DC Resistivity Imaging

Electrical resistivity tomography (ERT) uses four-electrode arrays at varying spacings to produce 2-D or 3-D images. Unlike GPR, resistivity can probe deeper without resolution loss, but it requires physical electrode contact and is sensitive to salinity and temperature. Surveys like Vertical Electrical Sounding (VES) assume no lateral variation and are used for sounding curves.

Characterization of the Deep Vadose Zone

As depth increases, surface methods lose resolution. Boreholes provide access for more accurate measurements.

Borehole GPR Configurations

• MOG (Multiple Offset Gather): Extensive 2D velocity tomograms; slow to acquire.

• ZOP (Zero Offset Profile): Transmitter and receiver are moved at equal depths; fast, providing a 1D inter-borehole profile.

• VRP (Vertical Radar Profile): One antenna at surface, one in a single borehole; provides a 1D vertical profile.

Case Study: Sherwood Sandstone, Eggborough, UK Natural infiltration was monitored monthly from 1999 to 2001 using ZOP. Gamma ray logs from 12 boreholes identified layered sandstone sequences. Modeling using the Richards' equation demonstrated that matching GPR-derived $\theta$ required accounting for an averaging window of 2–3 m ($Fresnel zone width$).

Case Study: Trecate, Italy (Oil Spill) VRP monitoring followed a 1994 crude oil blowout. Using 100 MHz antennas, surveys every two weeks identified reflections at 2 m and 6.5 m and tracked water table oscillations of 5–6 m. Occam inversion was used for travel time data, preserving physical limits despite travel time picking errors (approx. 0.5 ns above water table, 4 ns below).

Case Study: Hatfield, UK (Tracer Injection) To determine saturated hydraulic conductivity ($K_s$), 2.1 m$^3$ of NaCl-enriched water ($660\,\mu S/cm$) was injected over 3 days. ZOP radar tracked the vertical movement of the tracer bulb’s center of mass. Matching this movement in a 3D Richards' equation model identified an effective field $K_s$ of $0.4\,m/d$.

Deep Borehole Resistivity

Cross-borehole ERT provides high-resolution images between boreholes. In the Hatfield tracer test, ERT confirmed GPR results but also suggested significant lateral spreading of the tracer. However, a mass balance test on 3-D resistivity images showed a 50% under-prediction of the known injected volume, highlighting sensitivity limitations in the center of the inter-borehole image.