COMSOL 6.4 - Variably Saturated Flow and Transport

Variably Saturated Flow and Transport — Sorbing Solute

In this example, water ponded in a ring on the ground moves into a relatively dry soil column and carries a chemical with it. As it moves through the variably saturated soil column, the chemical attaches to soil particles, slowing the solute transport relative to the water. Additionally the chemical concentrations decay from biodegradation in both the liquid and the solid phase. The inspiration for the problem comes from Ref. 1.

Figure 1: Geometry of the infiltration ring and soil column.

The example uses the Richards’ Equation interface to define nonlinear relationships with retention and permeability properties according to van Genuchten (Ref. 2). In the Richards’ Equation interface you can also define these material properties with Brooks and Corey’s analytic permeability and retention formulas (Ref. 3) or by interpolation from experimental data.

Richards’ equation accounts for changes in the fluid volume fraction with time, and also for changes in the storage related to variations in the pressure head according to Bear (Ref. 4). With the storage terms, the transport of diluted species in porous media equations in the Subsurface Flow Module also account explicitly for the time changes in liquid and air content.

Model Definition

The water moves from a ring on the ground into the subsurface. The 0.25-m radius ring ponds the water to a depth of 0.01 m but is open to the ground surface. Permeable soils exist to a depth of 1.3 m. The soil in the uppermost layer is slightly less permeable than the bottom one. The lower layer sits above relatively impermeable soil, so only a very small amount of leakage exits from the base. The flow is symmetric about the line r = 0. The initial distribution of pressure heads is known.

The water in the ring contains a dissolved solute at a constant concentration, c0. The solute enters the ground with the water and moves through the subsurface by advection and dispersion. Additionally, the solute adsorbs or attaches to solid surfaces, which reduces the aqueous concentrations and also slows solute movement relative to the water. Microbial degradation also reduces both the liquid-phase and solid-phase concentrations. The sorption and the biodegradation are linearly proportional to aqueous concentrations. The fluid in the ring is the only chemical source, and the solute is free to leave the soil column with the fluid flux. Initially the soil is free of the solute. You track its transport for five days.

Fluid Flow

Richards’ equation governs the saturated-unsaturated flow of water in the soil. The soil air is open to the atmosphere, so you can assume that pressure changes in air do not affect the flow and use Richards’ equation here for single-phase flow. Given by Ref. 4, Richards’ equation in pressure head reads

where Cm denotes specific moisture capacity (m−1); Se is the effective saturation of the soil (dimensionless); Sp is a storage coefficient (m−1); Hp is the pressure head (m), which is proportional to the dependent variable, p (Pa); t is time; K equals the hydraulic conductivity (m/s); D is the direction (typically, the z direction) that represents vertical elevation (m).

To be able to combine boundary conditions and sources with the Darcy’s Law formulation, COMSOL Multiphysics converts Richards’ equation to SI units and solves for the pressure (SI unit: Pa). Hydraulic head, H, pressure head, Hp, and elevation D are related to pressure p as

Also, the permeability κ (SI unit: m2) and hydraulic conductivity K (SI unit: m/s) are related to the viscosity μ (SI unit: Pa·s) and density ρ (SI unit: kg/m3) of the fluid, and the acceleration of gravity g (SI unit: m/s2) by

In this problem, Sp = (θs − θr)/(1 m·ρg) where θs and θr denote the volume fraction of fluid at saturation and after drainage, respectively. For more details see The Richards’ Equation Interface in the Subsurface Flow Module User’s Guide.

With variably saturated flow, fluid moves through but may or may not completely fill the pores in the soil, and θ denotes the volume fraction of fluid within the soil. The coefficients Cm, Se, and K vary with the pressure head, Hp, and with θ, making Richards’ equation nonlinear. The specific moisture capacity, Cm, relates variations in soil moisture to pressure head as Cm = ∂θ/∂Hp. In the governing equation, Cm defines storage changes produced by varying fluid content because Cm∂Hp/∂t = ∂θ/∂t. Because Cm goes to zero at saturation, time change in storage relates to compression of the aquifer and water under saturated conditions. The saturated storage comes about with the effective saturation, as represented by the second term in the time-coefficient. Furthermore, K is a function that defines how readily the porous media transmits fluid. The relative permeability of the soil, κr, increases with fluid content giving K = Ksκr, where Ks (m/d) is the constant hydraulic conductivity at saturation.

