Influence of a Fault on Steady-State Groundwater Flow

William C. Haneberg, bill@haneberg.com, www.haneberg.com

Introduction

This example uses analytical solutions from Haneberg ("Steady-State Groundwater Flow across Idealized Faults," Water Resources Research 31, 1815-1820, 1995) to calculate the influence of an idealized vertical fault on steady-state horizontal groundwater flow. The solutions are derived by specifying general solutions to the 1D steady-state flow equation for each of three subdomains (right-hand aquifer, fault, and left-hand aquifer). Constants of integration were eliminated by specifying head at the left and right boundaries, and then requiring continuity of head and flux across each of the two internal boundaries. Solution derivations are included at the end of this notebook. 

All variables are dimensionless. The horizontal axis shows distance from the fault (normalized to the problem domain width), and the vertical axis shows normalized hydraulic head changed (normalized to the total head drop across the problem domain).

Evaluate the FaultFlowPlot Function

Syntax is FaultFlowPlot[TL, TR, Lf], where TL is the transmissivity of the left-hand aquifer relative to that of the fault, TR is the transmissivity of the right-hand aquifer relative to that of the fault, and Lf is the width of the fault relative to the total width of the problem domain.

HL[TL_, TR_, 
        Lf_] := ((-TL + Lf TL + TR + Lf TR - 2 Lf TL TR + 4 TR x)/(2 (TL - Lf TL + TR - Lf TR + 2 Lf TL TR)))
HF[TL_, TR_, 
        Lf_] := ((-TL + Lf TL + TR - Lf TR + 4 TL TR x)/(2 (TL - Lf TL + TR - Lf TR + 2 Lf TL TR)))
HR[TL_, TR_, 
        Lf_] := ((-TL - Lf TL + TR - Lf TR + 2 Lf TL TR + 4 TL x)/(2 (TL - Lf TL + TR - Lf TR + 2 Lf TL TR)))
FaultFlowPlot[TL_, TR_, LF_] := 
      Block[{p1, p2, p3, x}, 
        p1 = Plot[HL[TL, TR, LF], {x, -0.5, -(LF/2)}, 
            DisplayFunction -> Identity]; 
        p2 = Plot[HF[TL, TR, LF], {x, -(LF/2), (LF/2)}, 
            DisplayFunction -> Identity]; 
        p3 = Plot[HR[TL, TR, LF], {x, (LF/2), 0.5}, 
            DisplayFunction -> Identity]; 
        Show[Graphics[{GrayLevel[0. ], Dashing[{0.02}], 
              Line[{{-0.5, -0.5}, {0.5, 0.5}}]}], p1, p2, p3, 
          AspectRatio -> 1, DefaultFont -> {"Helvetica", 12}, 
          Frame -> True, 
          FrameTicks -> {{-0.5, -0.25, 0, 0.25, 
                0.5}, {-0.5, -0.25, 0, 0.25, 0.5}, None, None}, 
          PlotRange -> {{-0.51, 0.51}, {-0.51, 0.51}}, 
          DisplayFunction -> $DisplayFunction]]

Create a Plot

In this case, the left-hand aquifer transmissivity is 10 times that of the fault and the right-hand aquifer transmissivity is 50 times that of the fault. The width of the fault relative to the problem domain is 0.05. The solid line shows the calculated hydraulic head profile across the fault, and the dashed line shows the unperturbed hydraulic gradient that would have existed without the fault. As discussed in Haneberg (1995), this very simple model accounts exceptionally well for water-level changes observed across faults in the Albuquerque basin.

FaultPlot = FaultFlowPlot[10, 50, 0.05]

Now, add some explanatory information and a background to highlight the location of the fault.

fault = 
        Graphics[{RGBColor[1, 1, 0.2], 
            Rectangle[{-0.025, -0.5}, {0.025, 0.5}]}];
leftaquifer = 
        Graphics[{RGBColor[0.5, 0.8, 1], 
            Rectangle[{-0.5, -0.5}, {-0.025, 0.5}]}];
