Variably Saturated Flow

This example uses the Richards’ Equation interface to assess how well geophysical irrigation sensors detect the true level of fluid saturation in variably saturated soils. Andrew Hinnell, Alex Furman, and Ty Ferre from the Department of Hydrology and Water Resources at the University of Arizona brought the example to us. They originally worked out the problem in COMSOL Multiphysics’ PDE interfaces, but this discussion shares their elegant model reformulated in the Richards’ Equation interface.

A major challenge when characterizing fluid movement in variably saturated porous media lies primarily in the need to describe how the capacity to transmit and store fluids changes as fluids enter and fill the pore spaces. Experimental data for these properties are difficult to obtain. Moreover, the properties that change value as the soil saturates happen to be equation coefficients, which makes the mathematics notoriously nonlinear. The Richards’ Equation interface provides interfaces that automate the van Genuchten (Ref. 1) as well as the Brooks and Corey (Ref. 2) relationships for fluid retention and material properties that vary with the solution.

Figure 1: A variably saturated porous medium.

This example uses the model of Hinnell, Furman, and Ferre to characterize how the distribution of water changes around three impermeable sensors inserted into two different blocks of uniform soil partially saturated with water. The question for the model to address is this: Does the saturation localized around the sensor give a valid picture of the saturation within the total block?

This example demonstrates how to use the Richards’ Equation interface including the van Genuchten (Ref. 1) as well as the Brooks and Corey (Ref. 2) retention models. The step-by-step instructions assist you in learning how to set up two separate equation systems in one COMSOL Multiphysics model file and then solve them simultaneously.

Model Definition

The problem setup is as follows. Two homogeneous columns of soil, each 2 m-by- 2 m on a face, are partially saturated with water. A plot of the hydraulic properties of the first soil (Soil Type 1) fits the van Genuchten retention and permeability formulas. The other soil (Soil Type 2) has material properties that suit the Brooks and Corey formulas. Within each soil column are three impermeable rods, each with a 0.1 m radius. The rods are spaced at 0.5 m increments so they run horizontally down the centerline of each block; see Figure 2. Just after the rods are emplaced, the pressure head is still uniform, but water begins to move vertically downward in steady drainage. Because all vertical slices down a block are identical, you can model a 2D cross section and observe the changes in the flow field for 900 s or 15 minutes.

Figure 2: Soil block with three rods. The shaded plane represents a vertical cross section.

Governing Equation

Richards’ equation describes the unsaturated-saturated flow of water in the soils. In this problem you can only use Richards’ equation for the water since the air in the soil is supposed to be at atmospheric pressure. The governing equation for the model is

Pressure head, Hp (m), is the dependent variable. C denotes specific moisture capacity (

), Se is the effective saturation, S is a storage coefficient (m−1), t is time, K denotes the hydraulic conductivity (m/s), and D is the coordinate (for example x, y, or z) for the vertical elevation (m). The equation does not show the volumetric fraction of water, θ, which is a constitutive relation that depends on Hp. Nonlinearities appear because C, Se, and K change with Hp and θ.

The first term in the equation explains that fluid storage can change with time during both unsaturated and saturated conditions. When the soil is unsaturated, the pores fill with (or drain) water. After the pore spaces completely fill, there is slight compression of the fluid and the pore space. The specific moisture capacity C = ∂θ/∂Hp describes the change in fluid volume fraction θ with pressure head. The storage coefficient addresses storage changes due to compression and expansion of the pore spaces and the water when the soil is fully wet. To model the storage coefficient, this example uses the specific storage option, which sets S = ρf g(χp + θχf). Here, ρf is the fluid density (kg/m3), g is the acceleration of gravity, while χp and χf are the compressibilities of the solid particles and fluid, respectively (m·s2/kg).

The van Genuchten and the Brooks and Corey formulas that describe the change in C, Se, K, and θ with Hp require data for the saturated and liquid volume fractions, θs and θr, as well as for other constants α, n, m, and l, which specify a particular type of medium. With the van Genuchten equations that follow, you consider the soil as being saturated when fluid pressure is atmospheric (that is, Hp = 0):

