Groundwater and analytic checks — OGR Slip2D documentation

The groundwater engine and several other tools of OGR Slip2D are checked against results that admit no approximation: a closed-form solution, an identity that has to hold exactly, or a requirement of consistency.

Each check has an exact target, or one that can be computed independently (Section 4.1.2). Required is the tolerance written in the test; obtained is the value of the quantity that the test asserts, measured by running the model of the test in version 0.1.235. The file of each test is linked in its section, and each can be run on its own (Section 4.1.7).

4.3.1The checks at a glance#

Table 4.3.1 The checks of this page, what each test requires and what version 0.1.235 gives.

Check	Reference	Required	Obtained
Confined flow	Darcy's law: linear head, $Q = K\,i\,A$	Heads within 10⁻⁹ m; flow within 10⁻⁶	8×10⁻¹⁴ m; 2×10⁻¹⁵
Layered media	Harmonic mean (series), arithmetic mean (parallel)	Flow within 10⁻⁴ and 10⁻³	1×10⁻¹²; 1×10⁻¹⁴
Anisotropy	Patch test with $K_{xy} \ne 0$	Heads within 10⁻⁹ m; velocities within 10⁻⁶	2×10⁻¹⁴ m; 4×10⁻¹⁴
Transient flow	Complementary error function	Within 1 % of the head step	0.11 %
Long time	The steady state	Heads within 0.05 m	3×10⁻⁷ m
Finite-element mesh	Area of the regions	Within 10⁻⁶ m²	0
Support back analysis	Active = passive at $F_t = 1$	Difference below 10⁻⁶	0
Hoek–Brown constants	GSI = 100: $m_b = m_i$ , $s = 1$ , $a = 1/2$	10⁻⁹, 10⁻⁹, 10⁻³	exact
DXF round trip	The shape is kept	Area within 10⁻⁹; vertices within 10⁻⁶ m	0; 4×10⁻¹⁵ m
Rectangular dam	Discharge of Dupuit and Charnyi	Error falls to below 5 % as $k_r$ sharpens	34 % to 1.8 %
Retention curve	Derivative of $\theta(\psi)$	Within 10⁻³	9×10⁻⁹
Rapid drawdown	Morgenstern (1963)	Within 6 % of 1.20; worst at the total drawdown	1.2103 (+0.86 %)

In OGR Slip2D

The checks run in the test suite, and the tests behind this page can be run on their own, from a clone of the repository: python tests/_runner.py seepage_v126 transient_v130 mesh_v125 for the first six rows, and the file named in each section for the others. Several groundwater checks are made twice: on the solver, and through the operations layer ogr_api, whose operations the MCP server publishes as tools (Section 5.3), so that the numbers an agent receives are the verified ones.

Source: tests/ · ogr_fem2d/

4.3.2Confined flow: Darcy's law in one dimension#

In a confined, saturated layer between two edges with prescribed heads and no flow across the other boundaries, the head is linear and the discharge follows the law of Darcy (1856):

\begin{aligned} H(x) &= H_0 + (H_1 - H_0)\,\dfrac{x}{L} \\ Q &= K\,\dfrac{H_0 - H_1}{L}\,h \\ u &= \gamma_w\,(H - y) \end{aligned}

where

$H_0,\ H_1$: total heads prescribed at $x = 0$ and $x = L$ (m)
$K$: conductivity of the layer (m/s)
$h$: thickness of the layer, across the flow (m)
$Q$: discharge per unit length of the section (m²/s)
$u,\ \gamma_w,\ y$: pore pressure (kPa), unit weight of water (kN/m³) and elevation (m)

A linear head belongs to the space of the linear triangles, so the finite-element solution of this problem is exact up to round-off whatever the mesh; any fault in the assembly, in the boundary conditions or in the computation of the flow shows as a difference.

Table 4.3.2 Darcy's law: the models of the two tests and what they require.

	Test of the solver	Test through the operations layer
