Variably Saturated flow

Numerous numerical models have been developed for simulating water movements in materials with various degrees of saturation. Clement et al (1994) reviewed many of them underlining the most common shortcomings. In most applications, the pressure-based form of the variably saturated flow equation through homogeneous and isotropic porous media is used (Cooley, 1983; Huyakorn et al., 1984):
                                                                                        (1)
where:
S_s is the specific storage of the medium [L^-1]
Ψ is the pressure head [L]
θ is the water content
η is the porosity
K(Ψ) is the hydraulic conductivity [L T^-1]
t is time [T]
This equation describes the movement of water in soils under the following assumptions:
1. The dynamics of the air phase do not affect those of the water phase;
2. The density of water is only a function of pressure;
3. The spatial gradient of the water density is negligible.
Finite elements solutions to this approach often incur mass-balance problems in unsaturated media (Celia et al., 1990).
The finite-difference algorithm developed by Clement et al. has the advantages of being computationally simple while remaining capable of modelling a wide variety of problems, including infiltration into dry soils. In this approach, a mixed form of Equation 1 is solved by calculating the pressure heads Ψ at each node for a given time step (n+1) using a modified Picard iteration, where at each iteration (m+1) we solve a system of linear equations defined as:
                                                                                                                                                                                                                                   (2)
where
 
 
 
                                                                                                              (3)
where
Θ is the moisture content
K is the hydraulic conductivity
C is the water capacity
Δx is the horizontal step
Δz is the vertical step
Δt is the time step
S_s is the specific storage
η is the porosity
As can be seen, all the coefficients used for the current iteration (m+1) are calculated from the results of the previous iteration (m). Θ, K and C vary for each node and are functions of Ψ, therefore need to be re-calculated after each iteration.
As discussed earlier, the link between Θ and Ψ is given by the Soil-Water Characteristic Curve (SWCC). The hydraulic conductivity of the materials can be expressed as a function of the pressure heads (Van Genuchten, 1980) by:
                                                                                               (4)
Finally, by definition, the water capacity represents the slope of the SWCC and therefore can be expressed as:
                                                                                                                                               (5)
Equation 2 applies to all interior nodes, while at boundary nodes it needs to be modified to consider the appropriate boundary conditions.
In the examples shown hereafter we consider a pavement structure like the one in Figure 1, where the water table depth is 2m from the top of the unbound layer. In normal conditions this system is in equilibrium and the water content does not vary with time.

Figure 1: Variably saturated flow – pavement structure
The boundaries considered are as follows (see also Figure 2):
► Top and right-hand side nodes: Neumann boundaries (no flow);
► Bottom nodes: Dirichlet boundaries (constant moisture content);
► Left-hand side nodes: symmetry.
Figure 2: Variably saturated flow – boundary conditions

Posted in: model, moisture