Chapter 5: Flow Nets

Flow Nets

5.1 Flow Nets by Graphical Construction

We have seen in Chapter 2 that a groundwater flow system can be represented by a three-dimensional set of equipotential surfaces and a corresponding set of orthogonal flowlines. If a meaningful two-dimensional cross section can be chosen through the three-dimensional system, the set of equipotential lines and flowlines so exposed constitutes a flow net. The construction of flow nets is one of the most powerful analytical tools for the analysis of groundwater flow.

In Section 2.11 and Figure 2.25, we saw that a flow net can be viewed as the solution of a two-dimensional, steady-state, boundary-value problem. The solution requires knowledge of the region of flow, the boundary conditions on the boundaries of the region, and the spatial distribution of hydraulic conductivity within the region. In Appendix III, an analytical mathematical method of solution is presented. In this section, we will learn that flow nets can also be constructed graphically, without recourse to the sophisticated mathematics.

Homogeneous, Isotropic Systems

Let us first consider a region of flow that is homogenous, isotropic, and fully saturated. For steady-state flow in such a region, three types of boundaries can exist: (1) impermeable boundaries, (2) constant-head boundaries, and (3) water-table boundaries. First, let us consider flow in the vicinity of an impermeable boundary [Figure 5.1(a)]. Since there can be no flow across the boundary, the flowlines adjacent to the boundary must be parallel to it, and the equipotential lines must meet the boundary at right angles. By invoking Darcy’s law and setting the specific discharge across the boundary equal to zero, we are led to the mathematical statement of the boundary condition. For boundaries that parallel the axes in an xz plane:

\frac{\partial h}{\partial x} = 0 \hspace{1cm} \text{or} \hspace{1cm} \frac{\partial h}{\partial z} = 0(5.1)

Figure 5.1 Groundwater flow in the vicinity of (a) an impermeable boundary, (b) a constant-head boundary, and (c) a water-table boundary.

In effect, any flowline in a flow net constitutes an imaginary impermeable boundary, in that there is no flow across a flowline. In flow-net construction, it is often desirable to reduce the size of the region of flow by considering only those portions of the region on one side or the other of some line of symmetry. If it is clear that the line of symmetry is also a flowline, the boundary condition to be imposed on the symmetry boundary is that of Eq. (5.1).

A boundary on which the hydraulic head is constant [Figure 5.1(b)] is an equipotential line. Flowlines must meet the boundary at right angles, and adjacent equipotential lines must be parallel to the boundary. The mathematical condition is

h = c (5.2)

On the water table, the pressure head, ψ, equals zero, and the simple head relationship, h = ψ + z, yields

h = z (5.3)

for the boundary condition. As shown in Figure 5.1(c), for a recharge case the water table is neither a flowline nor an equipotential line. It is simply a line of variable but known h.

If we know the hydraulic conductivity K for the material in a homogenous, isotropic region of flow, it is possible to calculate the discharge through the system from a flow net. Figure 5.2 is a completed flow net for the simple case first presented in Figure 2.25(a). The area between two adjacent flowlines is known as a streamtube or flowtube. If the flowlines are equally spaced, the discharge through each streamtube is the same. Consider the flow through the region ABCD in Figure 5.2. If the distances AB and BC are ds and dm, respectively, and if the hydraulic-head drop between AD and BC is dh, the discharge across this region through a cross-sectional area of unit depth perpendicular to the page is

dQ = K\frac{dh}{ds}dm (5.4)

Figure 5.2 Quantitative flow net for a very simple flow system.

Under steady-state conditions, the discharge across any plane of unit depth (say, at AD, EH, or FG) within the streamtube must also be dQ. In other words, the discharge through any part of a streamtube can be calculated from a consideration of the flow in just one element of it.

If we arbitrarily decide to construct the flow net in squares, with ds = dm, then Eq. (5.4) becomes

dQ = K dh (5.5)

For a system with m streamtubes, the total discharge is

Q = mK dh (5.6)

If the total head drop across the region of flow is H and there are n divisions of head in the flow net (H = n dh), then

Q = \frac{mKH}{n} (5.7)

For Figure 5.2, m = 3, n = 6, H = 60 m, and from Eq. (5.7), Q = 30 K. For K = 10-4 m/s, Q = 3 × 10-3 m3/s (per meter of section perpendicular to the flow net).

Equation (5.7) must be used with care. It is applicable only to simple flow systems with one recharge boundary and one discharge boundary. For more complicated systems, it is best to simply calculate dQ for one streamtube and multiply by the number of streamtubes to get Q.

Figure 5.3 is a flow net that displays the seepage beneath a dam through a foundation rock bounded at depth by an impermeable boundary. It can be used to make three additional points about flow-net construction.

Figure 5.3 Seepage beneath a dam through homogenous, isotropic foundation rocks.
  1. The “squares” in all but the simplest flow nets are actually “curvilinear” squares; that is, they have equal central dimensions; or viewed another way, they enclose a circle that is tangent to all four bounding lines.
  2. It is not necessary that flow nets have finite boundaries on all sides; regions of flow that extend to infinity in one or more directions, like the horizontally infinite layer in Figure 5.3, are tractable.
  3. A flow net can be constructed with a “partial” streamtube on the edge.

