Soil moisture velocity equation

The soil moisture velocity equation[1] or SMVE is a Lagrangian reinterpretation of the Eulerian Richards' equation wherein the dependent variable is the position z of a wetting front of a particular moisture content ${\displaystyle \theta }$ with time.

${\displaystyle \left.{\frac {dz}{dt}}\right\vert _{\theta }={\frac {\partial K(\theta )}{\partial \theta }}\left[1-\left({\frac {\partial \psi (\theta )}{\partial z}}\right)\right]-D(\theta ){\frac {\partial ^{2}\psi /\partial z^{2}}{\partial \psi /\partial z}}}$

where:

${\displaystyle z}$ is the vertical coordinate [L] (positive downward),

${\displaystyle \theta }$ is the water content of the soil at a point [-]

${\displaystyle K(\theta )}$ is the unsaturated hydraulic conductivity [L T⁻¹],

${\displaystyle \psi (\theta )}$ is the capillary pressure head [L],

${\displaystyle D(\theta )}$ is the soil water diffusivity, which is defined as: ${\displaystyle K(\theta )\partial \psi /\partial \theta }$, [L² T]

${\displaystyle t}$ is time [T].

The first term on the right-hand side of the SMVE is called the "advection-like" term, while the second term is called the "diffusion-like" term. The advection-like term of the Soil Moisture Velocity Equation is particularly useful for calculating the advance of wetting fronts for a liquid invading an unsaturated porous medium under the combined action of gravity and capillarity because it is convertible to an ordinary differential equation by neglecting the diffusion-like term.[5] and it avoids the problem of representative elementary volume by use of a fine water-content discretization and solution method.

This equation was converted into a set of three ordinary differential equations (ODEs)[5] using the method of lines[8] to convert the partial derivatives on the right-hand side of the equation into appropriate finite difference forms. These three ODEs represent the dynamics of infiltrating water, falling slugs, and capillary groundwater, respectively.

This derivation of the 1-D soil moisture velocity equation[1] for calculating vertical flux ${\displaystyle q}$ of water in the vadose zone starts with conservation of mass for an unsaturated porous medium without sources or sinks:

${\displaystyle {\frac {\partial \theta }{\partial t}}+{\frac {\partial q}{\partial z}}=0.}$

We next insert the unsaturated Buckingham–Darcy flux:[9]

${\displaystyle q=-K(\theta ){\frac {\partial \psi (\theta )}{\partial z}}+K(\theta ),}$

yielding Richards' equation[2] in mixed form because it includes both the water content ${\displaystyle \theta }$and capillary head ${\displaystyle \psi (\theta )}$:

${\displaystyle {\frac {\partial \theta }{\partial t}}={\frac {\partial }{\partial z}}\left[K(\theta )\left({\frac {\partial \psi (\theta )}{\partial z}}-1\right)\right]}$.

Applying the chain rule of differentiation to the right-hand side of Richards' equation:

${\displaystyle {\frac {\partial \theta }{\partial t}}={\frac {\partial }{\partial z}}K(\theta (z,t)){\frac {\partial }{\partial z}}\psi (\theta (z,t))+K(\theta ){\frac {\partial ^{2}}{\partial z^{2}}}\psi (\theta (z,t))-{\frac {\partial }{\partial z}}K(\theta (z,t))}$.

Assuming that the constitutive relations for unsaturated hydraulic conductivity and soil capillarity are solely functions of the water content, ${\displaystyle K=K(\theta )}$and ${\displaystyle \psi =\psi (\theta )}$, respectively:

${\displaystyle {\frac {\partial \theta }{\partial t}}=K'(\theta )\psi '(\theta )\left({\frac {\partial \theta }{\partial z}}\right)^{2}+K(\theta )\left[\psi ''(\theta )\left({\frac {\partial \theta }{\partial z}}\right)^{2}+\psi '(\theta ){\frac {\partial ^{2}\theta }{\partial z^{2}}}\right]-K'(\theta ){\frac {\partial \theta }{\partial z}}}$.

This equation implicitly defines a function ${\displaystyle Z_{R}(\theta ,t)}$that describes the position of a particular moisture content within the soil using a finite moisture-content discretization. Employing the Implicit function theorem, which by the cyclic rule required dividing both sides of this equation by ${\displaystyle {-\partial \theta }/{\partial z}}$ to perform the change in variable, resulting in:

${\displaystyle {\frac {\partial Z_{R}}{\partial t}}=-K'(\theta )\psi '(\theta ){\frac {\partial \theta }{\partial z}}-K(\theta )\psi ''(\theta ){\frac {\partial \theta }{\partial z}}-K(\theta )\psi '(\theta ){\frac {\partial ^{2}\theta /\partial z^{2}}{\partial \theta /\partial z}}+K'(\theta )}$,

which can be written as:

${\displaystyle {\frac {\partial Z_{R}}{\partial t}}=-K'(\theta )\left[{\frac {\partial \psi (\theta )}{\partial z}}-1\right]-K(\theta )\left[\psi ''(\theta ){\frac {\partial \theta }{\partial z}}+\psi '(\theta ){\frac {\partial ^{2}\theta /\partial z^{2}}{\partial \theta /\partial z}}\right]}$.