Model	Rectangle L = 20 m by h = 10 m; K = 1×10⁻⁵ m/s; H₀ = 10 m, H₁ = 4 m	Rectangle 20 m by 10 m; K = 1×10⁻⁵ m/s; H₀ = 15 m, H₁ = 12 m
Mesh	Linear triangles, target size 1.2 m	Target of 400 elements
Closed form	Q = 3×10⁻⁵ m²/s	Q = 1.5×10⁻⁵ m²/s
Required	Head within 10⁻⁹ m at every node; discharge through the middle section within 10⁻⁶; the same discharge, within 10⁻⁶, through the sections x = 2, 6, 10, 14 and 18 m; pressure head equal to H − y to 10⁻¹² m	Head within 10⁻⁶ m at four points and u within 10⁻⁵ kPa; discharge through three sections within 10⁻⁶; the sign of the flow follows the orientation of the section
Obtained	Head 7.9×10⁻¹⁴ m; discharge 2.0×10⁻¹⁵; sections 7.3×10⁻¹⁴	Passes

In OGR Slip2D

The solver imposes the prescribed heads by elimination and computes the discharge through any segment of the model (Section 3.7.3). The second test builds the same box through the operations mesh_generate, seepage_bc_set and groundwater_run, and shows, besides, that the conductivity set by hydraulic_set is the one the solver reads.

Source: tests/test_seepage_v126.py · tests/test_groundwater_ops_v1200.py

4.3.3Layered media: harmonic and arithmetic means#

Two layers of conductivities $K_1$ and $K_2$ give two limiting cases, both consequences of Darcy's law and of continuity. When the flow crosses the layers, in series, the same discharge passes through each of them and the head losses add up, which gives the harmonic mean of the conductivities; when it runs along them, in parallel, every layer sees the same head and the discharges add up, which gives the arithmetic mean:

\begin{aligned} K_{\mathrm{eq}} &= \dfrac{L}{\sum_i L_i/K_i} \\ H_{\mathrm{int}} &= H_0 - (H_0 - H_1)\,\dfrac{r_1}{r_1 + r_2}, \qquad r_i = \dfrac{L_i}{K_i} \end{aligned}

\begin{aligned} K_{\mathrm{eq}} &= \dfrac{\sum_i K_i\,b_i}{\sum_i b_i} \\ Q &= K_{\mathrm{eq}}\,\dfrac{H_0 - H_1}{L}\,h \end{aligned}

where

$L_i$: length of layer $i$ along the flow, with $L = \sum L_i$ (m)
$b_i$: thickness of layer $i$ across the flow, with $h = \sum b_i$ (m)
$r_i$: hydraulic resistance of layer $i$ (s)
$H_{\mathrm{int}}$: head at the interface of two layers in series (m)

The check uses a contrast of 100, which puts the numerical method under stress: almost the whole head loss has to take place in the layer of low conductivity, and the flow has to be conserved across the interface.

Table 4.3.3 Layered media: the models of the two tests and what they require.

	Test of the solver	Test through the operations layer
Model	Box 20 m by 10 m; K₁ = 1×10⁻⁴ and K₂ = 1×10⁻⁶ m/s; H₀ = 10 m, H₁ = 4 m; mesh of 0.9 m	Box 20 m by 10 m; K = 1×10⁻⁵ and 4×10⁻⁵ m/s; H₀ = 15 m, H₁ = 12 m
In series	Two layers 10 m wide, side by side; discharge through x = 5 m; Q = 5.94×10⁻⁶ m²/s	Two layers 10 m wide; discharge through x = 5 m, and head at the interface
In parallel	Two layers 5 m thick, one above the other; discharge through x = 10 m; Q = 1.515×10⁻⁴ m²/s	Two layers 5 m thick; discharge through x = 10 m
Required	Discharge within 10⁻⁴ in series and 10⁻³ in parallel; the flow through six sections on both sides of the interface equal within 10⁻⁴; the head loss in the layer of low conductivity more than 90 times that in the other (the exact ratio is K₁/K₂ = 100)	Discharge within 10⁻⁶, in series and in parallel; head at the interface within 10⁻⁶ m of H_int
Obtained	Series 1.0×10⁻¹²; parallel 1.2×10⁻¹⁴; sections 1.3×10⁻¹²; ratio of head losses 100.0000000000	Passes