For the flow net shown in Figure 5.3, m = 3\frac{1}{2}. If H = 100 m and K = 10-4 m/s, then, since n = 6, we have Q = 5.8 × 10-3 m3/s (per meter section perpendicular to the flow net).

In homogeneous, isotropic media, the distribution of hydraulic head depends only on the configuration of the boundary conditions. The qualitative nature of the flow net is independent of the hydraulic conductivity of the media. The hydraulic conductivity comes into play only when quantitative discharge calculations are made. It is also worth noting that flow nets are dimensionless. The flow nets of Figures 5.2 and 5.3 are equally valid whether the regions of flow are considered to be a few meters square or thousands of meters square.

The sketching of flow nets is something of an art. One usually pursues the task on a trial-and-error basis. Some hydrologists become extremely talented at arriving at acceptable flow nets quickly. For others, it is a source of continuing frustration. For a flow net in homogeneous, isotropic media, the rules of graphical construction are deceptively simple. We can summarize them as follows: (1) flowlines and equipotential lines must intersect at right angles throughout the system; (2) equipotential lines must meet impermeable boundaries at right angles; (3) equipotential lines must parallel constant-head boundaries; and (4) if the flow net is drawn such that squares are created in one portion of the field, then, with the possible exception of partial flow tubes at the edge, squares must exist throughout the entire field.

Heterogeneous Systems and the Tangent Law

When groundwater flowlines cross a geologic boundary between two formations with different values of hydraulic conductivity, they refract, much as light does when it passes from one medium to another. However, in contradistinction to Snell’s law, which is a sine law, groundwater refraction obeys a tangent law.

Consider the streamtube shown in Figure 5.4. Flow proceeds from a medium with hydraulic conductivity K1 to a medium with hydraulic conductivity K2, where K2 > K1.

Figure 5.4 Refraction of flowlines at a geologic boundary.

The streamtube has a unit depth perpendicular to the page, and the angles and distances are as indicated on the figure. For steady flow, the inflow Q1 must equal the outflow Q2; or, from Darcy’s law,

K_1a\frac{dh_1}{dl_1} = K_2c\frac{dh_2}{dl_2} (5.8)

where dh1 is the head drop across the distance dl1, and dh2, is the head drop across the distance dl2. In that dl1 and dl2 bound the same two equipotential lines, it is clear that dh1 = dh2; and from geometrical considerations, a = b cos θ1 and c = b cos θ2. Noting that b/dl1 = 1/sinθ1, and b/dl2 = 1/sinθ2, Eq. (5.8) becomes

K_1\frac{\cos \theta_1}{\sin \theta_1} = K_2\frac{\cos \theta_2}{\sin \theta_2} (5.9)


\frac{K_1}{K_2} = \frac{\tan \theta_1}{\tan \theta_2} (5.10)

Equation (5.10) constitutes the tangent law for the refraction of groundwater flowlines at a geologic boundary in heterogeneous media. Knowing K1, K2, and θ1, one can solve Eq. (5.10) for θ2. Figure 5.5 shows the flowline refractions for two cases with K1/K2 = 10. Flowlines, as if they had a mind of their own, prefer to use high-permeability formations as conduits, and they try to traverse low-permeability formations by the shortest route. In aquifer-aquitard systems with permeability contrasts of 2 orders of magnitude or more, flowlines tend to become almost horizontal in the aquifers and almost vertical in the aquitards. When one considers the wide range in hydraulic conductivity values exhibited in Table 2.2, it is clear that contrasts of 2 orders of magnitude and more are not at all uncommon.

Figure 5.5 Refraction of flowlines in layered systems (after Hubbert, 1940).

If one attempts to draw the equipotential lines to complete the flow systems on the diagrams of Figure 5.5, it will soon become clear that it is not possible to construct squares in all formations. In heterogeneous systems, squares in one formation become rectangles in another.

We can summarize the rules for graphical flow net construction in heterogeneous, isotropic systems as follows: (1) flowlines and equipotential lines must intersect at right angles throughout the system; (2) equipotential lines must meet impermeable boundaries at right angles; (3) equipotential lines must parallel constant-head boundaries; (4) the tangent law must be satisfied at geologic boundaries; and (5) if the flow net is drawn such that squares are created in one portion of one formation, squares must exist throughout that formation and throughout all formations with the same hydraulic conductivity. Rectangles will be created in formations of different conductivity.

The last two rules make it extremely difficult to draw accurate quantitative flow nets in complex heterogeneous systems. However, qualitative flow nets, in which the orthogonality is preserved but no attempt is made to create squares, can be of great help in understanding a groundwater flow system. Figure 5.6 is a qualitatively sketched flow net for the dam seepage problem first introduced in Figure 5.3, but with a foundation rock that is now layered.

Figure 5.6 Seepage beneath a dam through heterogeneous, isotopic foundation rocks.