Inserting the definition of the soil water diffusivity:

${\displaystyle D(\theta )\equiv K(\theta ){\frac {\partial \psi }{\partial \theta }}}$

into the previous equation produces:

${\displaystyle {\frac {\partial Z_{R}}{\partial t}}=-K'(\theta )\left[{\frac {\partial \psi (\theta )}{\partial z}}-1\right]-D(\theta ){\frac {\partial ^{2}\psi /\partial z^{2}}{\partial \psi /\partial z}}}$

If we consider the velocity of a particular water content ${\displaystyle \theta }$, then we can write the equation in the form of the Soil Moisture Velocity Equation:

${\displaystyle \left.{\frac {dz}{dt}}\right\vert _{\theta }={\frac {\partial K(\theta )}{\partial \theta }}\left[1-\left({\frac {\partial \psi (\theta )}{\partial z}}\right)\right]-D(\theta ){\frac {\partial ^{2}\psi /\partial z^{2}}{\partial \psi /\partial z}}}$

Written in moisture content form, 1-D Richards' equation is[10]

${\displaystyle {\frac {\partial \theta }{\partial t}}={\frac {\partial }{\partial z}}\left(D(\theta ){\frac {\partial \theta }{\partial z}}\right)+{\frac {\partial K(\theta )}{\partial z}}}$

Where D(θ) [L²/T] is 'the soil water diffusivity' as previously defined.

Note that with ${\displaystyle \theta }$ as the dependent variable, physical interpretation is difficult because all the factors that affect the divergence of the flux are wrapped up in the soil moisture diffusivity term ${\displaystyle D(\theta )}$. However, in the SMVE, the three factors that drive flow are in separate terms that have physical significance.

The primary assumptions used in the derivation of the Soil Moisture Velocity Equation are that ${\displaystyle K=K(\theta )}$ and ${\displaystyle \psi =\psi (\theta )}$ are not overly restrictive. Analytical and experimental results show that these assumptions are acceptable under most conditions in natural soils. In this case, the Soil Moisture Velocity Equation is equivalent to the 1-D Richards' equation, albeit with a change in dependent variable. This change of dependent variable is convenient because it reduces the complexity of the problem because compared to Richards' equation, which requires the calculation of the divergence of the flux, the SMVE represents a flux calculation, not a divergence calculation. The first term on the right-hand side of the SMVE represents the two scalar drivers of flow, gravity and the integrated capillarity of the wetting front. Considering just that term, the SMVE becomes:

${\displaystyle {\frac {\partial Z_{R}}{\partial t}}=-K'(\theta )\left[{\frac {\partial \psi (\theta )}{\partial z}}-1\right]}$

where ${\displaystyle {\partial \psi (\theta )}/{\partial z}}$ is the capillary head gradient that is driving the flux and the remaining conductivity term ${\displaystyle K'(\theta )}$ represents the ability of gravity to conduct flux through the soil. This term is responsible for the true advection of water through the soil under the combined influences of gravity and capillarity. As such, it is called the "advection-like" term.

Neglecting gravity and the scalar wetting front capillarity, we can consider only the second term on the right-hand side of the SMVE. In this case the Soil Moisture Velocity Equation becomes:

${\displaystyle {\frac {\partial Z_{R}}{\partial t}}=-D(\theta ){\frac {\partial ^{2}\psi /\partial z^{2}}{\partial \psi /\partial z}}}$

This term is strikingly similar to Fick's second law of diffusion. For this reason, this term is called the "diffusion-like" term of the SMVE.

This term represents the flux due to the shape of the wetting front ${\displaystyle -D(\theta ){\partial ^{2}\psi /\partial z^{2}}}$, divided by the spatial gradient of the capillary head ${\displaystyle {\partial \psi /\partial z}}$. Looking at this diffusion-like term, it is reasonable to ask when might this term be negligible? The first answer is that this term will be zero when the first derivative ${\displaystyle <\partial \psi /\partial z=C}$, because the second derivative will equal zero. One example where this occurs is in the case of an equilibrium hydrostatic moisture profile, when ${\displaystyle \partial \psi /\partial z=-1}$ with z defined as positive upward. This is a physically realistic result because an equilibrium hydrostatic moisture profile is known to not produce fluxes.

Another instance when the diffusion-like term will be nearly zero is in the case of sharp wetting fronts, where the denominator of the diffusion-like term ${\displaystyle \partial \psi /\partial z\to \infty }$, causing the term to vanish. Notably, sharp wetting fronts are notoriously difficult to resolve and accurately solve with traditional numerical Richards' equation solvers.[11]

Finally, in the case of dry soils, ${\displaystyle K(\theta )}$ tends towards ${\displaystyle 0}$, making the soil water diffusivity ${\displaystyle D(\theta )}$ tend towards zero as well. In this case, the diffusion-like term would produce no flux.

Comparing against exact solutions of Richards' equation for infiltration into idealized soils developed by Ross & Parlange (1994)[12] revealed[1] that indeed, neglecting the diffusion-like term resulted in accuracy >99% in calculated cumulative infiltration. This result indicates that the advection-like term of the SMVE, converted into an ordinary differential equation using the method of lines, is an accurate ODE solution of the infiltration problem. This is consistent with the result published by Ogden et al.[5] who found errors in simulated cumulative infiltration of 0.3% using 263 cm of tropical rainfall over an 8-month simulation to drive infiltration simulations that compared the advection-like SMVE solution against the numerical solution of Richards' equation.

The advection-like term of the SMVE can be solved using the method of lines and a finite moisture content discretization. This solution of the SMVE advection-like term replaces the 1-D Richards' equation PDE with a set of three ordinary differential equations (ODEs). These three ODEs are:

Infiltration fronts

Infiltration fronts in finite water-content domain

With reference to Figure 1, water infiltrating the land surface can flow through the pore space between ${\displaystyle \theta _{d}}$ and ${\displaystyle \theta _{i}}$. Using the method of lines to convert the SMVE advection-like term into an ODE:

${\displaystyle {\frac {\partial K(\theta )}{\partial \theta }}={\frac {K(\theta _{d})-K(\theta _{i})}{\theta _{d}-\theta _{i}}}.}$

Given that any ponded depth of water on the land surface is ${\displaystyle h_{p}}$, the Green and Ampt (1911)[13] assumption is employed,

${\displaystyle {\frac {\partial \psi (\theta )}{\partial z}}={\frac {|\psi (\theta _{d})|+h_{p}}{z_{j}}},}$

represents the capillary head gradient that is driving the flow in the ${\displaystyle j^{th}}$ discretization or "bin". Therefore, the finite water-content equation in the case of infiltration fronts is:

${\displaystyle \left({\frac {dz}{dt}}\right)_{j}={\frac {K(\theta _{d})-K(\theta _{i})}{\theta _{d}-\theta _{i}}}\left({\frac {|\psi (\theta _{d})|+h_{p}}{z_{j}}}+1\right).}$

Falling slugs

Falling slugs in the finite water-content domain. The water in each bin is considered a separate slug.

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 ${\displaystyle j^{\text{th}}\ \Delta \theta }$ bin:

${\displaystyle \left({\frac {dz}{dt}}\right)_{j}={\frac {K(\theta _{j})-K(\theta _{j-1})}{\theta _{j}-\theta _{j-1}}}}$.

This approach to solving the capillary-free solution is very similar to the kinematic wave approximation.

Capillary groundwater fronts

Groundwater capillary fronts in finite water-content domain

In this case, the flux of water to the ${\displaystyle j^{\text{th}}}$ bin occurs between bin j and i. Therefore, in the context of the method of lines:

${\displaystyle {\frac {\partial K(\theta )}{\partial \theta }}={\frac {K(\theta _{j})-K(\theta _{i})}{\theta _{j}-\theta _{i}}},}$

and

${\displaystyle {\frac {\partial \psi (\theta )}{\partial z}}={\frac {|\psi (\theta _{j})|}{H_{j}}}}$

which yields:

${\displaystyle \left({\frac {dH}{dt}}\right)_{j}={\frac {K(\theta _{j})-K(\theta _{i})}{\theta _{j}-\theta _{i}}}\left({\frac {|\psi (\theta _{j})|}{H_{j}}}-1\right).}$

Note the "-1" in parentheses, representing the fact that gravity and capillarity are acting in opposite directions. The performance of this equation was verified,[7] using a column experiment fashioned after that by Childs and Poulovassilis (1962).[6] Results of that validation showed that the finite water-content vadose zone flux calculation method performed comparably to the numerical solution of Richards' equation. The photo shows apparatus. Data from this column experiment are available by clicking on this hot-linked DOI. These data are useful for evaluating models of near-surface water table dynamics.

It is noteworthy that the SMVE advection-like term solved using the finite moisture-content method completely avoids the need to estimate the specific yield. Calculating the specific yield as the water table nears the land surface is made cumbersome my non-linearities. However, the SMVE solved using a finite moisture-content discretization essentially does this automatically in the case of a dynamic near-surface water table.

Column experiment used to observe moisture response in a fine sand above a moving water table. Note stepper-motor controlled constant head reservoir (white bucket).

The paper on the Soil Moisture Velocity Equation was highlighted by the editor in the issue of J. Adv. Modeling of Earth Systems when the paper was first published, and is in the public domain. The paper may be freely downloaded here by anyone. The paper describing the finite moisture-content solution of the advection-like term of the Soil Moisture Velocity Equation was selected to receive the 2015 Coolest Paper Award by the early career members of the International Association of Hydrogeologists.

• YouTube video of SMVE-based solution slowed during rainfall to highlight behavior, with fixed water table at 1.0 m and evapotranspiration from a 0.5 m root zone