Figure 4.3.1 Head along the middle of the series model of the test of the solver, computed with OGR Slip2D 0.1.235 (dots, at 0.5 m). The layer of conductivity

10^{-6}

m/s takes 5.94 m of the 6 m of head loss; the first 10 m lose 0.06 m.

In OGR Slip2D

Each layer is a material region with its own hydraulic properties, and the mesh is built region by region, so that the elements do not cross the interface and its nodes are shared (Section 3.7.4).

Source: tests/test_seepage_v126.py · tests/test_groundwater_ops_v1200.py

4.3.4Anisotropy: the patch test with a full tensor#

The conductivity of a material is a tensor whose principal values $K_1$ and $K_2 = (K_2/K_1)\,K_1$ are rotated by the angle $\theta_1$ of $K_1$ from the $+x$ axis, which gives $K_{xy} = (K_1 - K_2)\sin\theta_1\cos\theta_1$ (Eq. (3.7.5)). A linear field of head satisfies the steady equation for any constant tensor, since its gradient is constant:

\begin{aligned} H &= A\,x + B\,y + C \\ \mathbf{v} &= -\mathbf{K}\,\nabla H = -\begin{pmatrix} K_{xx}A + K_{xy}B \\ K_{xy}A + K_{yy}B \end{pmatrix} \end{aligned}

If that field is prescribed on the whole boundary, the finite-element solution has to reproduce it at every node, and the Darcy velocity has to be the constant vector of Eq. (4.3.4) in every element. It is the patch test of an element. With $K_{xy} \ne 0$ it exercises the rotation of the tensor and the cross terms of the element matrix, which isotropic cases cannot, and so it checks both the rotation and the Galerkin assembly.

Table 4.3.4 Patch test: the five tensors of the test. The model is a rectangle 20 m by 10 m, with a mesh of 1.2 m, and

A = 0.3

B = 0.2

and

C = 5.0

Principal conductivity	Ratio	Angle	Cross term
K₁ = 1×10⁻⁵ m/s	K₂/K₁ = 1	0°	$K_{xy} = 0$
K₁ = 1×10⁻⁵ m/s	K₂/K₁ = 0.25	0°	$K_{xy} = 0$
K₁ = 1×10⁻⁵ m/s	K₂/K₁ = 0.25	45°	$K_{xy} \ne 0$
K₁ = 2×10⁻⁴ m/s	K₂/K₁ = 0.1	30°	$K_{xy} \ne 0$
K₁ = 5×10⁻⁶ m/s	K₂/K₁ = 4	−60°	$K_{xy} \ne 0$

Required. The head within 10⁻⁹ m at every node, for the five tensors; the velocity of every element within 10⁻⁶ of the exact one, relative to its larger component; and, for the isotropic case, the magnitude of the gradient equal to $\sqrt{A^2 + B^2}$ within 10⁻⁹. Obtained. Heads within 2.5×10⁻¹⁴ m and velocities within 3.6×10⁻¹⁴, for all five.

In OGR Slip2D

The tensor is built from Saturated permeability Ks, K2 / K1 and K1 angle (deg from +X) of the material, and in the unsaturated zone the whole tensor is multiplied by $k_r$ , so that the ratio and the direction are kept (Section 3.7.2). The tests also check that the tensor is positive definite and that its trace does not change with the rotation.

Source: tests/test_seepage_v126.py · ogr_core/hydraulic/hydraulic_properties.py

4.3.5Transient flow: the response to a step in head#

In a confined, fully saturated aquifer the transient equation reduces to the diffusion equation, with the hydraulic diffusivity $D = K/S_s$ , where $S_s$ is the specific storage. For a semi-infinite aquifer at the uniform head $H_0$ , whose boundary at $x = 0$ is raised to $H_1$ at $t = 0$ , the solution is the classical profile of the complementary error function (Carslaw & Jaeger, 1959):