Anisotropic Systems and the Transformed Section

In homogeneous but anisotropic media, flow-net construction is complicated by the fact that flowlines and equipotential lines are not orthogonal. Maasland (1957), Bear and Dagan (1965), and Liakopoulos (l965b) provide discussions of the theoretical principles that underlie this phenomenon, and Bear (1972) presents an extensive theoretical review. In this section, we shall look primarily at the practical response that has been devised to circumvent the conditions of nonorthogonality. It involves flow-net construction in the transformed section.

Consider the flow in a two-dimension region in a homogeneous, anisotropic medium with principal hydraulic conductivities Kx and Kz. The hydraulic-conductivity ellipse (Figure 5.7) will have semiaxes and \sqrt{K_x} and \sqrt{K_z}.

Figure 5.7 Hydraulic conductivity ellipse for an anisotropic medium with Kx/Kz = 5. The circles represent two possible isotropic transformations.

Let us transform the scale of the region of flow such that coordinates in the transformed region with coordinate directions X and Z will be related to those in the original xy system by

X = x

Z = \frac{z\sqrt{K_x}}{\sqrt{K_z}}

For Kx > Kz, this transformation will expand the vertical scale of the region of flow. It will also expand the hydraulic conductivity ellipse into a circle of radius \sqrt{K_x} (the outer circle in Figure 5.7); and the fictitious, expanded region of flow will then act as if it were homogenous with conductivity Kx.

The validity of this transformation can be defended on the basis of the steady state equation of flow. In the original xz coordinate system, for an anisotropic media, we have, from Eq. (2.69),

\frac{\partial}{\partial x}\left(K_x\frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial z}\left(K_z\frac{\partial h}{\partial x}\right) = 0

Dividing through by Kx, yields

\frac{\partial^2 h}{\partial x^2} + \frac{\partial}{\partial z} \left(\frac{K_z}{K_x}\frac{\partial h}{\partial x}\right) = 0

For the transformed section, we have from the second expression of Eq. (5.11),

\frac{\partial}{\partial z} + \frac{\sqrt{K_x}}{\sqrt{K_z}}\frac{\partial}{\partial Z}

Noting the first expression of Eq. (5.13), yields

\frac{\partial^2 h}{\partial X^2} + \frac{\partial^2 h}{\partial Z^2} = 0 (5.15)

which is the equation of flow for a homogenous, isotropic medium in the transformed section.

An equally valid transformation could be accomplished by contracting the region in the x direction according to the relations

X = \frac{x\sqrt{K_z}}{\sqrt{K_x}}

Z = z

In this case, the conductivity ellipse will be transformed into the small circle in Figure 5.7, and the fictitious, transformed medium will act as if it were homogenous with hydraulic conductivity Kz.

With the concept of transformed section in hand, the steps in the graphical construction of a flow net in a homogeneous, anisotropic medium become self evident: (1) carry out a transformation of coordinates using either Eqs. (5.11) or Eqs. (5.16); (2) construct a flow net in the fictitious, transformed section, according to the rules for a homogenous, isotropic media; and (3) invert the scaling ratio.

Figure 5.8 (a) Flow problem in a homogenous anisotropic region with \sqrt{K_x}/\sqrt{K_z} = 4. (b) Flow net in the transformed isotropic section. (c) Flow net in the actual anisotropic section. T, transformation; I, inversion.

Figure 5.8 is an example of the technique. The boundary-value problem illustrated in Figure 5.8(a) is a vertical section that represents flow from a surface pond at h = 100 toward a drain at h = 0. The drain is considered to be one of many parallel drains set at a similar depth oriented perpendicular to the page. The vertical impermeable boundaries are “imaginary”; they are created by the symmetry of the overall flow system. The lower boundary is a real boundary; it represents the base of the surficial soil, which is underlain by a soil or rock formation with a conductivity several orders of magnitude lower. If the vertical axis is arbitrarily set with z = 100 at the drain and z = 100 at the surface, then from h = ψ + z, and the h values given, we have ψ = 0 at both boundaries. At the surface, this condition implies that the soil is just saturated. The “pond” is incipient; it has zero depth. At the drain, ψ = 0 implies free-flowing conditions. The soil in the flow field has an anisotropic conductivity of Kx/Kz = 16. The transformed section of Figure 5.8(b) therefore has a vertical expansion of \sqrt{K_x}/\sqrt{K_z} = 4. Figure 5.8(c) shows the result of the inverse transformation, wherein the homogeneous, isotropic flow net from the transformed section is brought back into the true-scale region of flow. Under the inversion, the hydraulic head at any point (X, Z) in Figure 5.8(b) becomes the hydraulic head at point (x, z) in Figure 5.8(c).

The size of the transformed section is obviously dependent on whether Eqs (5.11) or Eqs. (5.16) are used for the transformation, but the shape of the region and the resulting flow net are the same in either case.

