FLAC3D Theory and Background • Fluid-Mechanical Interaction

Semi-confined Aquifer


The project file for this example is available to be viewed/run in FLAC3D. The project’s main data files are shown at the end of this example.

Fluid leakage into a shallow semi-confined aquifer can be modeled with FLAC3D using the zone face apply leakage command. This is demonstrated for the example defined by the sketch in Figure 1. The aquifer has a length \(L\), height \(H\), and rests on an impermeable base. Fluid flow obeys Darcy’s law; the mobility coefficient \(k\) is homogeneous and isotropic. The semi-permeable top layer has permeability \(k_*\), and thickness \(H_*\). The effect of gravity is neglected in this example. Fluid pressure at the top of the leaky layer is constant and equal to \(p_*\). The lateral fluid-flow conditions correspond to a constant pressure \(p_0\) at the left boundary, and \(p_1\) at the right.

The objective is to determine the steady-state pore-pressure profile and total leakage into the aquifer. The general solution of pore pressure for a shallow semi-confined aquifer has the form (see Strack 1989)

(1)\[p - p_* = Ae^{x/\lambda} + Be^{-x/\lambda}\]

where \(\lambda\) is the seepage factor, which has the dimension of length and is defined as \(\lambda = \sqrt{k H H_*/k_*}\); \(A\) and \(B\) are constants determined from the pressure boundary conditions.


Figure 1: Shallow semi-confined aquifer.

The boundary conditions for this problem are

\(p = p_0\)


\(x = 0\)

\(p = p_1\)


\(x = L\)

The pore pressure solution is

(2)\[p = A e^{x/\lambda} + Be^{-x/\lambda} + p_*\]


(3)\[A = {{(p_1 - p_*)e^{L/\lambda} - (p_0 - p_*)} \over {e^{2L/\lambda} - 1}}\]
(4)\[B = {{(p_0 - p_*)e^{2L/\lambda} - (p_1 - p_*)e^{L/\lambda}} \over {e^{2L/\lambda} - 1}}\]

The steady state discharge over the height of the aquifer is obtained from Darcy’s law:

(5)\[Q_x = k H {dp \over dx}\]

After differentiation of Equation (2) with respect to \(x\), and substitution into Equation (5), we obtain

(6)\[Q = {{kH} \over \lambda}\Bigl\lbrack A e^{x/\lambda} - B e^{-x/\lambda} \Bigr\rbrack\]

The total amount of leakage into the aquifer is, by continuity of flow, equal to the difference between the discharge leaving at \(x\) = \(L\) and that entering at \(x\) = 0. Using Equation (6), we obtain, after some manipulation,

(7)\[Q_x = 2{{kH} \over \lambda} {{e^{L/\lambda} - 1} \over {e^{L/\lambda} + 1}}\Bigl\lbrack {p_0+p_1 \over 2} - p_* \Bigr\rbrack\]

Equation (6) and Equation (7) are used for comparison to the FLAC3D solution.

The FLAC3D data file for this problem is listed in “SemiConfinedAquifer.dat”. The analytical solution is programmed in FISH as part of the data file. The FLAC3D model is a 20 zone by 2 zone mesh with a constant pore pressure of \(p_0\) = 20 kPa applied at the left boundary, \(x\) = 0, and a constant pore pressure of \(p_1\) = 10 kPa applied at the right boundary, \(x\) = 20 m. A leaky aquifer boundary condition is applied along the top boundary of the model, \(y\) = 1 m, using the zone face apply leakage command. The pore pressure at the top is \(p_*\) = 1.8 kPa, and the leakage coefficient \(h\) (see this equation), is evaluated to be \(k_*/H_*\) = 2.98 × 10-9 m3 /(N sec), based on the properties of the leaky layer. The properties for this problem are listed in the setup function. The fluid-flow calculation mode is turned on, the mechanical calculation mode is turned off, and the simulation is run until steady-state flow is reached.

The FISH function checkit compares the amount of leakage calculated by FLAC3D to the solution of (7) at steady-state flow. The difference is printed (in a FISH dialog message) to be 0.03%. The analytical and numerical pore pressure profiles recorded along the base of the model, from \(x\) = 0 to \(x\) = 20, are compared in Figure 2.


Figure 2: Pore pressure profile.


Strack, O. D. L. Groundwater Mechanics. New Jersey: Prentice Hall (1989).


Data File


; file for 'Semiconfined aquifer'
model new
model large-strain off
model title 'Semiconfined aquifer'
fish automatic-create off
model configure fluid
; Geometry
zone create brick point 0 (0,0,0) point 1 (20,0,0) point 2 (0,1,0) ...
                  point 3 (0,0,1) size 20 2 2
zone face skin
; Constitutive model and properties
zone fluid cmodel assign isotropic
zone fluid property porosity=0.04 permeability=2e-8
zone gridpoint initialize fluid-tension -1e10
zone gridpoint initialize fluid-modulus 2e9
; --- boundary conditions ---
zone face apply leakage 1.8e4 3e-9 range group 'Top'
zone gridpoint fix pore-pressure 2e4 range group 'West'
zone gridpoint fix pore-pressure 1e4 range group 'East'
; --- fluid flow solution ---
model fluid active on
model mechanical active off
model solve fluid time-total 1
model save 'confaquifer'