\begin{aligned} \dfrac{\partial H}{\partial t} &= D\,\dfrac{\partial^2 H}{\partial x^2} \\ H(x,t) &= H_0 + (H_1 - H_0)\,\operatorname{erfc}\!\left(\dfrac{x}{2\sqrt{D\,t}}\right) \end{aligned}

The front advances as $2\sqrt{D\,t}$ . The solution holds only while the front is far from the other end of the aquifer, and it is the analytical solution, not the solver, that stops applying beyond that: the test compares only where the semi-infinite assumption holds.

Table 4.3.5 Step response of a confined aquifer: model and requirement of the test of the solver.

Item	Data
Model	Strip 50 m by 5 m; $K = 1\times10^{-4}$ m/s; $S_s = 1\times10^{-3}$ m⁻¹, so that $D = 0.1$ m²/s; no flow across the long edges
Heads	Initially $H_0 = 100$ m; the left edge is raised to $H_1 = 110$ m, the right edge stays at 100 m (high enough to stay saturated everywhere)
Numerics	Mesh of 1.0 m; 150 time steps; Picard iteration with a tolerance of 10⁻⁹ and at most 12 iterations; times $t = 200$ and $t = 1000$ s
Compared nodes	At mid-height (within 0.6 m of it), for $0.5\ \mathrm{m} \lt x \lt \min(2\sqrt{Dt},\ L - 3\sqrt{Dt})$ : up to 8.94 m at $t = 200$ (16 nodes) and up to 20 m at $t = 1000$ (38 nodes)
Required	The largest difference from the closed form below 1 % of the head step $H_1 - H_0$ , that is, 0.1 m
Obtained	0.109 % at $t = 200$ ; 0.097 % at $t = 1000$

Figure 4.3.2 Head at mid-height of the aquifer after the step, computed with OGR Slip2D 0.1.235 (dots) and the complementary error function (lines). The dashed lines mark the end of the compared range at each time.

The same file checks that the front advances with time and that a storage ten times larger delays it. The same aquifer is built through the operations layer, with the same bound of 1 % at the same two times, and a third time of 10⁶ s has to reproduce the linear steady profile within 10⁻³ of the head step.

In OGR Slip2D

The transient solver integrates the mixed form of Richards' equation in time with a modified Picard iteration, by stages (Section 3.7.8). In a saturated zone the storage is the elastic specific storage $S_s$ , which is what keeps the system from collapsing onto the steady state.

Source: tests/test_transient_v130.py · tests/test_groundwater_ops_v1200.py

4.3.6Long time: the steady state#

As time goes on the stored water stops changing, and the transient solution has to approach the steady solution of the same boundary conditions, which is that of Eq. (3.7.4):

\lim_{t\to\infty} H(\mathbf{x}, t) = H_\infty(\mathbf{x}), \qquad \nabla\cdot\left[\,k_r(\psi)\,\mathbf{K}\,\nabla H_\infty\,\right] = 0

The unsaturated steady problem has no closed form, so this is a check of consistency and not of a closed form: it ties the transient scheme to the steady solver, which is itself checked against the exact discharge of a rectangular dam (Section 4.3.11).

Model. A rectangle 20 m by 12 m, with a mesh of 0.9 m. The permeability function is Gardner's, with $K_s = 1\times10^{-5}$ m/s, $a = 1$ and $n = 3$ ; the retention curve is van Genuchten's, with $\alpha = 0.1$ m⁻¹ and $n = 2$ , $\theta_s = 0.4$ and $\theta_r = 0.05$ ; the specific storage is $1\times10^{-4}$ m⁻¹. At the start the head is 10 m on both sides; the final conditions are 10 m on the left and 3 m on the right, prescribed up to those levels and free above them. The solution after a single stage of 10⁹ s is compared with the steady solution of the final conditions, computed by the unsaturated steady solver (relaxation 0.4, at most 200 iterations, tolerance 10⁻⁶).