If discharge quantities or flow velocities are required, it is often easiest to make these calculations in the transformed section. The question then arises as to what hydraulic conductivity value ought to be used for such calculations. It is clear that it would be incorrect to use Kx, for a vertically expanded section and Kz, for a horizontally contracted one, as might be inferred from Figure 5.7, for this would produce two different sets of quantitative calculations for he two equivalent representations of the same problem. In fact, the correct value to use is

K' = \sqrt{K_x \cdot K_z} (5.17)

The validity of Eq. (5.17) rests on the condition that flows in each of the two equivalent transformed representations of the flow region must be equal. The proof requires an application of Darcy’s law to a single flowtube in each of the two transformations.

The influence of anisotropy on the nature of groundwater flow nets is illustrated in Figure 5.9 for the same boundary-value problem that was brought into play in Figure 5.8 The most important feature of the anisotropic flow nets [Figure 5.9(a) and 5.9(c)] is their lack of orthogonality. It seems to us that the transformation techniques introduced in this section provide an indirect but satisfying explanation of this phenomenon.

Figure 5.9 Flow nets for the flow problem of Figure 5.8(a) for \sqrt{K_x}/\sqrt{K_z} = (a) ¼, (b) 1, (c) 4 (after Maasland, 1957).

There are many situations where one may wish to construct a flow net on the basis of piezometric data from the field. If the geologic formations are known to be anisotropic, great care must be exercised in the inference of flow directions from the equipotential data. If the complex flow net is desired, a transformed section is in order, but if flow directions at specific points are all that is required, there is a graphical construction that can be useful. In Figure 5.10 the dashed line represents the directional trend of an equipotential line at some point of interest within an xz field. An inverse hydraulic-conductivity ellipse is then constructed about the point. This ellipse has principal semiaxes 1/\sqrt{K_x} and 1/\sqrt{K_z} (rather than \sqrt{K_x} and \sqrt{K_z}, as in Figure 5.7). A line drawn in the direction of the hydraulic gradient intersects the ellipse at the point A. If a tangent is drawn to the ellipse at A, the direction of flow is perpendicular to this tangent line. As an example of the application of this construction, one might compare the results of Figure 5.10 with the flowline/equipotential line intersections in the right-central portion of Figure 5.9(c).

Figure 5.10 Determination of direction of flow in an anisotropic region with \sqrt{K_x}/\sqrt{K_z} = 5.

5.2 Flow Nets by Analog Simulation

For flow in homogeneous, isotropic medium in an xz coordinate system, the equipotential lines in a flow net are a contoured reflection of the solution, h(x, z), of the boundary-value problem that describes steady-flow state in the region. The construction of the flow net is an indirect solution to Laplace’s equation:

\frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial z^2} = 0 (5.18)

This equation is one of the most commonly occurring partial differential equations in mathematical physics. Among he other physical phenomena hat it describes are the flow of heat through solids and the flow of electrical current through conductive media. For the latter case, Laplace’s equation takes the form

\frac{\partial^2 V}{\partial x^2} + \frac{\partial^2 V}{\partial z^2} = 0 (5.19)

where V is the electrical potential or voltage.

The similarity of Eqs. (5.18) and (5.19) reveals a mathematical and physical analogy between electrical flow and groundwater flow. The two equations are both developed on the basis of a linear flow law, Darcy’s law in one case and Ohm’s law in the other; and a continuity relationship, the conservation of fluid mass in one case and the conservation of electrical charge in the other. A comparison of Ohm’s law,

l_x = -\sigma\frac{\partial V}{\partial x} (5.20)

and Darcy’s law,

v_x = -K\frac{\partial h}{\partial x} (5.21)

clarifies the analogy at once. The specific discharge, vx (discharge per unit area), is analogous to the current density, Ix (electrical current per unit area); the hydraulic conductivity, K, is analogous to the specific electrical conductivity, σ; and the hydraulic head, h, is analogous to the electrical potential, V.

The analogy between electrical flow and groundwater flow is the basis for two types of analog model that have proven useful for the generation of quantitative flow nets. The first type involves the use of conductive paper and the second type utilizes resistance networks.

Conductive-Paper Analogs

Let us consider once again the hydraulic problem first shown in Figure 5.8 and now reproduced in Figure 5.11(a). The electric analog [Figure 5.11 (b)] consists of a sheet of conductive paper cut in the same geometrical shape as the groundwater flow field. A power supply is used to set up a voltage differential across the boundaries, and a sensing probe connected to the circuit through a voltmeter is used to measure the potential distribution throughout the conductive sheet. Constant-head boundaries, such as the V = 100 boundary on Figure 5.11(b), are created with highly conductive silver paint; impermeable boundaries are simulated by the unconnected edges of the paper model. It is usually possible to search out the equipotential lines rather efficiently, so that a full equipotential net can be quickly generated.

The method is limited to homogeneous, isotropic systems in two dimensions, but it is capable of handling complex region shapes and boundary conditions. Variations in the conductivity of commercially available paper may lead to random errors that limit the quantitative accuracy of the method. Two of the most detailed applications of the method are Childs’ (1943) theoretical analysis of near surface flow systems in drained land, and Tóth’s (1968) consideration of regional groundwater flow nets for a field area in Alberta.