With the Brooks and Corey approach, an air-entry pressure distinguishes saturated (Hp > −1/α) and unsaturated (Hp < −1/α) soil as in

To find a unique solution to this problem, you must specify initial and boundary conditions. Initially, the column has uniform pressure head of Hp0. The boundary conditions are

where n is the unit vector normal to the boundary.

Model Data

The following table gives the data needed to complete the two example problems:


Variable	Unit	Description	van Genuchten	Brooks & Corey
g	m/s2	Gravity	9.82	9.82
ρf	kg/m3	Fluid density	1000	1000
χp	m·s2/kg	Compressibility solid particles	10-8	10-8
χf	m·s2/kg	Compressibility of fluid	4.4·10-10	4.4·10-10
Ks	m/s	Saturated hydraulic conductivity	8.25·10-5	5.83·10-5
θs		Porosity/void fraction	0.43	0.417
θr		Residual saturation	0.045	0.02
α	m-1	alpha parameter	14.5	13.8
n		n parameter	2.68	0.592
m		m parameter	1-1/n	n/a
l		Pore connectivity parameter	0.5	1
Hp0	m	Specified pressure	-0.06	-0.2
Hp0	m	Initial pressure	-0.06	-0.2

Results and Discussion

Figure 3 and Figure 4 are solutions to the Richards’ equation problem of Prof. Ty Ferre, Andrew Hinnell, and Alex Furman from the University of Arizona’s Department of Hydrology and Water Resources. Each figure gives results for similar variably saturated flow problem posed for different soil types. Each snapshot shows effective fluid saturation (surface plot), pressure head (contours), and fluid velocities (arrows). The flow field varies around the rods but remains largely uniform over the remainder of each block.

Figure 3: Solution for effective saturation (surface plot), pressure head (contours), and velocity (arrows) at 15 minutes for Soil Type 1 (van Genuchten).

Figure 4: Solution for effective saturation (surface plot), pressure head (contours), and velocity (arrows) at 15 minutes for Soil Type 2 (Brooks and Corey).

Figure 5 shows the effective saturation evolving over time at the rod-soil boundary. The intervals (−180°, 0°) and (0°, 180°) for the angular coordinate along the horizontal axis correspond to the boundary’s lower and upper halves, respectively. In the figure, the solid lines denote the solution for Soil Type 1, and the dashed lines correspond to Soil Type 2. The soil is wetter just above the rods than below them.

Figure 5: Effective saturation around the center rod in Soil Type 1 (solid lines) and Soil Type 2 (dashed lines).

Figure 6 compares the average fluid saturations at the rod boundary with the average within the two blocks of soil. The range of effective saturation estimates at the rod boundary appears as a scatter plot for different time steps. The solid line is the average of the effective saturation at the rod boundary. The dashed line is the average of effective saturation for the soil. Clearly the average effective saturation at the rods increases with time, but the average for the soil does not change. While the effects shown here are more pronounced in Soil Type 1 than Soil Type 2, the results call to question whether the sensors can accurately assess soil moisture if kept in situ for long times.

Figure 6: Average effective saturation at sensor circumference and overall soil block for Soil Type 1 (blue) and Soil Type 2 (green). Also shown is the range of effective saturation at the rod circumference for Soil Type 1 (circles) and Soil Type 2 (squares).

References

1. M.Th. van Genuchten, “A Closed-form Equation for Predicting the Hydraulic of Conductivity of Unsaturated Soils,” Soil Sci. Soc. Am. J., vol. 44, pp. 892–898, 1980.

2. R.H. Brooks and A.T. Corey, “Properties of Porous Media Affecting Fluid Flow,” J. Irrig. Drainage Div., ASCE Proc, vol. 72 (IR2), pp. 61–88, 1966.

Subsurface_Flow_Module/Fluid_Flow/variably_saturated_flow

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

