COMSOL 6.4 - 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.

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. For the detailed formulas, see the section Retention and Permeability Relationships in the Subsurface Flow Module User’s Guide.

To find a unique solution to this problem, you must specify initial and boundary conditions. Initially, the column has uniform pressure head of Hp0. No flow conditions are applied to the rings and to the sides. At the surface and bottom of the geometry the pressure head Hp0is specified.

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).

The results can also be displayed in 3D using an extrusion dataset. As an example, the effective saturation is plotted in slices within a 3D domain in Figure 7 and Figure 8.

Figure 7: Slice plot of the effective saturation and streamlines of the velocity at 15 minutes for Soil Type 1 (van Genuchten), displayed within a 3D domain.

Figure 8: Slice plot of the effective saturation and streamlines of the velocity at 15 minutes for Soil Type 2(Brooks and Corey), displayed within a 3D domain.

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 .

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.

Click to select the

toggle button for .

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 checkbox.

To investigate the two different retention models, it is useful to work with parameter cases. First, load a set of parameters and create two different parameter sets.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file variably_saturated_flow_parameters.txt.

In the toolbar, click

In the window for , type Case 1: van Genuchten in the text field.

In the toolbar, click

In the window for , type Case 2: Brooks and Corey in the text field.

Locate the section. In the table, enter the following settings:


Name	Expression	Description
poro	0.417	Porosity
theta_r	0.02	Residual liquid volume fraction
K	5.8333e-5[m/s]	Hydraulic conductivity
alpha	13.8[1/m]	Constitutive relation constant
n	0.592	Constitutive relation constant
l	1	Constitutive relation constant
Hp0	-0.2[m]	Pressure head

Continue with setting up the physics using the parameters.

Richards’ Equation (dl)

Fluid 1

In the window, under > > click .

In the window for , locate the section.

From the list, choose .

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

From the χf list, choose . In the associated text field, type chi_f.

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the εp list, choose . In the associated text field, type poro.

In the χp text field, type chi_p.

From the list, choose .

In the Ks text field, type K.

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

In the l text field, type l.

In the n text field, type n.

In the θr text field, type theta_r.

Unsaturated Porous Medium (van Genuchten)

In the window, under > click .

In the window for , type Unsaturated Porous Medium (van Genuchten) in the text field.

Unsaturated Porous Medium (Brooks and Corey)

Right-click and choose .

In the window for , type Unsaturated Porous Medium (Brooks and Corey) in the text field.

Porous Matrix 1

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Initial Values 1

In the window, under > click .

In the window for , locate the section.

Click the button.

In the Hp text field, type Hp0.

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 Hp0.

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 checkbox. In the associated text field, type 0.025.

Click

Study 1

In the window, click .

In the window for , locate the section.

Clear the checkbox.

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.

Disable the feature that uses the Brooks and Corey retention model for this study.

Locate the section. Select the checkbox.

In the tree, select > > .

Click

Add a to run this study for the set of the van Genuchten parameters.

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Click

In the table, enter the following settings:


Switch	Cases	Case numbers
Parameters 1	User defined	1

In the window, click .

In the window for , type Study 1: van Genuchten in the text field.

In the toolbar, click

Results

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

Soil Type 1

In the toolbar, click

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 checkbox.

In the toolbar, click

Soil Type 1

Set up a new study for the Brooks and Corey parameter set.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Study 2

Step 1: Time Dependent

In the window for , locate the section.

From the list, choose .

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

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Click

In the table, enter the following settings:


Switch	Cases	Case numbers
Parameters 1	User defined	2

In the window, click .

In the window for , type Study 2: Brooks and Corey in the text field.

Locate the section. Clear the checkbox.

In the toolbar, click

Results

Soil Type 2

In the window, right-click and choose .

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

Locate the section. From the list, choose .

To generate Figure 5, continue with the steps below.

Effective Saturation

In the toolbar, click

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 .

In the window for , locate the section.

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

In the toolbar, click