Figure 5.11 Flow nets by electric-analog simulation. (a) Steady-State hydrogeologic boundary-value problem; (b) conducting-paper analog; (c) resistance-network analog.

Resistance Network Analogs

The use of a resistance network as an electrical analog is based on the same principles as the conductive-paper analog. In this approach [Figure 5.11(c)] the flow field is replaced by a network of resistors connected to one another at the nodal points of a grid. The flow of electricity through each resistor is analogous to the flow of groundwater through a flow tube parallel to the resistor and having a cross-sectional area reflected by the resistor spacing times a unit depth. For electrical flow through an individual resistor, the I in Eq. (5.20) must now be viewed as the current, and the σ is equal to 1/R, where R is the resistance of the resistor. As in the paper analog, a potential difference is set up across the constant-head boundaries of the model. A sensing probe is used to determine the voltage at each of the nodal points in the network, and these values, when recorded and contoured create the equipotential net.

By varying the resistances in the network it is possible to analyze heterogeneous and anisotropic systems with resistance-network analogs. They have an accuracy and versatility that is superior to the paper models, but they are not as flexible as the numerical methods introduced in the next section.

Karplus (1958) provides a detailed handbook for analog simulation Resistance-network analogs have been used to generate groundwater flow nets by Luthin (1953) in a drainage application, and by Bouwer and Little (1959) for a saturated-unsaturated system. Bouwer (1962) utilized the approach to analyze the configuration of groundwater mounds that develop beneath recharge ponds.

The most widespread use of electrical analog methods in groundwater hydrology is in the form of resistance-capacitance networks for the analysis of transient flow in aquifers. This application will be discussed in Section 8.9.

5.3 Flow Nets by Numerical Simulation

The hydraulic-head field, h(x, z), that allows construction of a flow net can be generated mathematically from the pertinent steady-state boundary-value problem in two ways. The first approach utilizes analytical solutions as discussed in Section 2.11 and Appendix III; the second approach uses numerical methods of solution. Analytical methods are limited to flow problems in which the region of flow, boundary conditions, and geologic configuration are simple and regular. As we shall see in this section, numerical methods are much more versatile, but their application is bound unequivocally to the use of a digital computer.

Numerical methods are approximate. They are based on a discretization of the continuum that makes up the region of flow. In the discretization, the region is divided into a finite number of blocks, each with its own hydrogeologic properties, and each having a node at the center at which the hydraulic head is defined for the entire block. Figure 5.12(a) shows a 7 × 5 block-centered nodal grid (i = 1 to i = 7 in the x direction, and j = 1 to j = 5 in the z direction) for a rectangular flow region.

Let us now examine the flow regime in the vicinity of one of the interior nodes—say, in the nodal block, i = 4, j = 3, and its four surrounding neighbors. To simplify the notation, we will renumber the nodes as indicated in Figure 5.12(b). If flow occurs from node 1 to node 5, we can calculate the discharge, Q15, from Darcy’s law:

Q_{15} = K_{15}\frac{h_1 - h_5}{\Delta z}\Delta x (5.22)

for flow through a cross section of unit depth perpendicular to the page. On the assumption that flow is directed toward the central node in each case, we can write down similar expressions for Q25, Q35, and Q45. For steady-state flow, consideration of the conservation of fluid mass decrees that the sum of these four flows must be zero. If the medium is homogeneous and isotropic, K15 = K25 = K35 = K45, and if we arbitrarily select a nodal grid that is square so that Δx = Δz, summation of the four terms leads to

h_5 = \frac{1}{4}(h_1 + h_2 + h_3 + h_4) (5.23)

This equation is known as a finite-deference equation. If we revert to the ij notation of Figure 5.12(a), it becomes

h_{i,j} = \frac{1}{4}(h_{i,j-1} + h_{i+1,j} + h_{i,j+1} + h_{i-1,j}) (5.24)

Figure 5.12 (a) Block-centered nodal grid for numerical flow-net simulation. (b) Interior node and its neighboring nodal blocks. (c) Finite-difference stars for an interior node and for nodes on a basal impermeable boundary and an impermeable corner.

In this form, Eq. (5.24) holds for all internal modes in the nodal grid. It states an elegantly simple truth: in steady flow, in a homogeneous, isotropic medium, the hydraulic head at any node is the average of the four surrounding values.

A similar exercise will reveal that the finite-difference equation for a node along the basal boundary, assuming that boundary to be impermeable, takes the form

h_{i,j} = \frac{1}{4}(h_{i-1,j} + h_{i+1,j} + 2h_{i,j+1}) (5.25)

and, at a corner node,

h_{i,j} = \frac{1}{4}(2h_{i-1,j} + 2h_{i,j+1}) (5.26)

Equations (5.24) through (5.26) are schematically represented, in a self-explanatory way, by the three finite-difference stars of Figure 5.12(c).