Required. The largest difference of total head between the two solutions below 0.05 m. Obtained. 2.8×10⁻⁷ m. Through the operations layer, the confined aquifer of the previous section at $t = 10^6$ s has to lie within 10⁻³ of the head step, 0.01 m, of its linear steady profile. A further test requires that draining lowers the stored water monotonically.

Test: tests/test_transient_v130.py, TestAsymptoticConsistency.

4.3.7The finite-element mesh keeps the area#

A conforming mesh of linear triangles that covers a region exactly has to add up to the area of the region:

\begin{aligned} \sum_e A_e &= A_\Omega \\ A_e &= \tfrac12\left|(x_2 - x_1)(y_3 - y_1) - (x_3 - x_1)(y_2 - y_1)\right| \end{aligned}

The identity fails if the mesh has gaps, overlaps or triangles that leave the domain, so none of these faults can pass it. It is checked on three models:

a square of 10 m by 10 m, with a target element size of 1.5 m: 100 m²;
two squares stacked and sharing an edge, which exercise the conformity of the interface between two regions: 200 m²;
a slope of the test suite with three material regions, at target numbers of 200, 800 and 2 000 elements: the total area must be the same at every refinement and equal to the area of the regions as the program computes it.

Required. A difference below 10⁻⁶ m² in each case, and a spread of the areas below 10⁻⁶ m² across the three refinements. Obtained. No difference at floating-point precision in any of them. The same file checks that no edge is shared by more than two elements, that no element crosses the interface between regions and that the nodes of a shared edge are shared, that every element has a positive area and is numbered anticlockwise, that the smallest angle stays above 20° on the square and above 15° on the slope, and that the shape functions reproduce a linear field exactly: the gradient of $3x - 2y + 7$ comes out right within 10⁻⁶.

In OGR Slip2D

The mesh is built region by region, and the nodes of an edge shared by two regions are shared by both, so the mesh conforms across interfaces by construction (Section 3.7.4).

Source: tests/test_mesh_v125.py · ogr_fem2d/mesh/mesh.py

4.3.8Support back analysis: active and passive forces coincide at a factor of 1#

The back analysis fixes the factor of safety at a target $F_t$ and solves for the horizontal force that has to be added to reach it (Section 3.10.6). With $R$ and $D$ the resisting and the driving sums of the method, evaluated at $F_t$ , and $a$ the arm factor of that section,

\begin{aligned} T_{\mathrm{act}} &= \dfrac{1}{a}\left(D - \dfrac{R}{F_t}\right) \\ T_{\mathrm{pas}} &= \dfrac{1}{a}\left(F_t\,D - R\right) = F_t\,T_{\mathrm{act}} \end{aligned}

so that the two assumptions give the same force at $F_t = 1$ , whatever the slope, the surface and the method, and the passive force is the larger one for any target above 1 when a force is needed. The identity would break under almost any algebraic slip, because $D - R/1$ and $1\cdot D - R$ are the same expression. The force is zero when $F_t$ is the factor of safety of the surface without support.

Models. In the test of the engine, a slope of the test suite with three material regions, the circle with centre (88.0, 70.5) m and radius 47.212 m, 25 slices, Bishop's simplified method and the force applied at an elevation of 25 m. In the test through the operations layer, the operation back_analysis_run on a slope 10 m high at 1V:2H ( $c' = 6$ kPa, $\phi' = 22^\circ$ , $\gamma = 19$ kN/m³) with a grid search of 3 × 3 intervals, a target of 1.0 and the force at an elevation of 5 m.

Required. At $F_t = 1.0$ , a difference between the active and the passive force below 10⁻⁶ (in the units of force per unit length), in both tests. In the test of the engine, also: a passive force above the active one for targets of 1.2, 1.5 and 2.0; a force that grows with the target (1.0, 1.2, 1.4, 1.6); and a force of zero, not negative, for a target the surface already exceeds, with less than 1 kN/m at the unsupported factor itself. Obtained. A difference of 0 at $F_t = 1$ , and a ratio $T_{\mathrm{pas}}/T_{\mathrm{act}}$ of 1.2, 1.5 and 2.0 for those targets, to within 4×10⁻¹⁶.