Square 1 (sq1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2.

Locate the section. From the list, choose .

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.1.

Locate the section. In the text field, type -0.5.

Array 1 (arr1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

In the text field, type 3.

Locate the section. In the text field, type 0.5.

Difference 1 (dif1)

In the toolbar, click

and choose .

Select the object only, to add it to the list.

In the window for , locate the section.

Find the subsection. Select the

toggle button.

Select the objects , , and only.

Click

Definitions

Next, create the selection to simplify the evaluation of the results.

Rod

In the toolbar, click

In the window for , type Rod in the text field.

Locate the section. From the list, choose .

Select Boundary 9 only.

Select the check box.

Richards’ Equation (dl)

Richards’ Equation Model 1

First, set up the model for Soil Type 1, which you build with the van Genuchten retention and permeability relationships.

In the window, under click .

In the window for , locate the section.

From the ρ list, choose . In the associated text field, type 1000.

Locate the section. From the list, choose .

In the Ks text field, type 8.25e-5.

Locate the section. From the χf list, choose . In the associated text field, type 4.4e-10.

In the χp text field, type 1e-8.

Locate the section. In the θs text field, type 0.43.

In the θr text field, type 0.045.

Locate the section. In the α text field, type 14.5.

In the n text field, type 2.68.

Initial Values 1

In the window, click .

In the window for , locate the section.

Click the button.

In the Hp text field, type -0.06.

Pressure Head 1

In the toolbar, click

and choose .

Select Boundaries 2 and 3 only.

In the window for , locate the section.

In the Hp0 text field, type -0.06.

Richards’ Equation 2 (dl2)

Richards’ Equation Model 1

Specify equations and boundary conditions for Soil Type 2 using the Brooks and Corey parameterization approach.

In the window, under click .

In the window for , locate the section.

From the ρ list, choose . In the associated text field, type 1000.

Locate the section. From the list, choose .

In the Ks text field, type 5.8333e-5.

Locate the section. From the χf list, choose . In the associated text field, type 4.4e-10.

In the χp text field, type 1e-8.

Locate the section. In the θs text field, type 0.417.

In the θr text field, type 0.02.

Locate the section. From the list, choose .

In the α text field, type 13.8.

In the n text field, type 0.592.

In the l text field, type 1.

Initial Values 1

In the window, click .

In the window for , locate the section.

Click the button.

In the Hp text field, type -0.2.

Pressure Head 1

In the toolbar, click

and choose .

Select Boundaries 2 and 3 only.

In the window for , locate the section.

In the Hp0 text field, type -0.2.

Mesh 1

Free Triangular 1

In the toolbar, click

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Boundaries 5–16 only.

Locate the section. Click the button.

Locate the section. Select the check box.

In the associated text field, type 0.025.

Click

Study 1

In the window, click .

In the window for , locate the section.

Clear the check box.

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the text field, type 0 1 5 10 15.

In the toolbar, click

Results

Follow the steps below to reproduce Figure 3 and Figure 4.

Soil Type 1

In the toolbar, click

and choose .

In the window for , type Soil Type 1 in the text field.

Surface 1

Right-click and choose .

In the window for , locate the section.

In the text field, type dl.Se.

Arrow Surface 1

In the window, right-click and choose .

Contour 1

Right-click and choose .

In the window for , locate the section.

In the text field, type dl.Hp.

Locate the section. From the list, choose .

From the list, choose .

Clear the check box.

In the toolbar, click

Soil Type 2

Right-click and choose .

In the window for , type Soil Type 2 in the text field.

Surface 1

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type dl2.Se.

Arrow Surface 1

In the window, click .

In the window for , locate the section.

In the text field, type dl2.u.

In the text field, type dl2.v.

Contour 1

In the window, click .

In the window for , locate the section.

In the text field, type dl2.Hp.

In the toolbar, click

To generate Figure 5, continue with the steps below.

Effective Saturation

In the toolbar, click

and choose .

In the window for , type Effective Saturation in the text field.

Line Graph 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type dl.Se.

Locate the section. From the list, choose .

In the text field, type atan2(y,x).

From the list, choose .

The function atan2(y,x) along the circle representing the rod produces unphysical values at one point, because there is 180°=-180°. To avoid this, the data are filtered.

Filter 1

Right-click and choose .