In short, it is possible to develop a finite-difference equation for every node in the nodal grid. If there are N nodes, there are N finite-difference equations. There are also N unknowns—the N values of h at the N nodes. We have therefore produced N linear, algebraic equations in N unknowns. If N were very small, we could solve the equations directly, using a technique such as Cramer’s rule, but where N is large, as it invariably is in numerical-flow-net simulation, we must introduce a more efficient method, known as relaxation.

Let us remain faithful to the flow problem of Figure 5.11(a) and assume that we wish to develop its flow net by numerical means. In the nodal grid of Figure 5.12(a), the nodal points at which the head is known are circled: h = 0 at i = 1, j = 3, and h = 100 at all the nodes in the j = 5 row. Relaxation involves repeated sweeps through the nodal net, top to bottom and left to right (or in any consistent fashion), applying the pertinent finite-difference equation at each node where the head is unknown. One must assume some starting h value at each node. For the problem at hand, h = 50 could be assigned as the starting value at all the uncircled nodes. In the application of the finite-difference equations, the most recently calculated h value is always used at every node. Each pass through the system is called an iteration, and after each iteration the calculated h values will approach more closely their final answers. The difference in h at any node between two successive iterations is called the residual. The maximum residual in the system will decrease as the iterations proceed. A solution has been reached when the maximum residual is reduced below a predetermined tolerance.

To test one’s understanding of the relaxation process, the reader might carry out a couple of iterations in the upper left-hand portion of the net. If the initial assigned value at the node, i = 2, j = 4, for example, is 50, then the value after the first iteration is 62.5 and after the second is 64.0. The final value, attained after several iterations, lies between 65 and 66.

Numerical simulation is capable of handling any shape of flow region and any distribution of boundary conditions. It is easy to redevelop the finite-difference equations for rectangular meshes where Δx ≠ Δz, and for heterogeneous, anisotropic conductivity distributions where the Kx and Kz values vary from node to node. In Eq. (5.22), the usual averaging technique would set K15 = (K1 + K5)/2, where K1, and K5, in this case refer to the vertical conductivities at nodes 1 and 5, and these might differ from each other and from the horizontal conductivities at these nodes. Numerical simulation allows flow-net construction in cases that are too complex for graphical construction or analytical solution. Numerical simulation is almost always programmed for the digital computer, and computer programs are usually written in a generalized form so that only new data cards are required to handle vastly differing flow problems. This is a distinct advantage over resistance-network analogs, which require a complete breakdown of the assembled hardware to effect a new simulation.

The development of the finite-equations presented in this section was rather informal. It is possible to begin with Laplace’s equation and proceed more mathematically toward the same result. In Appendix VI, we present a brief development along these lines. Perhaps it is worth noting in passing that the informal development utilizes Darcy’s law and the continuity relationship to reach the finite-difference expressions. These are the same two steps that led to the development of Laplace’s equation in Section 2.11.

The method we have called relaxation (after Shaw and Southwell, 1941) has several aliases. It is variously known as the Gauss-Seidel method, the Liebmann method, and the method of successive displacements. It is the simplest, but far from the most efficient, of many available methods for solving the set of finite-difference equations. For example, if the computed heads during relaxation are corrected according to

h_{corr}^k = \omega h^k + (1 - \omega)h_{corr}^{k-1} (5.27)

where hk is the computed head at the kth iteration and h_{corr}^{k-1} is the corrected head from the previous iteration, then the method is known as successive overrelaxation and the number of iterations required to reach a converged solution is significantly reduced. The parameter ω is known as the overrelaxation parameter, and it must lie in the range 1 ≤ ω ≤ 2.

There are many texts that will serve the numerical modeler well. McCracken and Dorn (1964) provide an elementary introduction to computer simulation techniques in their Fortran manual. Forsythe and Wasow (1960) deliver their message at a more advanced mathematical level. Remson et al. (1971) discuss a broad spectrum of numerical techniques with particular reference to groundwater flow. Pinder and Gray (1977) treat the subject at a more advanced level.

Numerical methods were introduced to the groundwater hydrology literature by Stallman (1956) in an analysis of regional water levels. Fayers and Sheldon (1962) were among the first to advocate steady-state numerical simulation in the study of regional hydrogeology. Remson et al. (1965) used the method to predict the effect of a proposed reservoir on groundwater levels in a sandstone aquifer. Freeze and Witherspoon (1966) generated many numerical flow nets in their theoretical study of regional groundwater flow. The method was in wide use much earlier in the agricultural drainage field (see Luthin and Gaskell, 1950) and in the derivation of seepage patterns in earth dams (Shaw and Southwell, 1941).

In recent years, the finite-difference method has been equaled in popularity by another numerical method of solution, known as the finite-element method. It, too, leads to a set of N equations in N unknowns that can be solved by relaxation, but the nodes in the finite-element method are the corner points of an irregular triangular or quadrilateral mesh that is designed by the modeler for each specific application, rather than the regular rectangular mesh of the finite-difference method. In many cases, a smaller nodal grid suffices and there are resulting economies in computer effort. The finite-element method is also capable of handling one situation that the finite-difference method cannot. The finite-difference method requires that the principal directions of anisotropy in an anisotropic formation parallel the coordinate directions. If there are two anisotropic formations in a flow field, each with different principal directions, the finite-difference method is stymied, whereas the finite-element method can provide a solution. The development of the finite-element equations requires a mathematical sophistication that is out of place in this introductory text. The interested reader is referred to Pinder and Gray (1977). Numerical methods, both finite-difference and finite-element, are widely used as the basis for digital computer simulation of transient flow in groundwater aquifers. This application is discussed in Section 8.8.

