MalcolmWoodruff/ TerraPercolatio

Finite-element groundwater seepage analysis and flow nets for 2D and axisymmetric problems

Contributed by: Malcolm Woodruff

TerraPercolatio solves steady-state groundwater seepage by the finite element method and draws true flow nets. It handles confined and unconfined problems in plane (2D) or axisymmetric geometry, with multi-layer, anisotropic soils. For unconfined problems it locates the phreatic surface and seepage faces iteratively, applies rainfall/surface infiltration through an unsaturated zone, and produces pore-pressure fields — including capped suction — ready for limit-equilibrium slope-stability analysis. Equipotentials come from the head field; flow lines come from the exact conjugate (stream-function) problem, so the flow net is correct in heterogeneous, anisotropic ground.

Installation Instructions

To install this paclet in your Wolfram Language environment, evaluate this code:
PacletInstall["MalcolmWoodruff/TerraPercolatio"]

To load the code after installation, evaluate this code:
Needs["MalcolmWoodruff`TerraPercolatio`"]

Details

Finite-element groundwater seepage analysis and flow nets for 2D and axisymmetric problems: multi-layer anisotropic soils, phreatic surfaces and seepage faces, rainfall infiltration, and pore-pressure output for slope stability.

See additional details and notes about the paclet in the User Manual.

Examples

Basic Examples (2)

Draw the flow net for an unconfined dam solution:

In[1]:=

Out[3]=

Plot the pore-pressure field:

In[4]:=

