2.1. Forchheimer Equation

The effective hydraulic conductivity (\(K'\)) in the groundwater flow model derives from the Forchheimer equation for two-regime flow (Darcian and non-Darcian).

For the general case of saturated groundwater flow, where the porous medium is anisotropic and flow is arbitrarily oriented within a rectangular coordinate system, the Forchheimer equation describes the relationship between discharge and hydraulic gradient.

For the simplified one-dimensional case:

(2.2)\[\frac{dh}{dl} = \frac{q}{K} + R \frac{q^2}{K^2}\]

Where:

  • \(q\) = groundwater discharge in the direction of the hydraulic gradient [LT⁻¹],

  • \(l\) = direction of groundwater flow [L],

  • \(K\) = hydraulic conductivity under Darcian flow conditions [LT⁻¹], and

  • \(R\) = ratio of non-Darcian to Darcian flow (a value of 0 indicates flow is completely Darcian) [L⁰].

For the general case of saturated groundwater flow, where the porous medium is anisotropic and flow is arbitrarily oriented within a rectangular coordinate system, the Forchheimer equation can be rewritten as follows:

(2.3)\[q_i = \sum_{j} \left[\frac{1}{\frac{1}{K_{ij}} + \frac{R}{K_{ij}^2} q}\right] \frac{\partial h}{\partial x_j}\]

Where:

  • \(q_i\) = component of groundwater discharge in the i-direction [LT⁻¹],

  • \(K_{ij}\) = hydraulic conductivity tensor for Darcian flow [LT⁻¹],

  • \(|q|\) = magnitude of the discharge vector [LT⁻¹], and

  • \(\frac{\partial h}{\partial x_j}\) = hydraulic gradient component in the j-direction [L⁰].

From the rewritten Forchheimer equation, the effective hydraulic conductivity tensor \(K'\) can be described as:

(2.4)\[K'_{ij} = \frac{1}{\frac{1}{K_{ij}} + b_{ij}q}\]

Where:

  • \(b_{ij}\) = tensor form of b in the Forchheimer equation, where

  • \(b_{ij} = \frac{R}{K_{ij}^2}\) [L²T⁻²]

MINEDW calculates the phreatic surface and potential seepage face (i.e., where the phreatic surface intersects an open-pit highwall) by simulating both saturated and unsaturated flow. Horizontal effective hydraulic conductivity is scaled by relative hydraulic conductivity, which is the ratio of unsaturated hydraulic conductivity to saturated hydraulic conductivity. Note that vertical hydraulic conductivity is not adjusted for unsaturated elements. For the unsaturated zone, this ratio is always less than unity because unsaturated hydraulic conductivity is less than saturated hydraulic conductivity.

For elements above the phreatic surface, the relative hydraulic conductivity equals the user-defined value (\(k_{r}'\)). Below the phreatic surface, the (\(k_r\)) is equal to unity. For each element which comprises the phreatic surface, the ratio of the element’s saturated thickness to its total thickness is calculated and used as the (\(k_r\)) in subsequent calculations if it is greater than the user-defined (\(k_r\)) value. Thus, the pressure head within an element is derived from the pressure head calculated on the nodes of that element.

(2.5)\[\begin{split}k_{r} = \begin{cases} k_{r}', & \text{for } P < 0,\\ 1, & \text{for } P \ge 0 . \end{cases}\end{split}\]
(2.6)\[P = h - z\]

Where:

  • \(P\) = pressure head of node associated with element e [L],

  • \(z\) = vertical elevation coordinate [L], and

  • \(k_{r}’\) = relative hydraulic conductivity in the domain above the phreatic surface [L⁰].