5.4 Saturated-Unsaturated Flow Nets

There is another type of flow net that is extremely difficult to construct by graphical means. For flow problems that involve both saturated and unsaturated flow, steady-state flow nets are usually derived by numerical simulation. Consider the flow net illustrated in Figure 5.13. It is similar to the problem that we have repeatedly analyzed in the past sections in that it involves flow to a drain in a system with impermeable boundaries on three sides, but it differs in that the vertical scale has been set up in such a way that the hydraulic head on the upper boundary now infers a pressure-head value that is less than atmospheric. This means that the soil is unsaturated at the surface, although if outflow to the drain is to take place, it must be saturated at depth. The qualitative flow net in Figure 5.13 has been developed for a soil whose unsaturated characteristic curves are those shown on the inset graphs.

Figure 5.13 Saturated-unsaturated flow net in a homogeneous, isotropic soil. The insets show the unsaturated characteristic curves for the soil.

These curves of hydraulic conductivity, K, and moisture content, θ as a function of ψ, are the wetting curves taken from Figure 2.13.

As in the one-dimensional saturated-unsaturated case that was schematically illustrated in Figure 2.12, there are three types of output from a numerically simulated, steady-state, saturated-unsaturated flow net in two dimensions. First there is the hydraulic-head pattern, h(x, z), that allows construction of the equipotential net (the dashed lines on Figure 5.13). Second, there is the pressure-head pattern, ψ(x, z) (the dotted lines in Figure 5.13), which is of particular value in defining the position of the water table (the ψ = 0 isobar). Third, there is the moisture content pattern, θ(x, z), which can be determined from the ψ(x, z) pattern with the aid of the θ(ψ) curve for the soil. For example, along the ψ = -50 cm dotted line in Figure 5.13, the moisture content, θ, is 27%.

The flowlines and equipotential lines form a continuous net over the full saturated-unsaturated region. They intersect at right angles throughout the system. A quantitative flow net could be drawn with curvilinear squares in the homogeneous, isotropic, saturated portion, but such flow tubes would not exhibit squares as they traverse the unsaturated zone, even in homogeneous, isotropic soil. As the pressure head (and moisture content) decrease, so does the hydraulic conductivity, and increased hydraulic gradients are required to deliver the same discharge through a given flow tube. This phenomenon can be observed in the flow tubes in the upper left-hand corner of the flow net in Figure 5.13, where the gradients increase toward the surface.

The concept of an integrated saturated-unsaturated flow system was introduced to the hydrologic literature by Luthin and Day (1955). They utilized numerical simulation and an experimental sand tank to derive their h(x, z) pattern. Bouwer and Little (1959) used an electrical resistance network to analyze tile drainage and subirrigation problems similar in concept to that shown in Figure 5.13. Saturated-unsaturated flow nets are required to explain perched water tables (Figures 2.15 and Figure 6.11), and to understand the hydrogeological regime on a hillslope as it pertains to streamflow generation (Section 6.5). Reisenauer (1963) and Jeppson and Nelson (1970) utilized numerical simulation to look at the saturated-unsaturated regime beneath ponds and canals. Their solutions have application to the analysis of artificial recharge of groundwater (Section 8.11). Freeze (1971b) considered the influence of the unsaturated zone on seepage through earth dams (Section 10.2).

5.5 The Seepage Face and Dupuit Flow

Seepage Face, Exit Point, and Free Surface

If a saturated-unsaturated flow system exists in the vicinity of a free-outflow boundary, such as a streambank or the downstream face of an earth dam, a seepage face will develop on the outflow boundary. In Figure 5.14(a), BC is a constant-head boundary and DC is impermeable. If there is no source of water at the surface, AB will also act like an impermeable boundary.

Figure 5.14 Development of a seepage face on a free-outflow boundary. (a) Saturated-unsaturated flow net; (b) free-surface flow net; (c) Dupuit-Forchheimer flow net.

The water table EF intersects the outflow boundary AD at the exit point E. All the flow must leave the system across the seepage face ED below the exit point E. Above E, along the line AE, the unsaturated pressure heads, ψ, are less than atmospheric, so outflow to the atmosphere is impossible. In effect, AE acts as an impermeable boundary. The condition on ED is h = z, the same as that on the water table. The problem in preparing a flow net for such cases lies in the fact that the position of the exit point that separates the two boundary conditions on the outflow boundary is not known a priori. In numerical simulation, it is necessary to provide an initial prediction for the position of the exit point. The correct exit point is then determined by a series of trial-and-error steady-state solutions.

