Finite water content vadose zone flow method

Derivation

The derivation of the 1-D finite water-content method equation for calculating vertical flux q of water in the vadose zone starts with conservation of mass for an unsaturated porous medium without sources or sinks:

∂ θ ∂ t + ∂ q ∂ z = 0.

Next, the cyclic chain rule and chain rule are used to perform a hodograph transformation and convert the conservation of mass equation into a kinematic equation with z as the dependent variable

∂ z ∂ t = ∂ q ∂ θ .

We next insert the unsaturated Buckingham–Darcy flux:

q = − K ( θ ) ∂ ψ ( θ ) ∂ z + K ( θ ) ,

which yields:

∂ q ∂ θ = ∂ ∂ θ [ K ( θ ) ( 1 − ∂ ψ ( θ ) ∂ z ) ] = ∂ K ( θ ) ∂ θ ( 1 − ∂ ψ ( θ ) ∂ z ) − K ( θ ) ∂ 2 ψ ( θ ) ∂ z ∂ θ .

The term on the far right-hand side of the above equation is negligible in the case of infiltration because we assume the soil has a delta-function diffusivity so that ψ is constant along the infiltration wetting front. Furthermore, evidence suggests that this term is small and negligible along the capillary groundwater wetting front when the groundwater table velocity is less than 0.92 K s . Therefore, the resulting finite water-content vadose zone flow equation is:

( d z d t ) θ = ∂ K ( θ ) ∂ θ [ 1 − ( ∂ ψ ( θ ) ∂ z ) ] .

One way to solve this equation is to solve it for q ( θ , t ) and z ( θ , t ) by integration:

∫ ∂ q ∂ θ d θ = ∫ ∂ z ∂ t d θ

Instead, a finite water-content discretization is used and the integrals are replaced with summations:

∑ j = 1 N [ ∂ q ∂ θ ] j Δ θ = ∑ j = 1 N [ ∂ z ∂ t ] j Δ θ

where N is the total number of finite water content bins.

Using this approach, the conservation equation for each bin is:

[ ∂ q ∂ θ ] j = [ ∂ z ∂ t ] j .

The method of lines is used to replace the partial differential forms on the right-hand side into appropriate finite-difference forms. This process results in a set of three ordinary differential equations that describe the dynamics of infiltration fronts, falling slugs, and groundwater capillary fronts using a finite water-content discretization.

Method essentials

The finite water-content vadose zone flux calculation method replaces the Richards' equation PDE with a set of three ordinary differential equations (ODEs). These three ODEs are developed in the following sections. Furthermore, because the finite water-content method does not explicitly include soil water diffusivity, it necessitates a separate capillary relaxation step. Capillary relaxation represents a free-energy minimization process at the pore scale that produces no advection beyond the REV scale.

Infiltration fronts

With reference to Figure 1, water infiltrating the land surface can flow through the pore space between θ d and θ i . In the context of the method of lines, the partial derivative terms are replaced with:

∂ K ( θ ) ∂ θ = K ( θ d ) − K ( θ i ) θ d − θ i .

Given that any ponded depth of water on the land surface is h p , the Green and Ampt (1911) assumption is employed,

∂ ψ ( θ ) ∂ z = | ψ ( θ d ) | + h p z j ,

represents the capillary head gradient that is driving the flow. Therefore the finite water-content equation in the case of infiltration fronts is:

( d z d t ) j = K ( θ d ) − K ( θ i ) θ d − θ i ( | ψ ( θ d ) | + h p z j + 1 ) .

Falling slugs

After rainfall stops and all surface water infiltrates, water in bins that contains infiltration fronts detaches from the land surface. Assuming that the capillarity at leading and trailing edges of this 'falling slug' of water is balanced, then the water falls through the media at the incremental conductivity associated with the j th   Δ θ bin:

( d z d t ) j = K ( θ j ) − K ( θ j − 1 ) θ j − θ j − 1

Capillary groundwater fronts

In this case, the flux of water to the j th bin occurs between bin j and i. Therefore in the context of the method of lines:

∂ K ( θ ) ∂ θ = K ( θ j ) − K ( θ i ) θ j − θ i ,

and,

∂ ψ ( θ ) ∂ z = | ψ ( θ j ) | H j

which yields:

( d z d t ) j = K ( θ j ) − K ( θ i ) θ j − θ i ( | ψ ( θ j ) | H j − 1 ) .

The performance of this equation was verified, using a column experiment fashioned after that by Childs and Poulovassilis (1962). Results of that validation showed that the finite water-content vadose zone flux calculation method performed comparably to the numerical solution of Richards' equation.

Capillary relaxation

Because the hydraulic conductivity rapidly increases as the water content moves towards saturation, with reference to Fig.1, right-most bins in both capillary groundwater fronts and infiltration fronts can "out-run" their neighbors to the left. In the finite water content discretization, these shocks are dissipated by the process of capillary relaxation, which represents a pore-scale free-energy minimization process that produces no advection beyond the REV scale Numerically, this process is a numerical sort that places the fronts in monotonically-decreasing magnitude from left-right.

Constitutive relations

The finite water content vadose zone flux method works with any monotonic water retention curve/unsaturated hydraulic conductivity relations such as Brooks and Corey Clapp and Hornberger and van Genuchten-Mualem. The method might work with hysteretic water retention relations- these have not yet been tested.

Limitations

The finite water content method lacks the effect of soil water diffusion. This omission does not affect the accuracy of flux calculations using the method because the mean of the diffusive flux is small. Practically, this means that the shape of the wetting front plays no role in driving the infiltration. The method is thus far limited to 1-D in practical applications. The infiltration equation was extended to 2- and quasi-3 dimensions. More work remains in extending the entire method into more than one dimension.

Demonstration video

8-month infiltration simulation from 2600 mm of rain, evapotranspiration in 50 cm root zone, loam soil, and water table fixed at 1m depth. In the video, blue color denotes capillary groundwater, red color denotes infiltration fronts, and pink color denotes falling slugs.

Awards

The paper describing this method was selected by the Early Career Hydrogeologists Network of the International Association of Hydrogeologists to receive the "Coolest paper Published in 2015" award in recognition of the potential impact of the publication on the future of hydrogeology. https://echn.iah.org/coolest-paper-of-2015-award