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 C denotes specific moisture capacity (m−1); Se is the effective saturation of the soil (dimensionless); S 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, S = (θ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 C, Se, and K vary with the pressure head, Hp, and with θ, making Richards’ equation nonlinear. The specific moisture capacity, C, relates variations in soil moisture to pressure head as C = ∂θ/∂Hp. In the governing equation, C defines storage changes produced by varying fluid content because C∂Hp/∂t = ∂θ/∂t. Because C, the first term in the time coefficient, 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, kr, increases with fluid content giving K = Kskr, 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 — θ, C, Se, and kr = K/Ks — vary with the solution Hp. The van Genuchten expressions read as follows:

Here α, n, m, and l are constants that specify the type of medium, with m = 1 − 1/n. In the equations, the system reaches saturation when fluid pressure is atmospheric (that is, Hp = 0). When the soil fully saturates, the four parameters reach constant values.

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 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,

and

denote the decay rates (d−1) for the dissolved and adsorbed solute concentrations, respectively.

Select the option for in the settings for the Partially Saturated 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 and Figure 3 give the solution to the fluid-flow problem at 0.3 days and 1 day, respectively. The images show effective saturation (surface plot), pressure head (contours), and velocities (arrows). The figures 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), pressure head (contours), and velocity (arrows) in variably saturated soil after 0.3 days.

Figure 3: Estimates of effective saturation (surface plot), pressure head (contours), and velocity (arrows) in variably saturated soil after 1 day.

Figure 4 and Figure 5 give the concentrations for 0.3 days and 1 day, respectively, along with the retardation factor. They illustrate 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.

Figure 4: Solution for dissolved concentrations (surface plot) and retardation factor (contours) at 0.3 days.

Figure 5: Solution for dissolved concentrations (surface plot) and retardation factor (contours) at 1 day.

Figure 6 shows an image of the retardation factor 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 6: Snapshot of the retardation factor (surface and contours) at 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. 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