## 1. Theoretical foundations

### 1.1. General considerations

Flow through saturated or unsaturated soils is governed by Darcy’s Law, which was originally proposed for saturated media. Research has demonstrated that this law is also applicable to the flow of water in unsaturated soils [1]. The main difference between these flows is that the hydraulic conductivity for saturated media is a constant value, but it varies as a function of the volumetric water content and also indirectly with changes in pore water pressure in unsaturated soils (**Figure 1**) [2, 3]. Darcy’s Law [4] is often written as follows:

where *v* = Darcy’s velocity, *k* = hydraulic conductivity, and *i* = total hydraulic head gradient.

The average velocity at which the water moves through a mass of soil is linear and is equal to the Darcy’s velocity divided by the porosity of the soil. In an unsaturated soil, the average velocity is equal to the Darcy’s velocity divided by the volumetric water content of the soil. The majority of analytical and numerical methods that are currently employed for solving water flow problems consider only the Darcy’s velocity.

### 1.2. Equation for steady-state flow (saturated porous media)

The equation that describes steady-state flow in a porous medium is based on Darcy’s Law [4] and on the principle of flow continuity (which states that the amount of water that enters the medium is equal to the amount that exits) and is known as the Laplace’s equation (for a homogeneous and isotropic medium with *k _{x}* =

*k*=

_{y}*k*):

_{z}where *h* = total hydraulic head, and *k _{x}*,

*k*, and

_{y}*k*= hydraulic conductivities in the

_{z}*x*-,

*y*-, and

*z*-directions, respectively.

The Laplace’s equation is satisfied under the following conditions: (a) the flow is steady-state, (b) the soil is saturated, (c) the water and the solid particles are incompressible, (d) the flow does not modify the soil structure, and (e) there are no sources (via injection or extraction of water).

### 1.3. Equation for transient flow (saturated or unsaturated porous media)

In transient flow, the water levels vary as a function of time, and thus, water is either stored in or discharged from the medium. In these cases:

or:

Equation (3) refers to the case of water drawdown, and Eq. (4) refers to the case of water filling. All the previous assumptions lead to the general mass balance equation:

where *h* = total hydraulic head, *k _{x}* = hydraulic conductivity in the

*x*-direction,

*k*= hydraulic conductivity in the

_{y}*y*-direction,

*Q*= source term (applied boundary flux: injection or extraction),

*θ*= volumetric water content, and

*t*= time.

Equation (5) is the so-called Richards’s equation, and it describes transient flow in unsaturated soils. The term on the right (*∂θ*/*∂t*, the rate of change of the volumetric water content with respect to time) is related to the change in the water level with time. When there is no variation with time (*∂θ*/*∂t* = 0) and no source term (*Q* = 0), Eq. (5) becomes the Laplace’s equation, which is used for steady-state flow in saturated soils. The Richards’s equation, or any of its modified forms, has constituted the basis for the development of most numerical models to calculate infiltration through unsaturated porous media under transient-state flow conditions [6].

## 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**):

#### 2.1.3. Method of fragments

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:

where *k* = hydraulic conductivity of the homogeneous and isotropic medium, Δ*h _{n}* = loss of hydraulic head through fragment

*n*, Δ

*h*= total loss of hydraulic head, and Φ

*n*= dimensionless form factor in fragment

*n*.

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:

where *k* = hydraulic conductivity of the homogeneous and isotropic medium, Δ*h* = total loss of hydraulic head, *$* = *n _{f}*/

*n*is the form factor,

_{e}*n*= number of flow intervals or flow channels, and

_{f}*n*= number of equipotential intervals.

_{e}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.

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 *a _{0}* is calculated with the formula

*a*=

_{0}*y*=

_{0}/2*[(d*+

^{2}*h*; (b)Draw a vertical line through

^{2})^{1/2}-d]/2*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.

### 2.3. Numerical solutions (approximate solutions)

#### 2.3.1. Finite elements

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:

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 hydraulic 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}* =

*δ*is drawn in the flow zone if the medium is homogenous and isotropic (similar to those in

_{y}**Figure 11a**and

**b**), and a rectangular mesh with dimensions of

*δ*≠

_{x}*δ*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

_{y}*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].

#### 2.3.3. Random walks

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*_{0} 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:

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

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]:

### 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]

## 3. Computer programs and general recommendations for different types of analysis