This example uses predefined interfaces for van Genuchten formulas (Ref. 2) to represent how the four retention and permeability properties — θ, Cm, Se, and κr = K/Ks — vary with the solution Hp. For more details about these relationships see the section Retention and Permeability Relationships in the Theory for the Richards’ Equation Interface chapter in the Subsurface Flow Module User’s Guide.

Boundary Conditions and Initial Conditions

The problem statement records all the boundary conditions you need for this model. The level of water in the ring is known at 0.01 m, giving a Dirichlet constraint on pressure head. Approximate the small leak from the base, N0, as 0.01Ks. Infinite Elements are used that apply a coordinate scaling in the r direction. The boundary condition applied to these infinite elements is effectively applied at a very large distance. Together with no flow crossing the surface outside of the pressure ring or the infinite elements, the following expressions summarize the boundary conditions:

In these expressions, n is the unit vector normal to the bounding surface, and u is Darcy’s velocity field. The initial conditions specify the pressure head in the modeling domain as

Transport of Diluted Species in Porous Media

Groundwater flow and solute transport are linked by fluid velocities. With the form of the transport equation that follows, the fluid velocities need to come from Darcy’s law:

In this expression, u is Darcy’s velocity field (SI unit: m/s).

The equation that governs advection, dispersion, sorption, and decay of solutes in groundwater is

It describes time rate of change in two terms: c denotes dissolved concentration (mol/m3), and cP is the mass of adsorbed contaminant per dry unit weight of solid (mg/kg). Further, θ denotes the volume fluid fraction (porosity), and ρb is the bulk density (kg/m3). Because ρb amounts to the dried solid mass per bulk volume of the solids and pores together, the term ρb cP gives solute mass attached to the soil as a concentration. In the equation, DL is the hydrodynamic dispersion tensor (m2/d); RL represents reactions in water (mol/(m3·d)); and RP equals reactions involving solutes attached to soil particles (mol/(m3·d)). Finally, Sc is solute added per unit volume of soil per unit time (mol/(m3·d)).

It is far more convenient to solve the above equation only for dissolved concentration. This requires expanding the left-hand side to

and inserting a few definitions.

In this problem, solute mass per solid mass, cP, relates to dissolved concentration, c, through a linear isotherm or partition coefficient kP (m3/kg) where cP = kP c. Because the relationship is linear, the derivative is kP = ∂cP/∂c. Making those substitutions gives the form of the solute transport problem you solve:

In the equation, ϕL and ϕP denote the decay rates (d−1) for the dissolved and adsorbed solute concentrations, respectively.

Select the option for in the settings for the Unsaturated Porous Media feature to employ results from the flow equation in the solute-transport model:

Note that COMSOL Multiphysics solves for pressure, p, and converts to Hp based on the fluid weight.

The hydrodynamic dispersion tensor, DL, describes mechanical spreading from groundwater movement in addition to chemical diffusion:

where DLii are the diagonal entries in the dispersion tensor; DLij are the cross terms; α is the dispersivity (m); the subscripts “1” and “2” denote longitudinal and transverse dispersivities, respectively; Dm denotes the coefficient of molecular diffusion (m2/d); and τL is a tortuosity factor that reduces impacts of molecular diffusion for porous media relative to free water. Here τL = θ −7/3 θs2.

Boundary Conditions and Initial Conditions

The boundary and initial conditions in the sorbing-solute problem are straightforward. The solute enters only with the water from the ring at a concentration c0. The solute is free to leave, but there is only minimal leakage from the lower boundary and no flow from the sides. Transport is symmetric about the line r = 0. The boundary conditions in this problem are:

where n is the unit vector normal to the boundary. Because the soil is pristine at the start of the experiment, the initial condition is one of zero concentration.