In OGR Slip2D

The back analysis supports Bishop's simplified method, on circular surfaces, and Janbu's simplified and corrected methods; every other method is refused. The elevation of the force changes the result of Bishop's method, whose moments depend on the arm, and does not change that of Janbu's, whose force balance does not, and both properties are tested. The calculation does not use or alter the results of the stability analysis.

Source: tests/test_back_analysis_v140.py · tests/test_f3b_ops_v1201.py · ogr_slip2d/back_analysis.py

4.3.9Hoek–Brown constants for intact rock#

The Parameter Calculator derives the constants of the generalised Hoek–Brown criterion from the Geological Strength Index, the intact-rock constant $m_i$ and the disturbance factor $D$ (Hoek, Carranza-Torres & Corkum, 2002; Eq. (3.5.14) and Eq. (3.5.15)):

\begin{aligned} m_b &= m_i\exp\!\left(\dfrac{\mathrm{GSI} - 100}{28 - 14D}\right) \\ s &= \exp\!\left(\dfrac{\mathrm{GSI} - 100}{9 - 3D}\right) \\ a &= \dfrac12 + \dfrac16\left(e^{-\mathrm{GSI}/15} - e^{-20/3}\right) \end{aligned}

For $\mathrm{GSI} = 100$ and $D = 0$ the exponents of $m_b$ and $s$ vanish, and the two exponentials of $a$ are the same number, so the rock mass reduces to the intact rock: $m_b = m_i$ , $s = 1$ and $a = 1/2$ , exactly.

Required. For $\mathrm{GSI} = 100$ , $m_i = 10$ and $D = 0$ : $m_b$ within 10⁻⁹ of 10, $s$ within 10⁻⁹ of 1 and $a$ within 10⁻³ of 0.5. Obtained. Exactly 10, 1 and 0.5. A second test uses a value of Hoek et al. (2002): for $\mathrm{GSI} = 50$ , $m_i = 10$ and $D = 0$ the constants are $m_b = 1.6768$ within 10⁻³, $s = 0.003866$ within 10⁻⁵ and $a = 0.5057$ within 10⁻³, and the program gives 1.67677, 0.0038659 and 0.50573. A further test goes through the button of the material dialog and checks, in closed form against $\exp[(\mathrm{GSI} - 100)/(9 - 3D)]$ , the value of $s$ that arrives in the material, and not only the value the calculator shows; a low GSI, whose $s$ is very small, has to arrive whole.

In OGR Slip2D

The calculator clamps GSI to 0–100, and $D$ to 1.999, where $m_b$ would become undefined, and says so in a note instead of returning an infinite constant; both clamps are tested (Section 2.3.2).

Source: tests/test_m6_v157.py · tests/test_gsi_calculator_write_v1192.py · ogr_core/materials/parameter_calculator.py

4.3.10DXF round trip#

A model exported to DXF and imported again has to come back as the same geometry. Stated precisely, the drawing is not vertex for vertex identical, and it should not be: the importer splits the boundaries at their crossings, so the external boundary gains nodes where a material boundary or a water line meets it. What has to hold is that the shape is unchanged:

the enclosed area is the same, with $A = \tfrac12\left|\sum_i (x_i\,y_{i+1} - x_{i+1}\,y_i)\right|$ ;
every original vertex is still present;
every added vertex lies on an original segment.

That is stronger than asking for identical lists, because it still catches a coordinate that is moved, dropped or rounded.

Model. A slope of the test suite whose external boundary has eight vertices, with two material boundaries and a water table from (0, 32) to (120, 22), exported in metres and imported with a welding tolerance of 0.05 % and no simplification; the area is also checked after a round trip in millimetres. Required. The same boundary types and counts; a relative difference of the enclosed area below 10⁻⁹; every original vertex present, matched to 10⁻⁶; every added vertex within 10⁻⁶ m of an original segment; the material boundaries unchanged to nine decimals; the closed or open character of each boundary kept; and the regions still buildable. Obtained. A difference of area of 0, in metres and in millimetres, and an external boundary that goes from 8 to 13 vertices, the five added vertices at up to 4×10⁻¹⁵ m from the original outline.