$InterpretationBox[FrameBox[TagBox[TooltipBox[PaneBox[GridBox[List[List[GraphicsBox[List[Thickness[0.0025`], List[FaceForm[List[RGBColor[0.9607843137254902`, 0.5058823529411764`, 0.19607843137254902`], Opacity[1.`]]], FilledCurveBox[List[List[List[0, 2, 0], List[0, 1, 0], List[0, 1, 0], List[0, 1, 0], List[0, 1, 0]], List[List[0, 2, 0], List[0, 1, 0], List[0, 1, 0], List[0, 1, 0], List[0, 1, 0]], List[List[0, 2, 0], List[0, 1, 0], List[0, 1, 0], List[0, 1, 0], List[0, 1, 0], List[0, 1, 0]], List[List[0, 2, 0], List[1, 3, 3], List[0, 1, 0], List[1, 3, 3], List[0, 1, 0], List[1, 3, 3], List[0, 1, 0], List[1, 3, 3], List[1, 3, 3], List[0, 1, 0], List[1, 3, 3], List[0, 1, 0], List[1, 3, 3]]], List[List[List[205.`, 22.863691329956055`], List[205.`, 212.31669425964355`], List[246.01799774169922`, 235.99870109558105`], List[369.0710144042969`, 307.0436840057373`], List[369.0710144042969`, 117.59068870544434`], List[205.`, 22.863691329956055`]], List[List[30.928985595703125`, 307.0436840057373`], List[153.98200225830078`, 235.99870109558105`], List[195.`, 212.31669425964355`], List[195.`, 22.863691329956055`], List[30.928985595703125`, 117.59068870544434`], List[30.928985595703125`, 307.0436840057373`]], List[List[200.`, 410.42970085144043`], List[364.0710144042969`, 315.7036876678467`], List[241.01799774169922`, 244.65868949890137`], List[200.`, 220.97669792175293`], List[158.98200225830078`, 244.65868949890137`], List[35.928985595703125`, 315.7036876678467`], List[200.`, 410.42970085144043`]], List[List[376.5710144042969`, 320.03370475769043`], List[202.5`, 420.53370475769043`], List[200.95300006866455`, 421.42667961120605`], List[199.04699993133545`, 421.42667961120605`], List[197.5`, 420.53370475769043`], List[23.428985595703125`, 320.03370475769043`], List[21.882003784179688`, 319.1406993865967`], List[20.928985595703125`, 317.4896984100342`], List[20.928985595703125`, 315.7036876678467`], List[20.928985595703125`, 114.70369529724121`], List[20.928985595703125`, 112.91769218444824`], List[21.882003784179688`, 111.26669120788574`], List[23.428985595703125`, 110.37369346618652`], List[197.5`, 9.87369155883789`], List[198.27300024032593`, 9.426692008972168`], List[199.13700008392334`, 9.203690528869629`], List[200.`, 9.203690528869629`], List[200.86299991607666`, 9.203690528869629`], List[201.72699999809265`, 9.426692008972168`], List[202.5`, 9.87369155883789`], List[376.5710144042969`, 110.37369346618652`], List[378.1179962158203`, 111.26669120788574`], List[379.0710144042969`, 112.91769218444824`], List[379.0710144042969`, 114.70369529724121`], List[379.0710144042969`, 315.7036876678467`], List[379.0710144042969`, 317.4896984100342`], List[378.1179962158203`, 319.1406993865967`], List[376.5710144042969`, 320.03370475769043`]]]]], List[FaceForm[List[RGBColor[0.5529411764705883`, 0.6745098039215687`, 0.8117647058823529`], Opacity[1.`]]], FilledCurveBox[List[List[List[0, 2, 0], List[0, 1, 0], List[0, 1, 0], List[0, 1, 0]]], List[List[List[44.92900085449219`, 282.59088134765625`], List[181.00001525878906`, 204.0298843383789`], List[181.00001525878906`, 46.90887451171875`], List[44.92900085449219`, 125.46986389160156`], List[44.92900085449219`, 282.59088134765625`]]]]], List[FaceForm[List[RGBColor[0.6627450980392157`, 0.803921568627451`, 0.5686274509803921`], Opacity[1.`]]], FilledCurveBox[List[List[List[0, 2, 0], List[0, 1, 0], List[0, 1, 0], List[0, 1, 0]]], List[List[List[355.0710144042969`, 282.59088134765625`], List[355.0710144042969`, 125.46986389160156`], List[219.`, 46.90887451171875`], List[219.`, 204.0298843383789`], List[355.0710144042969`, 282.59088134765625`]]]]], List[FaceForm[List[RGBColor[0.6901960784313725`, 0.5882352941176471`, 0.8117647058823529`], Opacity[1.`]]], FilledCurveBox[List[List[List[0, 2, 0], List[0, 1, 0], List[0, 1, 0], List[0, 1, 0]]], List[List[List[200.`, 394.0606994628906`], List[336.0710144042969`, 315.4997024536133`], List[200.`, 236.93968200683594`], List[63.928985595703125`, 315.4997024536133`], List[200.`, 394.0606994628906`]]]]]], List[Rule[BaselinePosition, Scaled[0.15`]], Rule[ImageSize, 10], Rule[ImageSize, 15]]], StyleBox[RowBox[List["SeepagePlot", " "]], Rule[ShowAutoStyles, False], Rule[ShowStringCharacters, False], Rule[FontSize, Times[0.9`, Inherited]], Rule[FontColor, GrayLevel[0.1`]]]]], Rule[GridBoxSpacings, List[Rule["Columns", List[List[0.25`]]]]]], Rule[Alignment, List[Left, Baseline]], Rule[BaselinePosition, Baseline], Rule[FrameMargins, List[List[3, 0], List[0, 0]]], Rule[BaseStyle, List[Rule[LineSpacing, List[0, 0]], Rule[LineBreakWithin, False]]]], RowBox[List["PacletSymbol", "[", RowBox[List["\"MalcolmWoodruff/TerraPercolatio\"", ",", "\"MalcolmWoodruff`TerraPercolatio`SeepagePlot\""]], "]"]], Rule[TooltipStyle, List[Rule[ShowAutoStyles, True], Rule[ShowStringCharacters, True]]]], Function[Annotation[Slot[1], Style[Defer[PacletSymbol["MalcolmWoodruff/TerraPercolatio", "MalcolmWoodruff`TerraPercolatio`SeepagePlot"]], Rule[ShowStringCharacters, True]], "Tooltip"]]], Rule[Background, RGBColor[0.968`, 0.976`, 0.984`]], Rule[BaselinePosition, Baseline], Rule[DefaultBaseStyle, List[]], Rule[FrameMargins, List[List[0, 0], List[1, 1]]], Rule[FrameStyle, RGBColor[0.831`, 0.847`, 0.85`]], Rule[RoundingRadius, 4]], PacletSymbol["MalcolmWoodruff/TerraPercolatio", "MalcolmWoodruff`TerraPercolatio`SeepagePlot"], Rule[Selectable, False], Rule[SelectWithContents, True], Rule[BoxID, "PacletSymbolBox"]][damSol, "PorePressure", "PlotPoints" -> 25]$

Out[4]=

Neat Examples (4)

A 10 m × 16 m excavation behind sheet-pile walls, dewatered to a floor 3 m below ground, in a 10 m stratum over an impermeable base. One quarter is modelled with symmetry planes at x = 0 and y = 0 (both natural no-flow boundaries). The walls and the excavated void are subtracted from the domain by constructive solid geometry, so wall faces are true no-flow boundaries. The far-field water table is held at ground level. All parameters are named below:

In[5]:=

$D3 = 10.; (* stratum thickness: impermeable base z = 0, ground z = 10 *) axq = 5.; (* excavation half-length in x (full length 10 m) *) byq = 8.; (* excavation half-width in y (full width 16 m) *) tw = 0.4; (* wall thickness, m *) ztip3 = 4.; (* wall toe elevation: 6 m walls, 3 m embedment below the floor *) zf = 7.; (* excavation floor = water level inside (dewatered), m *) hFar = 10.; (* far-field head: water table at ground level, m *) LxQ = 40.; LyQ = 40.; (* plan extent of the quarter model, m *) k3 = 1.*^-5; (* isotropic permeability, m/s *) wallX = Cuboid[{axq, 0., ztip3}, {axq + tw, byq + tw, D3}]; wallY = Cuboid[{0., byq, ztip3}, {axq + tw, byq + tw, D3}]; void = Cuboid[{0., 0., zf}, {axq, byq, D3}]; reg3 = RegionDifference[Cuboid[{0., 0., 0.}, {LxQ, LyQ, D3}], RegionUnion[wallX, wallY, void]]; mesh3 = ToElementMesh[reg3, MaxCellMeasure -> 30., MeshRefinementFunction -> Function[{v, vol}, Module[{c = Mean[v]}, vol > If[c[[1]] < 10 && c[[2]] < 13, 0.5, 30.]]], "MeshOrder" -> 2]; h3 = NDSolveValue[{Laplacian[h[x, y, z], {x, y, z}] == 0, DirichletCondition[h[x, y, z] == hFar, x >= LxQ - 10.^-4 || y >= LyQ - 10.^-4], DirichletCondition[h[x, y, z] == zf, z >= zf - 10.^-4 && z <= zf + 10.^-4 && x <= axq + 10.^-4 && y <= byq + 10.^-4]}, h, {x, y, z} \[Element] mesh3]$

Out[6]=

Total inflow (4 × the quarter-model flux through the floor) and exit gradients on the floor — note the 3D result that the pit corner is more critical than mid-wall:

In[7]:=

$dzh = Derivative[0, 0, 1][h3]; dxh = Derivative[1, 0, 0][h3]; Qtot3 = 4. NIntegrate[ k3 dzh[x, y, zf - 10.^-3], {x, 0., axq}, {y, 0., byq}, AccuracyGoal -> 4, PrecisionGoal -> 4]; {Abs[Qtot3], Abs[dzh[axq - 0.3, 1., zf - 10.^-3]], (* exit gradient, mid long wall *) Abs[dzh[axq - 0.3, byq - 0.3, zf - 10.^-3]], (* exit gradient, corner *) Abs[dzh[0.1, 0.1, zf - 10.^-3]] (* exit gradientcentre *)}$

Out[8]=

Section on the symmetry plane y = 0: equipotentials (red), streamlines (blue), water levels (▽) at ground level far away and at floor level in the pit:

In[9]:=

$soilY0 = Function[{x, z}, ! (axq <= x <= axq + tw && z >= ztip3) && ! (x <= axq && z >= zf)]; eqp = ContourPlot[h3[x, 0., z], {x, 0., 20.}, {z, 0., D3}, Sequence[ RegionFunction -> soilY0, Contours -> Range[7.25, 9.75, 0.25], ContourShading -> None, ContourStyle -> Directive[ RGBColor[0.75, 0.2, 0.15], AbsoluteThickness[1.2]], PlotPoints -> 70]]; stp = StreamPlot[{-dxh[x, 0., z], -dzh[x, 0., z]}, {x, 0., 20.}, {z, 0., D3}, Sequence[ RegionFunction -> soilY0, StreamPoints -> Fine, StreamColorFunction -> None, StreamStyle -> Directive[ RGBColor[0.1, 0.35, 0.7], AbsoluteThickness[1.]]]]; geo = Graphics[{ {EdgeForm[{GrayLevel[0.4]}], RGBColor[0.96, 0.93, 0.85], Rectangle[{0, 0}, {20, D3}]}, {GrayLevel[0.55], Rectangle[{0, zf}, {axq, D3}]}, {GrayLevel[0.25], Rectangle[{axq, ztip3}, {axq + tw, D3}]}, {RGBColor[0.1, 0.45, 0.85], AbsoluteThickness[1.6], Line[{{0, zf}, {axq, zf}}], Text[Style["\[EmptyDownTriangle]", RGBColor[0.1, 0.45, 0.85], 13], {2., zf}, {0, -0.9}], Line[{{axq + tw, D3}, {20, D3}}], Text[Style["\[EmptyDownTriangle]", RGBColor[0.1, 0.45, 0.85], 13], {14., D3}, {0, -0.9}]}}]; Show[geo, eqp, stp, Sequence[ Frame -> True, AspectRatio -> Automatic, ImageSize -> 800, PlotLabel -> "3D excavation - section on symmetry plane y = 0"]]$

Out[10]=

Plan view of the piezometric drawdown at mid-depth below the floor — rectangular near the pit, becoming circular with distance:

In[11]:=

$soilPlan = Function[{x, y}, ! (x <= axq + tw && y <= byq + tw)]; Show[ ContourPlot[h3[x, y, 5.5], {x, 0., 30.}, {y, 0., 30.}, Sequence[ RegionFunction -> soilPlan, Contours -> Range[7.5, 9.75, 0.25], ColorFunction -> "TemperatureMap", PlotLegends -> Automatic, PlotPoints -> 60]], Graphics[{{EdgeForm[{GrayLevel[0.25]}], GrayLevel[0.55], Rectangle[{0, 0}, {axq + tw, byq + tw}]}, Text[ Style["excavation\n(quarter)", White, 11], {2.7, 4.}]}], Sequence[ Frame -> True, AspectRatio -> Automatic, ImageSize -> 520, PlotLabel -> "plan: head at z = 5.5 m"]]$