The construction of a quantitative flow net in a saturated-unsaturated regime requires knowledge of both the saturated hydraulic conductivity, K, and the unsaturated characteristic curve, K(ψ), for the soil. In many engineering applications, including the analysis of seepage through earth dams, the latter data are seldom available. In such cases, the assumption is usually made that flow through the unsaturated portion of the system is negligible, or viewed another way, that the hydraulic conductivity of the soil at moisture contents less than saturation is negligible compared to the saturated hydraulic conductivity. In this case, the upper boundary of the flow net becomes the water table, and the water table itself becomes a flowline. Under these special circumstances, this upper boundary is known as a free surface. Flow nets in saturated systems bounded by a free surface can be constructed in the usual way, but there is one complication. The position of the entire free surface (not just the exit point) is unknown a priori. The boundary conditions on a free surface must satisfy both those of a water table h = z, and those of a flowline (equipotential lines must meet it at right angles). Its position is usually determined by graphical trial and error. Texts on engineering seepage, such as Harr (1962) or Cedergren (1967), provide hints on graphical construction and include many examples of steady-state, free-surface flow nets.

Figure 5.14(b) is the equivalent free-surface flow net to the saturated-unsaturated flow net shown in Figure 5.14(a). A glance at the two diagrams confirms that our decision to specify the water table as a flowline is quite a good approximation for this particular flow system. The outflow boundary ED is still known as a seepage face. We will encounter seepage faces in a practical sense when we examine hillslope hydrology (Section 6.5) and when we consider seepage through earth dams (Section 10.2).

Dupuit-Forchheimer Theory of Free-Surface Flow

For flow in unconfined systems bounded by a free surface, an approach pioneered by Dupuit (1863) and advanced by Forchheimer (1930) is often invoked. It is based on two assumptions: (1) flowlines are assumed to be horizontal and equipotentials vertical and (2) the hydraulic gradient is assumed to be equal to the slope of the free surface and to be invariant with depth. Figure 5.14(c) shows the equipotential net for the same problem as in Figure 5.14(a) but with the Dupuit assumptions in effect. The construction of rigorous flowlines is no longer possible. This paradoxical situation identifies Dupuit-Forchheimer theory for what it is, an empirical approximation to the actual flow field. In effect, the theory neglects the vertical flow components. In practice, its value lies in reducing the two-dimensional system to one dimension for the purposes of analysis. Calculations based on the Dupuit assumptions compare favorably with those based on more rigorous methods when the slope of the free surface is small and when the depth of the unconfined flow field is shallow. The discharge Q through a cross section of unit width perpendicular to the page in Figure 5.14(c) is given by

Q = Kh(x)\frac{dh}{dx} (5.28)

where h(x) is the elevation of the free surface above the base of the flow system at x, and the gradient dh/dx is given by the slope of the free surface Δhx at x. For steady-state flow, Q must be constant through the system and this can only be true if the free surface is a parabola.

The equation of flow for Dupuit-Forchheimer theory in a homogeneous, isotropic medium can be developed from the continuity relationship, dQ/dx = 0. From Eq. (5.28), this leads to

\frac{d^2(h^2)}{dx^2} = 0 (5.29)

If a three-dimensional unconfined flow field is reduced to a two-dimensional xy horizontal flow field by invocation of Dupuit-Forchheimer theory, the equation of flow in a homogeneous, isotropic medium becomes

\frac{\partial^2(h^2)}{\partial x^2} + \frac{\partial^2(h^2)}{\partial y^2} = 0 (5.30)

In other words, h2 rather than h must satisfy Laplace’s equation. It is possible to set up steady-state boundary-value problems based on Eq. (5.30) and to solve for h(x, y) in shallow, horizontal flow fields by analog or numerical simulation. It is also possible to develop a transient equation of flow for Dupuit free-surface flow in unconfined aquifers whereby h2 replaces h in the left-hand side of Eq. (2.77).

Harr (1962) discusses the practical aspects of Dupuit-Forchheimer theory in some detail. Bear (1972) includes a very lengthy theoretical treatment. Kirkham (1967) examines the paradoxes in the theory and provides some revealing explanations. The approach is widely used in engineering applications.

Suggested Readings

CEDERGREN, H. R. 1967. Seepage, Drainage and Flow Nets. Chapter 4: Flow Net Construction, John Wiley & Sons, New York, pp. 148-169.

HARR, M. E. 1962. Groundwater and Seepage. Chapter 2: Application of the Dupuit Theory of Unconfined Flow, McGraw-Hill, New York, pp. 40-61.

KIRKHAM, D. 1967. Explanation of paradoxes in Dupuit-Forchheimer seepage theory. Water Resources Res., 3, pp. 609-622.

PRICKETT, T. A. 1975. Modeling techniques for groundwater evaluation. Adv. Hydrosci., 10, pp. 42-45, 66-75.

REMSON, I., G. M. HORNBERGER, and F. J. MOLZ. 1971. Numerical Methods in Subsurface Hydrology. Chapter 4: Finite-Difference Methods Applied to Steady-Flow Problems, Wiley Interscience, New York, pp. 123-156.