rightaquifer = 
        Graphics[{RGBColor[0.5, 1, 1], 
            Rectangle[{0.025, -0.5}, {0.5, 0.5}]}];

Show[fault, leftaquifer, rightaquifer, FaultPlot, Frame -> True, 
      AspectRatio -> 1, DefaultFont -> {"Helvetica", 10}, 
      FrameTicks -> {{-0.5, -0.25, 0, 0.25, 
            0.5}, {-0.5, -0.25, 0, 0.25, 0.5}}, 
      FrameLabel -> {"Dimensionless Distance", "Dimensionless Head \
Change"}, 
      Epilog -> {Text["fault", {0, -0.3}], 
          Text["T = 1", {0, -0.4}], 
          Text["left aquifer", {-0.3, 0}], 
          Text["T = 10", {-0.3, -0.08}], 
          Text["right aquifer", {0.3, 0}], 
          Text["T = 50", {0.3, -0.08}]}]

Derivation of fault flow solutions

From Haneberg, W.C., "Steady-State Groundwater Flow across Idealized Faults,"
Water Resources Research 31 (1995) 1815-1820.

Define General Solutions

Define general solutions for heads in the left aquifer, fault, and right aquifer in which c1, ... ,c6 are the constants of integration. The governing equation for steady-state 1D flow in each of the three subdomains is (d^2)h/((dx)^2) = 0. DSolve can also be used to solve the entire system of differential equations and boundary conditions, and it is included as an alternative at the end of this notebook.

HL = c1 + c2 x

HF = c3 + c4 x

HR = c5 + c6 x

Evaluate the Boundary Conditions

At the left and right boundaries x  =  +/- 1/2. Head is fixed at +/- one half of the total head drop, Δ, across the system.

bc1 = 
      HL == (Δ/2) /. x -> -(1/2)

bc2 = 
      HR == -(Δ/2) /. x -> 
          (1/2)

Evaluate the internal boundary conditions that link the aquifers to the fault by requiring continuity of head and flux across the boundaries (x  =  +/-Lf/2). Flux is given in a vertically averaged form for 1D horizontal flow, q = T (dh/dx), where T is transmissivity (hydraulic conductivity times aquifer thickness).

bc3 = HL == HF /. x -> -(Lf/2)

bc4 = HF == HR /. x -> (Lf/2)

bc5 = 
      TL (∂_x) HL == 
          TF (∂_x) 
              HF /. x -> -(Lf/2)

bc6 = 
      TF (∂_x) HF == 
          TR (∂_x) HR /. x -> 
          (Lf/2)

Solve and Substitute

Solve the system of six boundary conditions in six variables (c1, ..., c6) to determine the constants of integration.

constants = 
      Solve[{bc1, bc2, bc3, bc4, bc5, bc6}, {c1, c2, c3, c4, c5, 
          c6}]

Finally, substitute the expressions for the constants back into the general solutions and simplify.

solHL  =  
      Simplify[HL /. constants[[1]]]

solHF  =  
      Simplify[HF /. constants[[1]]]

solHR  =  
      Simplify[HR /. constants[[1]]]

Using DSolve

This provides an alternative method to solve the system.

Clear[HL, HF, HR]

eq1=HL''[x]==0

eq2=HF''[x]==0

eq3=HR''[x]==0

eq4 = 
      HL[x] == (Δ/2) /. x -> -(1/2)

eq5 = 
      HR[x] == -(Δ/2) /. x -> 
          (1/2)

eq6 = 
      HL[x] == HF[x] /. x -> -(Lf/2)

eq7 = HF[x] == HR[x] /. x -> (Lf/2)

eq8=TL HL'[x]==TF HF'[x]/.x -> -(Lf/2)

eq9=TF HF'[x]==TR HR'[x]/.x -> (Lf/2)

DSolve[{eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9}, {HL[x], HF[x], 
        HR[x]}, x]

Simplify[%]

Compare

Using DSolve:

%

Solving for boundary conditions:

solHL

solHF

solHR