**2. Methods of water flow analysis**

#### **2.1. Analytical solutions**

#### *2.1.1. Exact solution*

An exact analytical solution is generally only possible to obtain when the geometry of the flux domain and the hydraulic and boundary conditions are simple (isotropic and homogeneous media). However, soil is a heterogeneous, anisotropic, and discontinuous medium that has different characteristics at each point. Thus, exact solutions are difficult and impractical to obtain in the case of the most complex flow problems, such as those in practical geotechnical engineering. In consequence, approximate solutions are usually sought.

#### *2.1.2. Conformal transformation or mapping*

One of the approximate analytical solutions that is used to solve two-dimensional water flow problems involves obtaining a function that can transform the problem from the complex geometric domain into a problem whose solution is known [7]. These functions transform geometric shapes from a complex plane *ω* into shapes in another complex plane *ζ*. Thus, the mapping procedure defines the correspondence of the points of one shape in a plane *ω* to the points of the respective figure in the plane *ζ*. A transformation or mapping function is called conformal (conformal mapping) when it does not change the angles of intersection or the approximate geometric shapes between the two planes of interest [8]. Based on this concept, the Laplace differential equation can be solved for a domain *G* (**Figure 2b**) if the transformation or conformal mapping of this domain with a simpler domain *G*1 (**Figure 2a**) is known [8–10]. The transformation is carried out by means of the analytical function of a complex variable. One of the best known functions that transform a system of uniform flux in the *ω* plane (**Figure 2a**) into a system of flux with confocal parabolas in the *ζ* plane (**Figure 2b**) is as follows (representing the well-known Kozeny's solution for unconfined flow through dams, **Figure 3**):

$$
\zeta = \alpha^2 \tag{6}
$$

**Figure 2.** Conformal mapping [7].

**Figure 3.** Kozeny's solution for water flow through an earth dam with a horizontal filter [7].

#### *2.1.3. Method of fragments*

obtain in the case of the most complex flow problems, such as those in practical geotechnical

One of the approximate analytical solutions that is used to solve two-dimensional water flow problems involves obtaining a function that can transform the problem from the complex geometric domain into a problem whose solution is known [7]. These functions transform geometric shapes from a complex plane *ω* into shapes in another complex plane *ζ*. Thus, the mapping procedure defines the correspondence of the points of one shape in a plane *ω* to the points of the respective figure in the plane *ζ*. A transformation or mapping function is called conformal (conformal mapping) when it does not change the angles of intersection or the approximate geometric shapes between the two planes of interest [8]. Based on this concept, the Laplace differential equation can be solved for a domain *G* (**Figure 2b**) if the transformation or conformal mapping of this domain with a simpler domain *G*1 (**Figure 2a**) is known [8–10]. The transformation is carried out by means of the analytical function of a complex variable. One of the best known functions that transform a system of uniform flux in the *ω* plane (**Figure 2a**) into a system of flux with confocal parabolas in the *ζ* plane (**Figure 2b**) is as follows (representing the well-known Kozeny's solution for unconfined flow through dams, **Figure 3**):

2

= (6)

z w

engineering. In consequence, approximate solutions are usually sought.

*2.1.2. Conformal transformation or mapping*

94 Groundwater - Contaminant and Resource Management

**Figure 2.** Conformal mapping [7].

The method of fragments is an approximate analytical method of solution for confined flow domains of finite depth. The fundamental assumption is that equipotential lines at different parts of the flow region can be approximated by straight vertical lines that divide the flow region into sections or fragments (**Figure 4**) [8, 11]. This method requires a form factor Φ that is obtained by solving definite integrals set up for each fragment of the flow region. Tables of expressions for different form factors for typical confined flow problems have been developed [8, 11]. The equation for calculating the discharge *q* through all fragments using this method is as follows:

$$q = k \frac{\sum \Delta h\_n}{\sum \Phi\_n} = \frac{k \cdot \Delta h}{\sum\_{n=1}^{n} \Phi\_n} \tag{7}$$

where *k* = hydraulic conductivity of the homogeneous and isotropic medium, Δ*hn* = loss of hydraulic head through fragment *n*, Δ*h* = total loss of hydraulic head, and Φ*n* = dimension‐ less form factor in fragment *n*.

**Figure 4.** Schematic representation of the method of fragments.

The method of fragments was originally proposed to study confined flow in homogeneous and isotropic media, but at present, it has been implemented in anisotropic media [12]. In combination with conformal mapping and the Kozeny's parabola, it has also been applied to solve problems of unconfined flow through homogeneous levees with horizontal filters [13].

#### **2.2. Graphical solutions**

#### *2.2.1. Confined flow*

The most popular approximation technique for solving water flow problems is known as the flow net. It is a graphical method that sets up two functions that satisfy the Laplace's equation and that geometrically constitute two families of orthogonal curves: (a) equipotential lines (constant potential *ϕ*) and (b) flow lines or streamlines (constant values of the stream function *ψ*). The graphical representation of these lines is the so-called flow net. A drawing that satisfies the boundary and orthogonality conditions allows water flow problems with homogeneous and isotropic soil to be solved simply and graphically. The expression that calculates the discharge or rate of seepage *q* using a flow net is as follows:

$$q = k \,\Delta h (\frac{n\_f}{n\_e}) = k \,\Delta h \,\text{S} \tag{8}$$

where *k* = hydraulic conductivity of the homogeneous and isotropic medium, Δ*h* = total loss of hydraulic head, *\$* = *nf* /*ne* is the form factor, *nf* = number of flow intervals or flow channels, and *ne* = number of equipotential intervals.

Flow nets are typically drawn on paper. However, nowadays, it is possible to draw flow nets by computer using several numerical techniques based on finite element method (FEM) or finite difference equations [14], such as successive over-relaxation (SOR) [15]. Examples of flow nets that were drawn using these techniques in confined homogeneous and stratified domains are shown in **Figures 5(a and b)**–**7**.

#### *2.2.2. Unconfined flow*

Free surface problems involve boundary value problems in which a portion of the boundary, the free surface, is unknown and must be determined as part of the solution. The presence of the free surface or water table makes the analysis methods more difficult. Dupuit's parabola [16] and Kozeny's parabola [17] are rigorous solutions for drawing the upper flow line and are only applicable for homogenous and isotropic media with specific geometries, such as vertical walls (Dupuit) or with filters (Kozeny). Other approximated methods, such as the *tangent* [18, 19] and *sine* methods [20], allow mainly calculate the discharge point of the upper flow line. Currently, this can be determined using numerical methods such as finite element method (FEM) and finite differences (FD), among others. The Baiocchi's method [21] and the extended pressure technique [22] are two variants of the successive over-relaxation (SOR) method (based on algebraic finite difference equations) that can be used to determine the position of the upper flow line in homogeneous media or media composed of materials with different permeability values, respectively.

The method of fragments was originally proposed to study confined flow in homogeneous and isotropic media, but at present, it has been implemented in anisotropic media [12]. In combination with conformal mapping and the Kozeny's parabola, it has also been applied to solve problems of unconfined flow through homogeneous levees with horizontal filters [13].

The most popular approximation technique for solving water flow problems is known as the flow net. It is a graphical method that sets up two functions that satisfy the Laplace's equation and that geometrically constitute two families of orthogonal curves: (a) equipotential lines (constant potential *ϕ*) and (b) flow lines or streamlines (constant values of the stream function *ψ*). The graphical representation of these lines is the so-called flow net. A drawing that satisfies the boundary and orthogonality conditions allows water flow problems with homogeneous and isotropic soil to be solved simply and graphically. The expression that calculates the

> () \$ *<sup>f</sup> e*

where *k* = hydraulic conductivity of the homogeneous and isotropic medium, Δ*h* = total loss

Flow nets are typically drawn on paper. However, nowadays, it is possible to draw flow nets by computer using several numerical techniques based on finite element method (FEM) or finite difference equations [14], such as successive over-relaxation (SOR) [15]. Examples of flow nets that were drawn using these techniques in confined homogeneous and stratified domains

Free surface problems involve boundary value problems in which a portion of the boundary, the free surface, is unknown and must be determined as part of the solution. The presence of the free surface or water table makes the analysis methods more difficult. Dupuit's parabola [16] and Kozeny's parabola [17] are rigorous solutions for drawing the upper flow line and are only applicable for homogenous and isotropic media with specific geometries, such as vertical walls (Dupuit) or with filters (Kozeny). Other approximated methods, such as the *tangent* [18, 19] and *sine* methods [20], allow mainly calculate the discharge point of the upper flow line. Currently, this can be determined using numerical methods such as finite element method (FEM) and finite differences (FD), among others. The Baiocchi's method [21] and the extended pressure technique [22] are two variants of the successive over-relaxation (SOR) method (based on algebraic finite difference equations) that can be used to determine the position of the upper

*<sup>n</sup> q kh kh <sup>n</sup>* =D =D (8)

/*ne* is the form factor, *nf* = number of flow intervals or flow channels,

discharge or rate of seepage *q* using a flow net is as follows:

**2.2. Graphical solutions**

96 Groundwater - Contaminant and Resource Management

of hydraulic head, *\$* = *nf*

*2.2.2. Unconfined flow*

and *ne* = number of equipotential intervals.

are shown in **Figures 5(a and b)**–**7**.

*2.2.1. Confined flow*

**Figure 5.** Flow nets in homogeneous and isotropic media obtained using the method of successive over-relaxations (SOR) [23].

**Figure 6.** Flow net in a homogeneous and isotropic medium obtained using the finite element method (FEM) [24].

**Figure 7.** Flow net in a stratified soil under an impermeable dam obtained using the method of successive over-relaxa‐ tions (SOR) [23].

A simple graphical procedure to draw the Kozeny's parabola or any parabolic upper flow line in a homogeneous and isotropic medium is as follows (**Figure 8a**): (a) The distance *a0* is calculated with the formula *a0* = *y0/2* = *[(d2* + *h2 ) 1/2-d]/2*; (b)Draw a vertical line through *O* and also a horizontal line through *M*; (c) Divide the *OB* segment in a number of equal parts, and the *MB* segment must also be divided into the same number of equal parts; (d) Draw straight lines joining the point *O* with the divisions made in the *MB* segment; (e) Draw horizontal lines passing through the divisions made in the *OB* segment; (f) The intersections of the previous lines are the points of the parabola sought. **Figures 8(b)** and **9** show upper flow lines obtained using this graphical procedure and other numerical methods. Additionally, **Figures 10** and **11** show flow nets that were numerically calculated in these types of free surface problems.

**Figure 8.** Homogeneous and isotropic embankment: (a) graphical procedure to draw the upper flow line, (b) compari‐ son of the upper flow lines obtained using different methods [14].

**Figure 9.** Comparison of the upper flow lines in a homogeneous and isotropic embankment, obtained using different methods [14].

**Figure 10.** Flow net in a homogeneous and isotropic embankment numerically calculated using the extended pressure technique [14].

**Figure 11.** Flow net in a dam with graded materials obtained using the extended pressure technique [23].

#### **2.3. Numerical solutions (approximate solutions)**

#### *2.3.1. Finite elements*

**Figure 8.** Homogeneous and isotropic embankment: (a) graphical procedure to draw the upper flow line, (b) compari‐

**Figure 9.** Comparison of the upper flow lines in a homogeneous and isotropic embankment, obtained using different

**Figure 10.** Flow net in a homogeneous and isotropic embankment numerically calculated using the extended pressure

son of the upper flow lines obtained using different methods [14].

98 Groundwater - Contaminant and Resource Management

methods [14].

technique [14].

The finite element method (FEM) is a numerical technique that provides approximate solutions to partial differential equations to solve a problem of a particular field. It is more versatile than other methods because it can consider anisotropy, heterogeneity, and multiple boundary conditions. The 2D finite elements that are generally used in water flow problems are triangles [25] or a combination of triangles and squares [26] whose nodes coincide with their vertices; in addition, triangles (2D) and tetrahedra (3D) can be used [27]. The hydraulic head is assumed to vary linearly within each finite element, and the Laplace's equation can be solved using a variational approach. Thus, the solution to this equation in a domain is found by obtaining the minimum of a function that is related to the equation and is defined for this domain. Based on these assumptions and after several mathematical manipulations, the following systems of homogeneous linear equations are set up:

$$\{S\}\{h\_r\} = 0\tag{9}$$

$$\{S\}\{\psi\_r\} = 0\tag{10}$$

The solution to Eq. (9) using a known method, such as Gaussian elimination, helps to determine the hydraulic head *h* at the nodes in the mesh of finite elements where it is unknown. Similarly, the solution to Eq. (10) provides nodal values of the stream function *ψ*. The flow net of the problem can be obtained by drawing the isovalue curves for this pair of families. Likewise, once the hydraulic heads *h* are calculated with Eq. (9), other results for the water flow problem can be obtained, such as the hydraulic gradients, flow velocities, pore pressure, degree of saturation, and flow rate, among others.

#### *2.3.2. Finite differences*

The Laplace's equation can be solved using finite difference equations, which are the same as those developed via truncated Taylor series or directly from Darcy's Law [4]. Several methods can be used to evaluate water flow problems that utilize finite differences, including: (a) the classical method of relaxations, (b) the method of successive over-relaxations, and (c) random walks.

#### *2.3.2.1. Classical method of relaxations*

The classical relaxation method is an iterative process in which solutions for water flow through porous media can be obtained by simply knowing the domain geometry and hy‐ draulic boundary conditions. This method can solve Laplace's equation for a point (node) relative to its surrounding points using an algebraic finite difference equation. For this procedure, a square mesh with dimensions of *δx* = *δy* is drawn in the flow zone if the medium is homogenous and isotropic (similar to those in **Figure 11a** and **b**), and a rectangular mesh with dimensions of *δ<sup>x</sup>* ≠ *δy* is drawn if the medium is anisotropic. The intersections of the squares or rectangles constitute the nodes of the mesh. For these nodes, approximate values of the hydraulic head or potential *h* (points where *h* requires to be calculated) must be assigned while respecting the known values of *h* in the flow boundaries. These values usually correspond to the upper and lower water levels or the upstream water level and downstream water level of the problem at hand. The values assigned in the nodes are arbitrary and can be zero or the result of a reasonable estimation. However, although there are several techniques that can be used to ensure that the value of the potential imposed on the nodes where *h* is not known is as accurate as possible (Young, 1950), it is important to verify the precision of the assigned data manually by calculating the residue in each node [28]. For example, the difference between the hydraulic potential of the four surrounding nodes is calculated with regard to the central or interior node and so on. Therefore, the relaxation procedure involves the systematic refinement of this residue throughout the grid until the residue in all mesh nodes of interest is zero or practically zero. This value indicates that the Laplace's equation in the study domain has been fulfilled, and therefore, the flow problem has been solved for a certain water flow problem. A disadvantage of the method of relaxations is that it is based on assigning arbitrary values to the nodes in the study mesh, which makes it difficult for the residuals to equal zero at an early stage of calculation. As a result, additional steps of reassigning values are generally necessary, which makes the method long and laborious in practice.

#### *2.3.2.2. Technique of successive over-relaxation*

The technique of successive over-relaxation (SOR) [15] is a modification of the classical method of relaxations in which the process of residual refinement at the nodes is automatic because it utilizes the Gauss-Seidel iterative method (**Figure 12**); this makes it possible to obtain residuals of zero or nearly zero at all the nodes, which allows water flow problems to be solved relatively quickly and easily [29, 30]. An additional significant advantage of the SOR method is that it can solve the so-called free surface problems (or unconfined flow problems), in which the position of the free (or phreatic) surface must be determined to solve the problem. Other improvements to the SOR method have been developed for this type of problem, including (a) the Baiocchi's solution [21] and (b) the extended pressure method [31, 22]. The former method helps to determine the position of the phreatic surface (upper flow line for steady-state flow) in a homogeneous medium (**Figure 10**), and the latter method helps to determine this line in both homogeneous and heterogeneous media (**Figure 11**) [14].

**Figure 12.** Schematic arrangement of nodes in the dam of graded materials in **Figure 11** using the extended pressure method [23].

#### *2.3.3. Random walks*

classical method of relaxations, (b) the method of successive over-relaxations, and (c) random

The classical relaxation method is an iterative process in which solutions for water flow through porous media can be obtained by simply knowing the domain geometry and hy‐ draulic boundary conditions. This method can solve Laplace's equation for a point (node) relative to its surrounding points using an algebraic finite difference equation. For this procedure, a square mesh with dimensions of *δx* = *δy* is drawn in the flow zone if the medium is homogenous and isotropic (similar to those in **Figure 11a** and **b**), and a rectangular mesh with dimensions of *δ<sup>x</sup>* ≠ *δy* is drawn if the medium is anisotropic. The intersections of the squares or rectangles constitute the nodes of the mesh. For these nodes, approximate values of the hydraulic head or potential *h* (points where *h* requires to be calculated) must be assigned while respecting the known values of *h* in the flow boundaries. These values usually correspond to the upper and lower water levels or the upstream water level and downstream water level of the problem at hand. The values assigned in the nodes are arbitrary and can be zero or the result of a reasonable estimation. However, although there are several techniques that can be used to ensure that the value of the potential imposed on the nodes where *h* is not known is as accurate as possible (Young, 1950), it is important to verify the precision of the assigned data manually by calculating the residue in each node [28]. For example, the difference between the hydraulic potential of the four surrounding nodes is calculated with regard to the central or interior node and so on. Therefore, the relaxation procedure involves the systematic refinement of this residue throughout the grid until the residue in all mesh nodes of interest is zero or practically zero. This value indicates that the Laplace's equation in the study domain has been fulfilled, and therefore, the flow problem has been solved for a certain water flow problem. A disadvantage of the method of relaxations is that it is based on assigning arbitrary values to the nodes in the study mesh, which makes it difficult for the residuals to equal zero at an early stage of calculation. As a result, additional steps of reassigning values are generally

necessary, which makes the method long and laborious in practice.

The technique of successive over-relaxation (SOR) [15] is a modification of the classical method of relaxations in which the process of residual refinement at the nodes is automatic because it utilizes the Gauss-Seidel iterative method (**Figure 12**); this makes it possible to obtain residuals of zero or nearly zero at all the nodes, which allows water flow problems to be solved relatively quickly and easily [29, 30]. An additional significant advantage of the SOR method is that it can solve the so-called free surface problems (or unconfined flow problems), in which the position of the free (or phreatic) surface must be determined to solve the problem. Other improvements to the SOR method have been developed for this type of problem, including (a) the Baiocchi's solution [21] and (b) the extended pressure method [31, 22]. The former method helps to determine the position of the phreatic surface (upper flow line for steady-state flow)

*2.3.2.2. Technique of successive over-relaxation*

walks.

*2.3.2.1. Classical method of relaxations*

100 Groundwater - Contaminant and Resource Management

The random walk method (RWM) consists of studying the movements of a particle that travels in a random way over the nodes of a mesh of a flow domain (**Figure 13a** and **b**), which allows the hydraulic head to be determined at points of interest by numerical solution of the Laplace's equation in terms of finite differences. This method relies on the so-called Monte Carlo techniques [32, 33], which are an alternative to the usual methods of water flow analysis. Specifically, the method generates a series of random trajectories (via random numbers) that start from node *p*<sup>0</sup> in the mesh. Thus, the particle moves randomly through the interior nodes of the mesh and stops when it reaches a boundary node, which is called an absorbent node, because the value of the hydraulic head at that node is known (imposed boundary condition). A complete trajectory is made up of a sequence of nodes, and the last node is an absorbent node. The hydraulic head is then determined by counting the number of trajectories that end at different boundaries, multiplying them by the value of the hydraulic head at the respective boundary and dividing the result by the total number of trajectories. This procedure is repeated several times, and the results are an unbiased measure of the hydraulic head at the node of interest:

$$h\_0 = \frac{n\_1 f\_1 + n\_2 f\_2}{n\_1 + n\_2} \tag{11}$$

where *h*0 = hydraulic head calculated at point *p*<sup>0</sup> (**Figure 13b**), *f*1 and *f*<sup>2</sup> = boundaries with known hydraulic heads, and *n*1 and *n*2 = number of trajectories that reach boundaries 1 and 2, respectively.

**Figure 13.** 2D and 3D meshes with which the homogeneous and isotropic flow regions are modelled by the random walk method [24].

The RWM has been utilized to solve confined water flow problems [33] and also to calculate the equivalent permeability *k*equivalent in simulated heterogeneous media (**Figure 14**) [24]:


**Figure 14.** Example of a random walk with the PASECA-2003 algorithm [24].

#### **2.4. Stochastic solutions**

The uncertainty due to the spatial variation of permeability is the most important factor that must be taken into account in the analysis of water flow through soils. **Figure 15** summarizes the main techniques that are used to evaluate the propagation of this uncertainty. It is common to utilize probabilistic techniques in combination with numerical methods, such as finite elements (FEM), finite differences (FDM), integral equations (BEM), and random walks (RWM). A stochastic analysis permits a more realistic evaluation of water flow problems and can be useful in defining zones of uncertainty, which can be used to identify the parts of the flow region that are prone to significant variations in properties such as the hydraulic head and hydraulic gradient, as is illustrated in **Figures 16** and **17**, respectively, [24]

**Figure 15.** A summary of the main techniques that are used to evaluate uncertainty in water flow analyses [24].

**Figure 13.** 2D and 3D meshes with which the homogeneous and isotropic flow regions are modelled by the random

The RWM has been utilized to solve confined water flow problems [33] and also to calculate the equivalent permeability *k*equivalent in simulated heterogeneous media (**Figure 14**) [24]:

The uncertainty due to the spatial variation of permeability is the most important factor that must be taken into account in the analysis of water flow through soils. **Figure 15** summarizes the main techniques that are used to evaluate the propagation of this uncertainty. It is common to utilize probabilistic techniques in combination with numerical methods, such as finite elements (FEM), finite differences (FDM), integral equations (BEM), and random walks (RWM). A stochastic analysis permits a more realistic evaluation of water flow problems and can be useful in defining zones of uncertainty, which can be used to identify the parts of the

**–** 1D water flow analysis→ *k*equivalent is the harmonic mean **–** 2D water flow analysis→ *k*equivalent is the geometric mean

102 Groundwater - Contaminant and Resource Management

**–** 3D water flow analysis→ *k*equivalent tends to the arithmetic mean

**Figure 14.** Example of a random walk with the PASECA-2003 algorithm [24].

walk method [24].

**2.4. Stochastic solutions**

**Figure 16.** Standard deviation of the hydraulic head (m) under a sheet pile in a medium of two stratified isotropic ma‐ terials [24].

**Figure 17.** Standard deviation of the (dimensionless) hydraulic gradient under a spillway structure in a medium of two stratified isotropic materials [24].