Model Data

The following table provides data for the fluid-flow model:


Variable	Unit	Description	upper layer	Lower layer
Ks	m/d	Saturated hydraulic conductivity	0.298	0.454
θs		Porosity/void fraction	0.399	0.339
θr		Residual saturation	0.001	0.001
α	m-1	alpha parameter	1.74	1.39
n		n parameter	1.38	1.60
m		m parameter	1-1/n	1-1/n
l		Pore connectivity parameter	n/a
Hp0	m	Pressure head in ring	0.01
Hp,init	m	Initial pressure head	-(z+1.2) -0.2(z+0.4)	-(z+1.2)

The inputs needed for the solute-transport model are:


Variable	Units	Description	Value
ρb	kg/m3	Bulk density	1400
kp	m3/kg	Partition coefficient	0.0001
Dm	m2/d	Coefficient of molecular diffusion	0.00374
αr	m	Longitudinal dispersivity	0.005
αz	m	Transverse dispersivity	0.001
	d-1	Decay rate in liquid	0.05
	d-1	Decay rate on solid	0.01
c0	mol/m3	Solute concentration in ring	1.0

Results and Discussion

Figure 2 gives the solution to the fluid-flow problem at 0.3 days, 1 day, and 5 days, respectively. The images show effective saturation (surface plot) and velocities (arrows). They illustrate the soil wetting with time. As the arrows indicate, the velocities just below the ring are high relative to the remainder of the soil column.

Figure 2: Estimates of effective saturation (surface plot) and velocity (arrows) in variably saturated soil after 0.3, 1, and 5 days.

Figure 3 gives the concentrations for 0.3 days, 1 day, and 5 days, respectively. It illustrates how the solute concentrations (surface plot) enter and move through the soil. Because the retardation factor depends on soil moisture, its value varies with the solution (see Figure 4).

Figure 3: Solution for dissolved concentrations at 0.3, 1, and 5 days.

Figure 4 shows an image of the retardation factor and the pressure head at the end of the simulation time interval. For variably saturated solute transport, the retardation factor changes with time. As shown in this image, the process of sorption has the greatest potential to slow the contaminant where the soils are relatively dry. The retardation coefficient here ranges from roughly 1.35 to 1.65, and the solute moves at approximately the velocity of the groundwater.

Figure 4: Retardation factor (surface) and pressure head (contours) after 5 days.

Notes About the COMSOL Implementation

This model makes use of the feature. It performs a coordinate scaling to the selected domain such that boundary conditions on the outside of the infinite element layer are effectively applied at a very large distance. Therefore unwanted effects of artificial boundary conditions on the region of interest are suppressed. This allows to model details in an area which is actually very large or infinite.

References

1. J. Simunek, T. Vogel, and M.Th. van Genuchten, “The SWMS_2D code for simulating water flow and solute transport in two-dimensional variably saturated media,” ver. 1.1., Research Report No. 132, U.S. Salinity Laboratory, USDA, 1994.

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

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

4. J. Bear, Hydraulics of Groundwater, McGraw Hill, 1978.

Subsurface_Flow_Module/Solute_Transport/sorbing_solute

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

Global Definitions

Load the parameters from file.

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 sorbing_solute_parameters.txt.

Geometry 1

The modeling domain is made up of the two permeable soil layers, each of which is represented by a rectangular domain in 2D axisymmetry.

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 1.5.

In the text field, type 0.9.

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

Click to expand the section. Select the checkbox.

Clear the checkbox.

In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	0.25

This additional layer to the right will later be used to define an . Read more about it in the Notes About the COMSOL Implementation section. Proceed with the second soil layer.

Rectangle 2 (r2)

Right-click and choose .

In the window for , locate the section.

In the text field, type 0.4.

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

To finish the model geometry, add a point on the top boundary marking the pond’s outer rim.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.25.

Click

Now, define the .

Definitions

Infinite Element Domain 1 (ie1)

In the toolbar, click

Select Domains 3 and 4 only.

In the window for , locate the section.

From the list, choose .

Richards’ Equation (dl)

