goVSF: Modelling multi-phase flow in variably saturated porous media

a fully-implicit solution to the one-dimensional, unsteady Richards equation of variably saturated flow in porous media with vapour exchange, written in Go.

Following the solution of Bittelli, M., Campbell, G.S., and Tomei, F., (2015). Soil Physics with Python. Oxford University Press.

dependencies

goHydro porousmedia - a Go struct used to define the properties of a variable porous media profile

Theory: multiphase flow through porous media

Liquid phase

Flux through a porous medium is defined as:

$$ f_l=-K\frac{d\psi}{dz} $$

where $f_l$ is liquid flux density [kg m $^{-2}$ s $^{-1}$], $\psi$ is the matric potential [J kg $^{-1}$] $\left(\text{note: }\psi=gh\right)$, and $K$ is the hydraulic conductivity [kg s m $^{-3}$]. Note: $K_{kg\cdot s\cdot m^{-3}}=\frac{\rho_l}{g}K_{m\cdot s^{-1}}$

When combined with the mass conservation equation, the result yields the transient vertical (1D) Richards (1931) equation:

$$ \rho_l\frac{\partial\theta}{\partial t}=\frac{\partial}{\partial z}\left[K(\psi)\left(\frac{\partial\psi}{\partial z}+g\right)\right], $$

which is often expressed in its so-called $\psi$--based form:

$$ \rho_lC(\psi)\frac{\partial\psi}{\partial t}=\frac{\partial}{\partial z}\left[K(\psi)\left(\frac{\partial\psi}{\partial z}+g\right)\right], $$

where $\rho_l$ is the density of the liquid [kg m $^{-3}$], $g$ is the acceleration to to gravity ($\approx9.80665$ m s $^{-2}$ or J kg $^{-1}$ m $^{-1}$). The capacity [kg J] term is solved using the Campbell (1974) model:

$$ C(\psi)=\frac{d\theta}{d\psi}=\frac{-\theta}{b\psi} $$

and conductivity:

$$ K(\psi)= \begin{cases} K_s\left(\frac{\psi_e}{\psi}\right)^{2+3/b} \qquad &\text{if $\psi_m\leq\psi_e$} \\ K_s &\text{otherwise} \end{cases} $$

where $K_s$ is the saturated conductivity [kg s m $^{-3}$], $b$ is a shape parameter and $\psi_e$ is the air-entry potential [J kg $^{-1}$]. Furthermore, water content ($\theta$) [m $^3$ m $^{-3}$] and $\psi$ are assumed related by the power-law relationship:

$$ \theta(\psi)= \begin{cases} \theta_s\left(\frac{\psi}{\psi_e}\right)^{-1/b} \qquad &\text{if $\psi_m\leq\psi_e$} \\ \theta_s &\text{otherwise} \end{cases} $$

or conversely,

$$ \psi(\theta)= \begin{cases} \psi_e\left(\frac{\theta}{\theta_s}\right)^{-b} \qquad &\text{if $\theta\leq\theta_s$} \\ \psi_e &\text{otherwise} \end{cases} $$

where $\theta_s$ is the water content at saturation. The water balance for cell $i$ can be solved using the its finite-difference form (linear method).

$$ \frac{\rho_l\left(\theta_i^{j+1}-\theta_i^j\right)\left(z_{i+1}-z_{i-1}\right)}{2\Delta t}=\frac{\overline{K}i\left(\psi{i+1}-\psi_i\right)}{z_{i+1}-z_i}-\frac{\overline{K}{i-1}\left(\psi_i-\psi{i-1}\right)}{z_i-z_{i-1}}+u_i $$

or in its $\psi$-based form:

$$ \frac{\rho_lC_i\left(\psi_i^{j+1}-\psi_i^j\right)\left(z_{i+1}-z_{i-1}\right)}{2\Delta t}=\frac{\overline{K}i\left(\psi{i+1}-\psi_i\right)}{z_{i+1}-z_i}-\frac{\overline{K}{i-1}\left(\psi_i-\psi{i-1}\right)}{z_i-z_{i-1}}+u_i $$

