Solution Tutorial

It is well-recognized that natural porous media exhibit significant small-scale spatial variability in hydraulic conductivity. This heterogeneity has been documented in several detailed field experiments, including the Borden and Cape Cod tracer tests:

• Mackay et al., Water Resour. Res., 22(13), 2017–2030, 1986
• LeBlanc et al., Water Resour. Res., 27(5), 895–910, 1991

These studies and subsequent research led to the development of stochastic models for subsurface heterogeneity. The simplest and most widely used model assumes that the hydraulic conductivity K(x) is a realization of a log-normally distributed stationary random field.

Field Decomposition

The spatially variable hydraulic conductivity is expressed as the exponential of a random log-conductivity field Y(x):

where μ_Y is the expected value or mean of Y, and y(x) is a zero-mean deviation from μ_Y. We assume that the mean value is a constant. The vector x represents a spatial location. We only consider a two-dimensional random field here, so x = (X₁, X₂). The coordinate X₁ is assumed to be the horizontal direction and X₂ is the transverse direction.

The fluctuation y(x) is modeled as a stationary Gaussian random field, completely characterized by its covariance function.

Covariance Function

The covariance between the field values at two arbitrary points x₁ and x₂ is defined as:

where E{} refers to the Expectation taken over an ensemble of random fields. Since the mean is assumed constant and y is zero-mean, expression (2) reduces to:

Here we further assume that the random field y(x) is second-order stationary so that the covariance function between two points x₁ and x₂ depends only on the separation distance ξ = x₂ − x₁. Thus R_YY(x₁, x₂) = R_YY(x₂ − x₁) = R_YY(ξ). Note that when the separation distance is zero, the covariance reduces to the variance of y: R_YY(0) = σ_Y².

Two covariance models are available in the applet. Both are anisotropic, with potentially different correlation lengths in x (λ₁) and y (λ₂).

Exponential Covariance

The exponential covariance model produces fields with a relatively rough (non-differentiable) spatial structure. It is commonly used in groundwater modeling:

where:

• σ_y² = variance
• λ₁, λ₂ = correlation lengths in the x₁ and x₂ directions
• ξ₁, ξ₂ = separation distances/components of ξ

Gaussian Covariance

The Gaussian (bell-shaped) covariance model produces fields that are infinitely differentiable (very smooth):

Gaussian fields appear visually smoother than exponential fields with the same λ, because the covariance decays more slowly near the origin.

Nugget Effect

The nugget variance σ_n² adds a discontinuity at the origin, representing small-scale variability below the grid resolution:

where δ(ξ) is the Dirac delta at the origin and Ĉ is the normalized covariance model. A nugget of 0 (default) means all variability is spatially correlated.

Overview

The Fgen92 algorithm uses the direct Fourier transform (spectral) method to generate a realization of the correlated Gaussian field Y(x). The key idea is that a stationary Gaussian random field is completely characterized by its spectral density function S(k₁, k₂), which is the 2D Fourier transform of the covariance function R_YY(ξ).

Spectral Density Function

For the exponential covariance model, the 2D spectral density is:

For the Gaussian covariance model:

where k₁ and k₂ are the angular wavenumbers and H_corv = λ₁ · λ₂ is the correlation volume.

Simulation Steps

The algorithm proceeds on a periodic grid of size nfull_x × nfull_y (powers of 2):

1. Compute spectral amplitudes: For each wavenumber pair (k₁, k₂), compute A(k) = √[S(k₁, k₂) · Δk₁ · Δk₂], where Δk are the wavenumber spacings.
2. Generate random phases: For each mode, draw two independent standard normal random numbers U₁, U₂ using the seeded pseudo-random generator. These form the real and imaginary parts of the complex Fourier coefficient: Z(k) = A(k) · (U₁ + i U₂).
3. Inverse FFT: Apply a 2D inverse FFT to the array of Z(k) values. This produces a realization of the spatially correlated Gaussian field Y(x) at each grid node. The FFT is unnormalized (following the original FORTRAN Fgen92 convention, matching the ISIGN=1 convention of the original FFTMD subroutine).
4. Extract sub-grid: Only the first Nx × Ny values are used (the remainder of the nfull grid is discarded). The mean μ_Y is added to obtain the final log K field: Y(x) = μ_Y + y(x).

Pseudo-Random Number Generator

The applet uses a multi-generator PRNG faithful to the original FORTRAN implementation (RANH/GAUSH routines). The generator produces standard normal deviates using a combination of multiplicative congruential generators with constants M1, M2, M3, ensuring the same field is produced from the same seed across platforms.

FFT Grid Size and Accuracy

The FFT grid (nfull_x, nfull_y) must be powers of 2 and must be at least as large as the desired output grid (Nx, Ny). Larger FFT grids improve the spectral resolution of the simulation (finer wavenumber sampling), which yields a more accurate representation of the target covariance. The minimum power-of-2 grid size is used by default.

Primary Reference (Algorithm)

• Robin, M.J.L., Gutjahr, A.L., Sudicky, E.A., and Wilson, J.L. (1993). “Cross-Correlated Random Field Generation with Direct Fourier Transform Method.” Water Resources Research, 29(7), 2385–2397.

Field Experiment Studies

• Mackay, D.M., Freyberg, D.L., Roberts, P.V., and Cherry, J.A. (1986). “A natural gradient experiment on solute transport in a sand aquifer: 1. Approach and overview of plume movement.” Water Resources Research, 22(13), 2017–2030.
• LeBlanc, D.R., Garabedian, S.P., Hess, K.M., Gelhar, L.W., Quadri, R.D., Stollenwerk, K.G., and Wood, W.W. (1991). “Large-scale natural gradient tracer test in sand and gravel, Cape Cod, Massachusetts, 1. Experimental design and observed tracer movement.” Water Resources Research, 27(5), 895–910.

Further Reading — Stochastic Groundwater Hydrology

• Dagan, G. (1989). Flow and Transport in Porous Formations. Springer-Verlag.
• Gelhar, L.W. (1993). Stochastic Subsurface Hydrology. Prentice-Hall.
• Zhang, D. (2002). Stochastic Methods for Flow in Porous Media: Coping with Uncertainties. Academic Press, San Diego, CA.