Begin by specifying the properties for the bottom soil layer in the default node, then duplicate this node and modify the domain selection and properties to match the top layer.

Unsaturated Porous Medium 1

In the window, under > click .

In the window for , locate the section.

From the list, choose . In the Sp text field, type Sp_1.

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the Ks text field, type Ks_1.

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

In the n text field, type n_1.

In the θr text field, type thetar_1.

Unsaturated Porous Medium 2

In the window, under > right-click and choose .

Select Domains 2 and 4 only.

In the window for , locate the section.

In the Sp text field, type Sp_2.

Porous Matrix 1

In the window, expand the node, then click .

In the window for , locate the section.

In the Ks text field, type Ks_2.

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

In the n text field, type n_2.

In the θr text field, type thetar_2.

Gravity 1

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Initial Values 1

In the window, click .

In the window for , locate the section.

Click the button.

In the Hp text field, type -(z+1.2).

Initial Values 2

Right-click > > and choose .

Select Domains 2 and 4 only.

In the window for , locate the section.

In the Hp text field, type -(z+1.2)-0.2*(z+0.4).

Pressure Head 1

In the toolbar, click

and choose .

Select Boundary 5 only.

In the window for , locate the section.

In the Hp0 text field, type Hp0.

Pervious Layer 1

In the toolbar, click

and choose .

Select Boundaries 2 and 8 only.

In the window for , locate the section.

In the Hb text field, type -2.

In the Rb text field, type 1/5[d].

Transport of Diluted Species in Porous Media (tds)

Now, set up the transport equation for an unsaturated porous medium, accounting for dispersion and adsorption.

Unsaturated Porous Medium 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Liquid 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the θl text field, type dl.theta_l.

This corresponds to the liquid volume fraction.

From the list, choose .

From the dHp/dt list, choose .

In the Cm text field, type dl.Cm.

Locate the section. From the u list, choose .

Locate the section. In the DL,c text field, type Dl.

Add adsorption as a subnode to the node.

Unsaturated Porous Medium 1

In the window, click .

Adsorption 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select the checkbox.

In the cP,c text field, type kp*c.

Unsaturated Porous Medium 1

In the window, click .

Dispersion 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Specify the α matrix as


alphar
alphaz

Reactions 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Rc text field, type (phip*kp*rhob-phil*dl.theta_l)*c.

Outflow 1

In the toolbar, click

and choose .

Select Boundaries 2, 8, 12, and 13 only.

Concentration 1

In the toolbar, click

and choose .

Select Boundary 5 only.

In the window for , locate the section.

Select the checkbox.

In the c0,c text field, type c0.

Materials

Some required material properties have not yet been defined. This is indicated by a small red cross at the material node. Continue as follows to add the missing properties.

Porous Material: Lower Layer

In the window, under right-click and choose > .

In the window for , type Porous Material: Lower Layer in the text field.

Locate the section. Click

Fluid 1 (pmat1.fluid1)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rhof	kg/m³	Basic

Solid 1 (pmat1.solid1)

In the window, click .

In the window for , locate the section.

In the θs text field, type 1-poro_1.

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


Property	Variable	Value	Unit	Property group
Density	rho	rhob	kg/m³	Basic

Porous Material: Upper Layer

In the window, under > right-click and choose .

In the window for , type Porous Material: Upper Layer in the text field.

Select Domains 2 and 4 only.

Solid 1 (pmat2.solid1)

In the window, expand the > > node, then click .

In the window for , locate the section.

In the θs text field, type 1-poro_2.

Mesh 1

Using a mapped mesh is a good idea for this geometry. It uses less mesh elements while keeping the accuracy compared to using a triangular mesh with the same mesh size.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Mapped 1

In the toolbar, click

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domain 2 only.

Locate the section. Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 0.02.

In the window, right-click and choose .

Study 1

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 range(0,0.1,0.9) range(1,1,5).

In the window, click .

In the window for , locate the section.

Clear the checkbox, to disable the automatic generation of default plots. Instead, choose the relevant plots from the Result Templates after computation.

In the toolbar, click