In OGR Slip2D

The geometry goes to layers that the importer recognises (OGR_EXTERNAL, OGR_MATERIAL, OGR_WATER_TABLE and the others) and the drawings of results go to layers whose name starts with OGR_X_, which the importer ignores, so that re-importing cannot turn a load arrow into a material boundary; a test checks the layer contract and the units recorded in the header (Section 2.2.5).

Source: tests/test_dxf_export_v148.py · ogr_core/dxf/

4.3.11A rectangular dam: the discharge of Dupuit and Charnyi#

For flow through a rectangular dam of width $L$ on an impervious base, with an upstream level $H_1$ and a downstream level $H_2$ above the base, the discharge per unit length is the following expression (Charnyi's result, which is also the discharge of Dupuit's approximation), provided that the zone above the free surface does not conduct:

q = \dfrac{K\,(H_1^2 - H_2^2)}{2L}

The unsaturated formulation of OGR Slip2D lets that zone conduct, so the computed discharge is expected to be higher. Making the permeability function sharper has to drive it, monotonically, down to the closed form; the test checks that convergence, which validates the solver and accounts for the excess.

Model. A dam 20 m wide and 12 m high on an impervious base, $K = 1\times10^{-5}$ m/s, $H_1 = 10$ m on the left face and $H_2 = 2$ m on the right; the crest and the faces above the water levels carry the condition of a possible seepage face; mesh of 0.6 m; the discharge through the vertical section in the middle of the dam. The closed form is $q = 10^{-5}(10^2 - 2^2)/(2\cdot 20) = 2.4\times10^{-5}$ m²/s. The permeability function is that of Gardner (1958), with a floor of $k_r = 10^{-8}$ , in four increasingly sharp versions, listed in Table 4.3.6.

Table 4.3.6 Rectangular dam: discharge for four permeability functions, computed with OGR Slip2D 0.1.235, against the closed form

2.4\times10^{-5}

m²/s.

Gardner parameters	Discharge	Difference
a = 0.05, n = 2	3.213×10⁻⁵ m²/s	+33.9 %
a = 1, n = 3	2.758×10⁻⁵ m²/s	+14.9 %
a = 50, n = 4	2.533×10⁻⁵ m²/s	+5.5 %
a = 10⁴, n = 5	2.443×10⁻⁵ m²/s	+1.8 %

Required. The difference decreases as the function sharpens (a rise of up to 10⁻³ between consecutive cases is tolerated) and the sharpest case falls below 5 %. The same file requires that the free surface leaves the downstream face above the tail-water level and below the upstream one, that it starts near the upstream head, that it descends monotonically, allowing for the noise of the discrete trace, and that the nodes of the seepage face are reported. Through the operations layer the same dam, with the sharpest function, has to give the closed form within 5 %.

In OGR Slip2D

The free surface is the line $P = 0$ of the solution, and the seepage face is found by nodal switching, with a hysteresis band and a budget of switches per node so that the iteration does not chatter (Section 3.7.6).

Source: tests/test_unsaturated_v127.py · tests/test_groundwater_ops_v1200.py

4.3.12The retention curve and its derivative#

The retention curve of van Genuchten (1980) gives the water content as a function of the suction head $\psi$ , and its derivative is the specific moisture capacity that appears as the storage term of the unsaturated zone (Eq. (3.7.15)):

\theta(\psi) = \theta_r + (\theta_s - \theta_r)\left[1 + (\alpha\psi)^n\right]^{-m}, \qquad m = 1 - \dfrac1n

C = -\dfrac{\mathrm{d}\theta}{\mathrm{d}\psi} = (\theta_s - \theta_r)\,m\,n\,\alpha\,(\alpha\psi)^{n-1}\left[1 + (\alpha\psi)^n\right]^{-m-1}

