Solution Tutorial — 2D Steady Flow in Heterogeneous Aquifers

This tutorial describes the 2D Steady Flow in Heterogeneous Aquifers applet, which simulates steady-state groundwater flow and particle transport in a two-dimensional rectangular domain with a randomly generated heterogeneous hydraulic conductivity field. This applet is an extended version of the original code developed by P. Hsieh, US Geological Survey. The model illustrates how spatial variability in hydraulic conductivity causes macro-scale spreading of fluid particles — analogous to solute dispersion.

For a guide to the applet controls, see the User Interface Tutorial.

1 Introduction

The two-dimensional model ParticleFlow simulates steady-state groundwater flow in a rectangular domain with spatially variable hydraulic conductivity generated by the FGEN random field generator. A key purpose of this model is to illustrate how heterogeneities in hydraulic properties cause the spatial spreading of fluid particles, which is analogous to macro-scale solute dispersion.

Figure 1. Rectangular flow domain. Specified-head boundaries on the left (AD) and right (BC) sides drive average flow in the +x direction. No-flow boundaries on the top (AB) and bottom (DC).

The flow domain is bounded on the left and right sides (AD and BC) by specified hydraulic head boundaries. The top and bottom sides (AB and DC) are no-flow boundaries. With the head along AD greater than the head along BC, the average flow direction is from left to right.

2 Governing Equation

The steady-state groundwater flow equation solved over the domain is:

$$\frac{\partial}{\partial x}\!\left(K\frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y}\!\left(K\frac{\partial h}{\partial y}\right) = 0$$

(1)

where $h$ is hydraulic head, $K$ is hydraulic conductivity (assumed isotropic at each point), and $x$ and $y$ are Cartesian coordinates. The finite element method with a rectangular mesh of square cells — each divided into two triangular elements with linear basis functions — is used to solve Equation (1) for the hydraulic head field.

3 Boundary Conditions

Along the left boundary AD, a constant specified head is imposed:

$$h = h_1$$

(2)

Along the right boundary BC, a second constant specified head is imposed:

$$h = h_2$$

(3)

where $h_1 > h_2$ so that flow proceeds from left to right on average. The top (AB) and bottom (DC) boundaries are no-flow:

$$\frac{\partial h}{\partial y} = 0$$

(4)

4 Velocity Field

After solving for the hydraulic head field, the average interstitial velocity components at each point are computed from Darcy's law divided by porosity:

$$v_x = -\frac{K}{n}\frac{\partial h}{\partial x} \qquad\qquad v_y = -\frac{K}{n}\frac{\partial h}{\partial y}$$

(5)

where $n$ is porosity. These velocity vectors govern the advective movement of fluid particles through the domain. Because $K$ varies spatially, the velocity field is nonuniform even for a simple one-dimensional average head gradient.

Note: Particle paths are purely advective — molecular diffusion and pore-scale dispersion are not included. All spreading observed is due solely to heterogeneity in $K$.

5 Spatial Variance & Macro-Dispersion

In a nonuniform velocity field, a cloud of fluid particles spreads over time. The spreading is quantified by the spatial variance of particle positions in the $x$ and $y$ directions:

$$S_{xx} = \frac{1}{N}\sum_{i=1}^{N}(x_i - x_c)^2 \qquad\qquad S_{yy} = \frac{1}{N}\sum_{i=1}^{N}(y_i - y_c)^2$$

(6)

where $N$ is the total number of fluid particles, $x_i$ and $y_i$ are the coordinates of the $i$-th particle, and $x_c$ and $y_c$ are the coordinates of the center of mass:

$$x_c = \frac{1}{N}\sum_{i=1}^{N} x_i \qquad\qquad y_c = \frac{1}{N}\sum_{i=1}^{N} y_i$$

(7)

If each particle carries a fixed mass of solute, particle spreading is analogous to macro-scale solute dispersion. When the spatial variance grows linearly with time (Fickian behavior), the longitudinal and transverse macro-dispersion coefficients can be estimated from the rate of growth:

$$D_{xx} = \frac{1}{2}\frac{dS_{xx}}{dt} \qquad\qquad D_{yy} = \frac{1}{2}\frac{dS_{yy}}{dt}$$

(8)

Tip: Increase the variance $\sigma_Y^2$ of the log-conductivity field and observe how the rate of particle spreading — and the estimated macro-dispersion coefficients — changes.

6 Symbol Definitions

Symbol	Definition
$h$	Hydraulic head [L]
$h_1$	Specified head on left boundary AD [L]
$h_2$	Specified head on right boundary BC [L]
$K$	Hydraulic conductivity, spatially variable [L/T]
$n$	Porosity [dimensionless]
$v_x,\, v_y$	Average interstitial velocity components in $x$ and $y$ [L/T]
$x,\, y$	Cartesian coordinates [L]
$N$	Total number of fluid particles
$x_i,\, y_i$	Position of the $i$-th particle [L]
$x_c,\, y_c$	Center of mass of the particle cloud [L]
$S_{xx},\, S_{yy}$	Spatial variance of particle positions in $x$ and $y$ [L²]
$D_{xx},\, D_{yy}$	Longitudinal and transverse macro-dispersion coefficients [L²/T]
$\sigma_Y^2$	Variance of the log-hydraulic-conductivity field $Y = \ln K$ [dimensionless]

7 References

Hsieh, P.A. (2001). ParticleFlow: A computer program to delineate pathlines and simulate solute transport in steady-state groundwater flow systems. U.S. Geological Survey Open-File Report 01-286.
Robin, M.J.L., Gutjahr, A.L., Sudicky, E.A., and Wilson, J.L. (1993). Cross-correlated random field generation with the direct Fourier transform method. Water Resources Research, 29(7), 2385–2397.
Gelhar, L.W., and Axness, C.L. (1983). Three-dimensional stochastic analysis of macrodispersion in aquifers. Water Resources Research, 19(1), 161–180.
Dagan, G. (1989). Flow and Transport in Porous Formations. Springer, Berlin.