The source term $u_i$ could include extraction by roots, evaporation, etc. In most cases, for vertical simulations, the source term represents gravitational flux, where:

$$ u_i=g\left(\overline{K}_{i-1}-\overline{K}_i\right). $$

Vapour phase within the soil matrix

Vapour flux density [kg m $^{-2}$ s $^{-1}$] (i.e., molecular diffusion flux of water vapour in the air) is defined using Fick's Law:

$$ f_v=-\rho_a D(\theta) \frac{dq}{dz}, $$

where $\rho_a$ is the density of air, [kg m $^{-3}$], $D(\theta)$ is the effective coefficient of turbulent diffusion of water vapour in pore-space air [m $^2$ s $^{-1}$], often given the form attributed to Penman (1940):

$$ D(\theta)=D_a\eta\left(\theta_s-\theta\right), $$

where $D_a$ is coefficient of molecular diffusion of water vapour in air $\approx2.12\times10^{-5}$ m $^2$ s $^{-1}$, and $\eta$ is a dimensionless coefficient characterizing the turbulent diffusion in porous media, $\eta\approx0.66$ (Penman, 1940). Note that the term $\left(\theta_s-\theta\right)$ is included correct for the volume of gas-filled porosity, where vapour flux is occurring. $q$ is the specific humidity (moisture content) of air [kg kg $^{-1}$]:

$$ q=\frac{\varepsilon e}{P}=wq^\ast(T,P), $$

where $e$ and $P$ are the vapour and total pressures, respectively [kg m $^{-1}$ s $^{-2}$], $\varepsilon$ is the ratio of the molecular weight of water to the molecular weight of air ($\approx 0.622$), $w$ is the relative humidity [--], and $q^\ast$ is the saturated specific humidity determined using the August-Roche-Magnus-Tetens relationship:

$$ q^\ast(T,P)=\frac{0.38}{P}\exp\left(\frac{17.625T}{T+243.04}\right), $$

where $P$ is in [kPa] and $T$ is in $[^\circ\text{C}]$.

Within a porous medium, assuming isobaric (i.e., $P_\text{pores}\approx P_a$) and isothermal (constant temperature) conditions, the specific humidity gradient can be related to pore-space relative humidity ($w_p$) and soil temperature ($T_s$) by:

$$ \frac{dq}{dz}=q^\ast(T_s,P)\frac{dw_p}{dz}. $$

By further assuming that the liquid and vapour phases are in equilibrium (i.e., have equal potential and the partial pressure of water has reached its vapour pressure, $e$), according to the ideal gas law, pore-space relative humidity can be determined from liquid potential by (Philip, 1957):

$$ w_p=\frac{e}{e^\ast(T_s)}=\frac{q}{q^\ast(T_s,P)}=\exp{\frac{\omega_w\psi}{RT_s}}, $$

where the gas constant $R=8.3143$ J mol $^{-1}$ K $^{-1}$, the molecular mass of water $\omega_w=0.01802$ kg mol $^{-1}$, and $e^\ast$ is the saturated vapour pressure at the temperature of the soil surface $T_s$ [K]. Next, using the chain rule,

$$ \frac{dw_p}{dz}=\frac{dw_p}{d\psi}\frac{d\psi}{dz}=w_p\frac{\omega_w}{RT_s}\frac{d\psi}{dz}, $$

Combining the above equations, results in:

$$ f_v=-K_v\frac{d\psi}{dz}, $$

where

$$ K_v=\frac{\omega_w\eta\rho_aD_aq}{RT_s}\left(\theta_s-\theta\right), $$