Both hold for $\psi \gt 0$ ; where the soil is saturated, $\theta = \theta_s$ and $C = 0$ . Model. $\theta_s = 0.45$ , $\theta_r = 0.08$ , $\alpha = 0.05$ m⁻¹ and $n = 2$ . Required. $\theta = 0.45$ within 10⁻¹² at zero suction and for positive pore pressure, and below 0.09 at a suction of 10⁶ m; and the analytical capacity within 10⁻³, in relative terms, of the central difference of the curve (relative step of 10⁻⁴) at $\psi = 2,\ 10,\ 40$ and 150 m. Obtained. Differences of 1.6×10⁻¹⁰, 2.0×10⁻⁹, 4.0×10⁻⁹ and 9.4×10⁻⁹, and $\theta = 0.08001$ at 10⁶ m. The same file checks that the storage content is the integral of the storage in both zones, to within 1 %, which is what keeps the modified Picard scheme valid across the water table.

In OGR Slip2D

The retention curve of van Genuchten is used with every permeability model: the permeability function governs how fast the water moves, and the retention curve how much of it is stored (Section 3.7.7).

Source: tests/test_transient_v130.py

4.3.13Rapid drawdown: the factor of safety against level and time#

Two properties of rapid drawdown are checked, one against a published result and one as a physical requirement.

The complete drawdown of a homogeneous slope. For a homogeneous slope of the kind charted by Morgenstern (1963), the reference values used by the test suite put the lowest factor of safety at the complete drawdown, 1.20, against 1.41 at half drawdown. The sweep of drawdown levels repeats the whole search at each level (Section 3.8.6). Model. A slope of 3 horizontal to 1 vertical, with $c'/\gamma H = 0.025$ , $\phi' = 30^\circ$ and $\gamma/\gamma_w = 2$ , undrained behaviour with $\bar B = 1$ , and a full reservoir up to the crest; Bishop's method, a grid of 3 × 4 intervals, 18 slices and four levels from the crest to the toe. Required. The critical level is the complete drawdown; its factor of safety is within 6 % of 1.20; the factor at the intermediate level closest to half drawdown is higher; and sweeping the levels reveals no unsafe margin above 1 %. Obtained. 1.2103 at the complete drawdown, 0.86 % above 1.20, and 1.304, 1.673 and 2.588 at one third, two thirds and the full height of the reservoir.

The factor of safety in time. With the transient analysis, the factor of safety of the stages flagged Calculate SF is lowest right after the level drops, when the pressures have not yet dissipated, and recovers as they do (Section 3.8.7 and Figure 3.8.4). The test builds a drawdown on a slope of the suite with three materials, whose heads fall from 45 m and 30 m on the two sides to 28 m and 26 m, with two flagged stages at 10⁴ and 10⁶ time units, and requires that the factor of safety of the second stage is higher than that of the first. It is a requirement of ordering and not a closed form, and it passes in version 0.1.235.

In OGR Slip2D

The sweep repeats the search at every level because the critical surface moves with the level: a surface held fixed would report a factor of safety that is too high, and a test guards that choice. The numbers of the sweep are the same through the interface and through the operations layer.

Source: tests/test_drawdown_sweep_v170.py · tests/test_groundwater_ops_v1200.py · tests/test_transient_v130.py

4.3.14References#

Darcy, H. (1856). Les fontaines publiques de la ville de Dijon. Victor Dalmont, Paris.
Gardner, W.R. (1958). Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Science, 85(4), 228–232.
Hoek, E., Carranza-Torres, C. & Corkum, B. (2002). Hoek–Brown failure criterion – 2002 edition. Proceedings of the 5th North American Rock Mechanics Symposium and 17th Tunnelling Association of Canada Conference (NARMS-TAC 2002), Toronto, Vol. 1, 267–273.
Morgenstern, N. (1963). Stability charts for earth slopes during rapid drawdown. Géotechnique, 13(2), 121–131.
van Genuchten, M.Th. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44(5), 892–898.
Carslaw, H.S. & Jaeger, J.C. (1959). Conduction of Heat in Solids, 2nd ed. Oxford University Press, Oxford.