Numerical techniques are currently preferred because of their ability to solve complex problems in which Eqs. (2) and (5) can be generalized to heterogeneous media with anisotropic materials and boundary conditions of variable complexity [34, 35]. The general methodology for performing a steady- or transient-state water flow analysis is illustrated in **Figure 18**, which shows that some of the most important data for performing water flow analysis are the hydraulic parameters of the materials, which must be obtained from field or laboratory tests that unfortunately are not always carried out because of time or cost. In addition, some tests require specialized personnel and equipment, such as the tests to determine the hydraulic parameters of unsaturated materials, including the soil–water characteristic curve (SWCC). The necessary hydraulic functions of the soil for analyzing unsaturated soils are as follows:

–

*The water retention curve*(**Figure 27**) is also known as the*soil–water characteristic curve*(SWCC) depending on whether the suction is expressed in terms of the degree of saturation or the volumetric water content, respectively. The SWCC is broadly defined as the relationship between the amount of water in the soil and soil suction.–

*The hydraulic conductivity function*(**Figure 28**) represents the suction as a function of the permeability.

**Figures 19** and **20** show several general considerations for different types of water flow analyses, which can be performed relatively easily and rapidly using any of the existing specialized algorithms. **Figure 21** summarizes some of the most popular programs that are used to numerically solve water flow problems. One benefit of computer programs is their ability to facilitate the study of transient flow and unsaturated soil conditions, which are difficult and laborious to solve analytically.

Several programs for water flow analysis consider soil classification systems that are different from the Unified Soil Classification System (USCS, which is commonly used in geology, soil mechanics, and geotechnical engineering) because they involve parameters that are used to study unsaturated soils, such as:

– HYPRES database = Hydraulic Properties of European Soils

– USDA = United States Department of Agriculture

– STARING = Dutch ‘Winand Staring Soil Series’

These systems imply that it is not advisable to use only the standard parameters that the computer programs include for certain types of materials (e.g., sand, clay, silt) because of the variations in the characteristics of soils (European, USA, or Dutch). It is preferable to assign the necessary hydraulic parameters based on the type of analysis and use values that are obtained from laboratory or field tests of the materials of the earth structure or soil that is being studied. In recent years, a comprehensive database that contains the hydraulic parameters of different types of soils from around the world has been developed [36]. Additionally, some algorithms [27] use granulometric curves as well as the index properties of the materials in the flow region and various mathematical expressions to estimate the hydraulic functions that are needed for the analyses (**Figures 27** and **28**). Some of the main mathematical models that are used to obtain the soil hydraulic functions are as follows[37]:

**Figures 25**–**31** show a case of the analysis of water flow through a cofferdam for *La Yesca* dam in Mexico composed of graduated materials assuming that the soil in one part of the cofferdam is partially saturated [37].

The PLAXFLOW algorithm [25] solves transient-state flow problems using the finite element method (FEM) by means of an approximate solution to Eq. (5) and by representing the flow in unsaturated soils with the *Van Genuchten* model. This algorithm performs analyses of transient-state flow in two different ways (**Figure 22**): (a) with step-wise conditions, in which each phase is defined by constant boundary conditions; that is, each time period is associated with a certain water level and (b) with time-dependent conditions, in which the continuous variation of the water level is explicitly considered as a function of time, which can be represented by linear functions, harmonic functions, or data in tables.

**Figures 23** and **24** show the results of water drawdown analyses that were carried out using the SEEP/W [26] and PLAXFLOW [25] algorithms, respectively [35, 43].

## 4. Conclusions

An exact analytical solution of a water flow problem 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). 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. Numerical techniques are currently preferred over other methods due to their ability to solve complex problems in which the equations for water flow analyses can be generalized to heterogeneous media with anisotropic materials and boundary conditions of varying complexity. However, more sophisticated analyses require the use of a greater number of material parameters, which involves laboratory and field testing that require specialized knowledge and personnel. A series of mathematical models is available to approximate these material parameters. However, the results of numerical analysis should be compared to measurements from monitoring (field instrumentation) of the structures being studied.

Some important comments about the numerical analyses of water flow are as follows:

– In 3D groundwater flow analyses, the use of finite difference equations requires less computational time than the 3D finite element method.

– Computer programs facilitate analyses of transient flow and unsaturated soil conditions, which are difficult and laborious to solve analytically.

– Computer programs cannot replace the good judgment of an engineer.