By combining the flux of water vapour, the Richards equation can be re-stated as (Nov'ak, 2012):

$$ \rho_l\frac{\partial\theta}{\partial t}=\frac{\partial}{\partial z}\left[K(\psi)\left(\frac{\partial\psi}{\partial z}+g\right)+\rho_a D(\theta)\frac{\partial q}{\partial z}\right], $$

and further:

$$ \rho_l\frac{\partial\theta}{\partial t}=\frac{\partial}{\partial z}\left[\left[K(\psi)+K_v(\theta,q)\right]\frac{\partial\psi}{\partial z}+gK(\psi)\right]. $$

Atmospheric exchange (evaporation-condensation)

At the soil surface--atmosphere interface, the flux density to the atmosphere remains (Nov'ak, 2012):

$$ f_v=\rho_a D \frac{dq}{dz} $$

Assuming that vertical fluxes do not change in the vicinity of the evaporating soil surface, integrating the above equation from the soil surface ($z_s$) and some effective height above the surface ($z_a$) yields: %pg.45 in Novak

$$ -\int\limits_{q_s}^{q_a}dq = \frac{f_v}{\rho_a}\int\limits_{z_s}^{z_a}\frac{dz}{D(z)}, $$

re-arranging:

$$ f_v=\frac{1}{\int_{z_s}^{z_a}\frac{dz}{D(z)}}\rho_a(q_s-q_a). $$

Next, by defining a water vapour turbulent transport coefficient ($k_v$) [m s $^{-1}$]:

$$ k_v=\frac{1}{\int_{z_s}^{z_a}\frac{dz}{D(z)}}, $$

vapour flux density to the atmosphere from the gas-filled pores can be written as:

$$ \begin{align*} f_{v,g} &=\rho_ak_v(q_s-q_a) \\ &=\rho_ak_v\left[q_s-w_aq^\ast(T_a,P_a)\right], \end{align*} $$

and

$$ f_{v,l}=\rho_ak_v(q^\ast(T_s,P_a)-q_a), $$

where $w_a$ is the atmospheric relative humidity, $T_a$ is air temperature [K], and $P_a$ is air pressure [Pa = kg m $^{-1}$ s $^{-2}$]. (Note that for water, the units for flux density [kg m $^{-2}$ s $^{-1}$] is equivalent to [mm s $^{-1}$].) The water vapour turbulent transport coefficient can be related to wind speed ($u(z)$) [m s $^{-1}$] using boundary-layer theory:

$$ k_v\approx\frac{\kappa^2u(z)}{\left[\ln\left(\frac{z}{z_0}\right)+\frac{\beta}{L_*}z\right]^2} \approx\frac{\kappa^2u(z)}{\ln\left(\frac{z-d_e}{z_0}\right)^2}, $$

where $\kappa$ is the von-K'arm'an constant ($\approx0.4$), $z_0$ is the roughness length [m], $\beta$ and $L_*$ are the Monin-Obukhov (1954) coefficient of atmospheric stability and characteristic length, respectively, and $d_e$ is the so-called zero-plane displacement height [m].

Alternatively, the Penman (1948)/Penman-Monteith (1965) relationship to aerodynamic resistance ($r_a$), soil surface resistance ($r_s$) and canopy resistance ($r_c$) [s m $^{-1}$] can be used, where:

$$ k_v=\frac{1}{r_a+r_s+r_c}. $$

Lastly, the effective mass flux evaporating from the soil surface must account for the portion of the soil surface area exposed to gas-filled pores and the area of direct liquid exposure:

$$ \begin{align*} f_e &=\left(\theta_s-\theta\right)f_{v,g}+\theta f_{v,l}. \\ &=\rho_ak_v\left[\theta_s(q_s-q_a)+\theta(q^\ast_s-q_s)\right] \end{align*} $$

Numerical solution to the Richards equation

One solution to the 1D Richards equation is the cell-centered finite volume solution scheme following Bittelli et.al. (2015). The cell-centered finite-volume form to the above equations gives the mass-balance at node $i$ as:

$$ \rho_lV_i\frac{\partial\theta_i}{\partial t}=\sum_{j=1}^{n}F_{ij}+u_i \qquad\forall i\neq j, $$

where $V_i$ is the finite volume of cell/node $i$, and the inter-nodal mass flux [kg s $^{-1}$]:

$$ F_{ij}=-A_{ij}\overline{K}i\frac{h_i-h_j}{L{ij}}, $$

and $h_i=\psi_i+gz$ is the total hydrostatic potential [J kg $^{-1}$]. In the 1D vertical, The water balance in cell $i$ (i.e., implicit Euler method---evaluating at the end of the time step) is:

$$ \frac{\rho_lV_i\overline{C}i}{\Delta t}\left(h_i^{k+1}-h_i^k\right)=F_i\left(h_i^{k+1}-h{i+i}^{k+1}\right)-F_{i-1}\left(h_{i-1}^{k+1}-h_i^{k+1}\right)+u_i, $$

where

$$ F_i=-A_i\frac{\overline{K}i}{z{i+1}-z_i}. $$

and

$$ \overline{C}_i=\frac{d\theta}{dh}\approx\frac{\theta_i^{k+1}-\theta_i^k}{h_i^{k+1}-h_i^k}. $$

The solution to the above equation is done implicitly, by solving a system equations, expanding the above and letting:

$$ \begin{align*} \varrho(i) &=\frac{\rho_l V_i\overline{C}_i}{\Delta t} \ f(i) &=\frac{A_i\overline{K}i}{z{i+1}-z_i} \ c(i) &= a(i+1) = -f(i) \ b(i) &= f(i-1)+f(i)+\varrho(i) \ d(i) &= \varrho(i)h_i^k+u_i \end{align*} $$

the finite-volume scheme for $4$ layers is put in matrix form:

$$ \begin{bmatrix} b(1) & c(1) & 0 & 0 \\ a(2) & b(2) & c(2) & 0 \\ 0 & a(3) & b(3) & c(3) \\ 0 & 0 & a(4) & b(4) \end{bmatrix} \begin{bmatrix} h^{k+1}(1) \\ h^{k+1}(2) \\ h^{k+1}(3) \\ h^{k+1}(4) \end{bmatrix} = \begin{bmatrix} d(1) \\ d(2) \\ d(3) \\ d(4) \end{bmatrix} $$

which is a tri-diagonal matrix that is solved using the Thomas algorithm.

Newton-Raphson transformation

A finite-difference solution to the 1D Richards equation is solved using the Newton-Raphson solution scheme following Bittelli et.al. (2015). The water balance in cell $i$ is:

$$ \begin{align*} F_i= & \frac{\overline{K}i\left(\psi{i+1}-\psi_i\right)}{z_{i+1}-z_i}-\frac{\overline{K}{i-1}\left(\psi_i-\psi{i-1}\right)}{z_i-z_{i-1}} \ & \quad-g\left(\overline{K}i-\overline{K}{i-1}\right)-\frac{\rho_l}{2\Delta t}\left(\theta_i^{k+1}-\theta_i^k\right)\left(z_{i+1}-z_{i-1}\right), \end{align*} $$

where $F$ is the mass balance residual and $\overline{K}$ is the mean elemental hydraulic conductivity. Applying integral transform methods (Bittelli et.al., 2015), $\overline{K}$ is defined by:

$$ \overline{K}i=\frac{K{i+1}\psi_{i+1}-K_i\psi_i}{\left(1+3/b\right)\left(\psi_{i+1}-\psi_i\right)}. $$

The solution to the above equation is then considered a minimization problem, whereby Newton-Raphson:

$$ \frac{\partial F}{\partial\psi}\cdot \left(\psi^{k+1}-\psi^k\right)=-F\rightarrow 0 \\ $$

or in a matrix-form for a $3$--layer profile is:

$$ \begin{bmatrix} \frac{\partial F_1}{\partial\psi_1} & \frac{\partial F_1}{\partial\psi_2} & \frac{\partial F_1}{\partial\psi_3} \\ \frac{\partial F_2}{\partial\psi_1} & \frac{\partial F_2}{\partial\psi_2} & \frac{\partial F_2}{\partial\psi_3} \\ \frac{\partial F_3}{\partial\psi_1} & \frac{\partial F_3}{\partial\psi_2} & \frac{\partial F_3}{\partial\psi_3} \end{bmatrix} \begin{bmatrix} \psi_1^{k+1}-\psi_1^k \\ \psi_2^{k+1}-\psi_2^k \\ \psi_3^{k+1}-\psi_3^k \end{bmatrix} = \begin{bmatrix} -F_1 \\ -F_2 \\ -F_3 \\ \end{bmatrix} $$

where $k$ indicates the $k^{th}$ iteration of the minimization attempt. The mass-balance function $F$ is solved by:

$$ F_i=f_i-f_{i-1}+u_i-u_{i-1}-\frac{\rho_l}{2\Delta t}\left(\theta_i^{k+1}-\theta_i^k\right)\left(z_{i+1}-z_{i-1}\right), $$

where

$$ f_i=-\frac{\overline{K}i\left(\psi{i+1}-\psi_i\right)}{z_{i+1}-z_i} = -\frac{K_{i+1}\psi_{i+1}-K_i\psi_i}{\left(1+3/b\right)\left(z_{i+1}-z_i\right)} $$

and letting

$$ \begin{align*} u_i &=-gK_i \\ \varrho_i &=-\rho_l\theta_i\frac{z_{i+1}-z_{i-1}}{2\Delta t}, \end{align*} $$

then substituting yields:

$$ \begin{align*} F_i= & \frac{K_i\psi_i-K_{i-1}\psi_{i-1}}{\left(1+3/b\right)\left(z_i-z_{i-1}\right)} -\frac{K_{i+1}\psi_{i+1}-K_i\psi_i}{\left(1+3/b\right)\left(z_{i+1}-z_i\right)} \\ & \quad-g\left(K_i-K_{i-1}\right)- \frac{\rho_l}{2\Delta t}\left(\theta_i^{k+1}-\theta_i^k\right)\left(z_{i+1}-z_{i-1}\right). \end{align*} $$

Differentiating $F_i$ with respect to $\psi_i$, $\psi_{i-1}$ and $\psi_{i+1}$ yields:

$$ \begin{align*} f'i &=\frac{\partial f_i}{\partial\psi_i}=\frac{K_i}{z{i+1}-z_i} \ u'i &=\frac{\partial u_i}{\partial\psi_i}=\frac{gK_i}{\psi_i}(2+3/b) \ \varrho'i &=\frac{\partial \varrho_i}{\partial\psi_i}=\rho_l\frac{\theta_i}{b\psi_i}\frac{z{i+1}-z{i-1}}{2\Delta t} \end{align*} $$

and thus

$$ \begin{align*} \frac{\partial F_i}{\partial\psi_i} &=\frac{K_i}{z_{i+1}-z_i}+\frac{K_{i-1}}{z_i-z_{i-1}}+u'i+\varrho'i \ \frac{\partial F_i}{\partial\psi{i-1}} &=\frac{-K{i-1}}{z_i-z_{i-1}}-u'{i-1} \ \frac{\partial F_i}{\partial\psi{i+1}} &=\frac{-K_i}{z_{i+1}-z_i} \end{align*} $$

Solving the Newton-Raphson system of equations is then accomplished like the linear methods described above.

Boundary conditions

Boundary conditions for the top of the vertical profile can either be a constant flux (e.g., evaporation) or constant potential:

• For constant flux, the source is added to $F_1$;
• For a constant potential, $\psi_1^{j+1}$ is known and is set at the start of the time step (since the Newton-Raphson method computes changes in $\psi$ to bring $F$ to zero, $F_1$ and $\frac{\partial F_1}{\partial\psi_i}$ are set to zero before solving; therefore the value of $\psi_1$ will remain constant throughout the time step. For infiltration, set $\psi_1=\psi_e$.

At the bottom of the profile, either a constant potential (i.e., water table) or a free drainage condition can be specified using a ghost below the profile that never becomes part of the solution. Free drainage is set by:

$$ \begin{align*} \psi_{n+1} &=\psi_n \\ \theta_{n+1} &=\theta_n \\ K_{n+1} &=K_n \end{align*} $$

after each outer iteration, where $n$ is the number of finite difference cells being solved.

Assumptions

Vapour Transport

• considers only vapour flux at the soil surface
• isothermal and isobaric throughout the time-step
• uniform soil temperature (model meant for soil surface)

Example

package main

import (
	"math"

	pm "github.com/maseology/goHydro/porousmedia"
	vsf "github.com/maseology/goVSF"
	"gonum.org/v1/plot"
	"gonum.org/v1/plot/plotter"
	"gonum.org/v1/plot/plotutil"
	"gonum.org/v1/plot/vg"
)

func main() {

	initialSe := 0.3        // initial degree of saturation
	simulationLength := 24. // hours

	dpth := []float64{0.2, 0.7, 1.} // profile depths
	pm := []pm.PorousMedium{
		{ // coarser (less clay) textured A horizon
			Ts: 0.44,
			Tr: 0.01,
			Ks: 0.000167,
			He: -2.08,
			B:  4.74,
		},
		{ // clay loam B horizon
			Ts: 0.45,
			Tr: 0.01,
			Ks: 0.000101,
			He: -5.15,
			B:  2.59,
		},
		{ // C horizon should be parent material; here sand is assumed
			Ts: 0.37,
			Tr: 0.01,
			Ks: 0.002814,
			He: -0.73,
			B:  1.69,
		},
	}

	var p vsf.Profile
	p.New(pm, dpth)
	var ps vsf.State
	ps.Initialize(p, initialSe, true)

	t, f, ok := ps.Solve(simulationLength)
	logt := make([]float64, len(t))
	for i, v := range t {
		logt[i] = math.Log10(v)
	}
	linePoints("flux.png", "infiltration rate", "log time (s)", "rate (m/s)", logt, f)
	if ok {
		w, z := ps.WaterContentProfile()
		linePoints("wcp.png", "water content", "water content (-)", "Depth below surface (m)", w, z)
	}
}

func linePoints(fp, nam, xnam, ynam string, x, y []float64) {
	p := plot.New()

	p.Title.Text = "3-layer soil moisture profile example"
	p.X.Label.Text = xnam
	p.Y.Label.Text = ynam

	err := plotutil.AddLinePoints(p, nam, points(x, y))
	if err != nil {
		panic(err)
	}

	// Save the plot to a PNG file.
	if err := p.Save(5*vg.Inch, 4*vg.Inch, fp); err != nil {
		panic(err)
	}
}

func points(x, y []float64) plotter.XYs {
	if len(x) != len(y) {
		panic("error: unequal points array sizes")
	}
	pts := make(plotter.XYs, len(x))
	for i := range pts {
		if !math.IsNaN(y[i]) {
			pts[i].X = x[i]
			pts[i].Y = y[i]
		} else {
			pts[i].X = x[i]
			pts[i].Y = 0.
		}
	}
	return pts
}

output:

References

Bittelli, M., Campbell, G.S., and Tomei, F., 2015. Soil Physics with Python. Oxford University Press.

Campbell, G.S., 1974. A simple method for determining unsaturated conductivity from moisture retention data. Soil Science, 117: 311--387.

Monteith, J.L.,1965. Evaporation and environment. Symposia of the Society for Experimental Biology 19: 205–-224.

Nov'ak, V., 2012. Evapotranspiration in the Soil-Plant-Atmosphere System. Springer. 253pp.

Penman, H.L., 1940. Gas and vapour movements in the soil: I. The diffusion of vapours through porous solids. Journal of Agricultural Science 30(3): 437--462.

Penman, H.L., 1948. Natural evaporation from open water, bare soil and grass. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 193(1032): 120--145.

Philip, J.R., 1957. Evaporation, and moisture and heat fields in the soil. Journal of Meteorology 14: 354--366.

Richards, L.A., 1931. Capillary conduction of liquids through porous media. Physics 1: 318--333.