Comparing and Contrasting Travelling Wave Behaviour for Groundwater Flow and Foam Drainage

Abstract

Liquid drainage within foam is generally described by the foam drainage equation which admits travelling wave solutions. Meanwhile, Richards’ equation has been used to describe liquid flow in unsaturated soil. Travelling wave solutions for Richards equation  are also available using soil material property functions which have been developed by Van Genuchten. In order to compare and contrast these solutions, the travelling waves are expressed as dimensionless height, ${\hat{\xi }}$, versus moisture content, $\varTheta$. For low moisture content, ${{\hat{\xi }}}$ exhibits an abrupt change away from the dry state in Richards equation compared to a much more gradual change in foam drainage. When moisture content nears saturation, ${\hat{\xi }}$ reaches large values (i.e. ${\hat{\xi }} \gg 1$) for both Richards equation and foam drainage, implying a gradual approach of $\varTheta$ towards the saturated state. The ${\hat{\xi }}$ values in Richards equation tend, however, to be larger than those in the foam drainage equation, implying an even more gradual approach towards saturation. The reasons for the difference between foam drainage and Richards equation solutions are traced back to soil material properties and in particular a soil specific parameter “m” which is determined from the soil-water retention curve.

Article Highlights

• Travelling waves for groundwater flow via Richards equation are compared with node- and channel-dominated foam drainage.

• In dry systems, liquid saturation varies abruptly with spatial location in soils (Richards equation) compared to foam.

• By contrast, in wet systems, full saturation is achieved more gradually with position in soils than in foams.

1 Introduction

Flow in soils (and foams) involves displacing air with water (or vice versa in fluid recovery processes) in the domain in which flow occurs. The laws governing fluid flow are analogous in both systems but may be expressed differently using different terminology depending on the medium considered. Despite similar physical laws (Celia et al. 1990; Richards 1931; Verbist and Weaire 1994; Verbist et al. 1996), unfamiliarity in terminology has limited exchange of ideas between the flow in soil and flow in foam research fields.

Based on work of Buckingham (1907) on capillary action in soils, Richards (1931) showed that capillary, viscous and gravity forces affected movement of water in porous media. He presented data and a general equation commonly referred to as Richards equation (hereafter RE) to describe capillary conduction of water through soil (and clay). This equation has been used extensively in groundwater hydrology, soil science, agricultural and environmental engineering (Broadbridge and White 1988; Celia et al. 1990; Patel and Pradhan 2015; Philip 1974; Parlange 1971; Raats and Van Genuchten 2006; Zlotnik et al. 2007).

Foam drainage which is governed by forces analogous to those in Richards equation was first described with an equation derived by Gol’dfarb et al. (1988). This was later advanced by Verbist and Weaire (1994); Verbist et al. (1996), who developed the so-called foam drainage equation (hereafter FDE). Their work was also based on research done by Princen and Kiss (1987) who worked on an analogous system (emulsions). The FDE as originally formulated was based on the assumption of viscous dissipation through the Plateau border channels (hereafter PB), which are channels in a foam along which three bubbles meet (Verbist and Weaire 1994; Verbist et al. 1996). Depending on surfactant type, dissipation through the foam may be predominantly via the PB (Verbist and Weaire 1994; Verbist et al. 1996), or through the nodes where PBs meet, a case treated by Koehler et al. (1999, 2000) leading to the formulation of another variant drainage equation. Both variants of the FDE are of interest in this research and both have been analysed extensively mathematically (Koehler et al. 2000; Verbist et al. 1996).

The FDE has been studied so extensively (Grassia et al. 2001; Koehler et al. 1999, 2000; Neethling et al. 2000; Verbist et al. 1996) in fact that, analysis done on it can now be used to gain insights into the behaviour of RE. Thus, the mathematics of flow of liquid in foams can be analysed alongside the flow of liquids through porous media. In particular, it is known that the FDEs have travelling wave solutions (sometimes loosely referred to as solitons (Weaire and Phelan 1996)) which can be used to predict the spatial and temporal transport of fluid within the foam (Koehler et al. 1999; Grassia et al. 2001). In the context of RE, such solutions have been deployed extensively in hydrology (Ahuja and Swartzendruber 1973; Broadbridge and White 1988; Caputo and Stepanyants 2008; Raats and Van Genuchten 2006; van Duijn et al. 2018).

As we will explain later on, the travelling wave solutions for foams and/or soils are relevant for systems with a constant liquid influx, but at sufficiently long times such that the liquid fraction at the top boundary has ceased to vary with time. They can be used to predict how much water is needed to irrigate a piece of land down to a certain depth, how rapidly water advances into the soil, as well as how rapidly an irrigation system might need to be switched off once a certain saturation is achieved at a specific depth.

A key novel contribution here will be to compare and contrast travelling wave solutions for foams and/or soils obtained over the full range of liquid contents (found numerically in the case of soils) with simpler asymptotic behaviours of those travelling waves, both for very dry systems and for systems close to saturation. In both asymptotic limits, analytical approximations for the shape of the travelling wave are available and having the analytical formula is convenient in that limit given the numerical solutions diverge there. These analytical approximations for foams and soils and for both dry and wet systems, are particularly appealing, being simple enough in their form to make it straightforward to interpret them physically. This paper thus seeks to bridge the gap between not only the mathematics, but also the physics of drainage in foams  compared to that in unsaturated soils.

 Some of the physical insights we will gain are that liquid content in soils can change very abruptly with position and/or time when soils are comparatively dry (contrasted with much more gradual changes of liquid content in foams). On the other hand, when soils become wetter, reaching the limit of full saturation is achieved very slowly indeed (more so than in foams). How abruptly or how gradually liquid content within soil changes in space and/or time may affect decisions on how to control a system used to irrigate that soil as we will explain.

This work is laid out as follows. In Sect. 2, we review concepts of liquid transport in soils and foams, and the underlying governing equations, namely Richards equation and the foam drainage equation variants. We look at the non-dimensional rescaling of the channel-dominated and node-dominated FDE variants, and the derivation and dimensionless rescaling of RE. Section 3 considers material properties of the porous media that are required within RE. Section 4 deals with equations that govern travelling wave solutions for RE. Sections 5–6 present results and discussion, and Sect. 7 concludes the paper.

2 Foam Drainage Equation and Richards Equation

As mentioned above, there is a close analogy between RE and the two FDEs. They are both continuity and mass conservation equations describing flow (in complex porous media). The forces controlling fluid flow are the same (capillary, gravity and viscous forces). In Or and Assouline (2013), an alternate approach to drainage in porous media was formulated using a variant of the FDE called the soil foam drainage equation (SFDE). Due to the complexity of obtaining accurate values for the hydraulic functions, namely the so-called hydraulic conductivity and diffusivity (Van Genuchten 1980; Van Genuchten and Nielsen 1985; Vogel and Cislerova 1988) in Richards equation, a variant of the FDE was considered treating the capillary network in soils as a network of foam Plateau borders (PBs). This is an appealing view even though unlike soils, the channels within foams expand with increasing moisture content until the foam breaks up into a bubbly liquid state, which is not the case with soils. Additionally, soils can have large local variations in their pore sizes unlike the typical situation with foams where capillary suction limits local variation of PB cross-sectional area.

Despite these structural differences, foams and soils have many similarities. In both cases, the capillary suction effects are strongest when the system is dry, but fall towards zero at full saturation. Likewise in both cases, hydraulic conduction is weakest when the system is dry, and strongest when the system is wet. Similarities like these can be exploited when comparing foam and soil systems. Solutions that have been determined for FDE (Koehler et al. 1999, 2000; Verbist and Weaire 1994; Verbist et al. 1996; Weaire and Phelan 1996) may therefore be extended to RE to gain insights into flow behaviours in soils or other porous media. This paper seeks to consider RE (Richards 1931) in close analogy to existing studies on the FDEs (Koehler et al. 1999, 2000; Verbist and Weaire 1994; Verbist et al. 1996), present travelling wave solutions for RE based on those for the FDE, and compare and contrast travelling wave solutions from RE and FDE.

We discuss the FDEs in Sects. 2.1 and 2.2 below, and RE in Sect. 2.3. To facilitate comparison, typically we consider not only the original equations, but also their dimensionless forms. We therefore consider non-dimensionalisation of the channel-dominated FDE and RE. The dimensionless node-dominated FDE is already available in literature (Koehler et al. 1999, 2000). The key dimensionless equations we obtain are Eqs. (4), (7) and (13). Readers familiar with derivations of these equations may wish to skip directly to Sect. 3.

2.1 Channel-Dominated FDE

The general form of the channel-dominated FDE proposed by Verbist and Weaire (1994); Verbist et al. (1996) is given as

$$\begin{aligned} \frac{\partial A}{\partial t} + \frac{1}{\eta }\frac{\partial }{\partial x}\left( {{\rho g}}{A^2} - \frac{C\gamma }{2} \sqrt{A} \frac{\partial A}{\partial x}\right) = 0, \end{aligned}$$

(1)

where A is area of the PB, t is time, x is position (measured downwards), $\eta = f \eta _l$ (f is a dissipation shape factor, and $\eta _{l}$ is the liquid viscosity), $\rho$ is density of liquid, g is gravity, C is a geometric shape factor and $\gamma$ is surface tension. It was found by Verbist et al. (1996) that $C\approx 0.4016$ and $f \approx 49$  (this value being an upper limit for $\eta /\eta _l$ based on the assumption of completely immobile film surfaces Leonard and Lemlich 1965).

In non-dimensionalising the channel-dominated FDE, a characteristic length-scale $x_0 = \sqrt{C\gamma /\rho {g}}$ and time-scale $t_0 = \eta /\sqrt{C \gamma \rho g}$ were chosen (Verbist et al. 1996). Although these scales can be used to non-dimensionalise the channel-dominated FDE, the dimensionless analogue of A (namely $A/x^{2}_{0}$) would not be identical to the moisture content of the foam, making it less straightforward to compare with RE. However, moisture content can be obtained from $\theta = \lambda {A}$ where $\lambda$ is the length of PB per unit volume (Brito-Parada et al. 2013; Neethling et al. 2001). Here, $\lambda$ is sensitive to bubble size (scaling as inverse square of bubble size), whereas $1/x_0^2$ is a continuum quantity (depending on $\gamma$ and $\rho$). When we introduce this factor $\lambda$, the channel-dominated FDE (1) becomes

$$\begin{aligned} \frac{\partial (\lambda {A})}{\partial t} + \frac{1}{\eta } \frac{\partial }{\partial x}\left( \frac{ {\rho g}{(\lambda {A})^2}}{\lambda } - \frac{C\gamma }{2\sqrt{\lambda }} \sqrt{\lambda {A}} \frac{\partial (\lambda {A})}{\partial x} \right) = 0. \end{aligned}$$

(2)

This form of the FDE given in Eq. (2) is the more useful form for our problem. Choosing to make the equation dimensionless using $x={x_0}\xi$ and $t = {t_0}\tau$, but now with $x_0 = ({C \gamma \sqrt{\lambda }})/({2\rho {g}})$ and $t_0 = ( C \gamma \eta {\lambda }^{3/2} )/( {2 {\rho }^2 {g}^2})$, we deduce

$$\begin{aligned} \frac{\partial \theta }{\partial \tau } - \frac{\partial }{\partial \xi }\cdot {\sqrt{\theta }}\frac{\partial \theta }{\partial \xi } + \frac{\partial \theta ^2}{\partial \xi } = 0. \end{aligned}$$

(3)

This is the dimensionless channel-dominated FDE used here. It differs slightly from a form derived by Verbist et al. (1996) (there is a factor of 2 appearing in their equation, but we have chosen $x_0$ and $t_0$ in such a fashion to scale this out to permit a fairer comparison with the node-dominated FDE (see Sect. 2.2 below)). Strictly speaking, the FDE is only valid in dry foam limit $\theta \ll 1$, but it is often considered to apply all the way up to some liquid fraction $\theta _{s}$ (saturated moisture content) at which foam breaks up into a bubbly liquid, although the geometric picture of long straight PBs leading (in the terminology of RE, see Sect. 2.3 below) to simple formulae for hydraulic conductivity and diffusivity, scaling like $\theta ^{2}$ and $\sqrt{\theta }$, respectively, does not apply for wet foams. If we assume these simple formulae apply, at least approximately up to $\theta = \theta _{s}$, we can define $\varTheta \equiv \theta /\theta _{s}$ as a rescaled liquid fraction, and after rescaling $\xi$ and $\tau$ to new variables (${\hat{\xi }} =\theta _{s}^{1/2}\,\xi$ and ${\hat{\tau }}=\theta _{s}^{3/2}\,\tau$), we obtain

$$\begin{aligned} \frac{\partial \varTheta }{\partial {\hat{\tau }}} - \frac{\partial }{\partial {\hat{\xi }}}\cdot {\sqrt{\varTheta }}\frac{\partial \varTheta }{\partial {\hat{\xi }}} + \frac{\partial \varTheta ^2}{\partial {\hat{\xi }}} = 0. \end{aligned}$$

(4)

In literature, this equation is often considered up to $\varTheta$ values as large as $\varTheta = 1$, even though formally (relying as it does upon there being a clear network of Plateau borders through which drainage can occur) the equation should be restricted to rather smaller $\varTheta$ (Weaire and Phelan 1996; Weaire and Hutzler 2001).

2.2 Node-Dominated FDE

The node-dominated FDE proposed by Koehler et al. (1999, 2000) is given as

$$\begin{aligned} \frac{\partial {\theta }}{\partial {t}} + \frac{ 2\delta _{a}L^2 }{ \eta _{l}\delta _{\theta }^{1/2}I } \left( \rho {g} \cdot \frac{\partial }{\partial {x}} {\theta }^{3/2} - \frac{\gamma \delta _{\theta }^{1/2}}{2L} \frac{\partial ^2{\theta }}{\partial {x}^2} \right) = 0, \end{aligned}$$

(5)

where $\theta$ is liquid fraction, t is time, L is length of a PB, $\delta _{a} \equiv C^2$ is a geometric shape factor, $\eta _{l}$ is fluid viscosity, $\delta _{\theta }$ is another geometric shape factor, I is a dimensionless number representative of the viscous force in the nodes and is assumed to be independent of $\theta$, $\rho$ is density of liquid, g is gravity, and $\gamma$ is surface tension. Koehler et al. (1999, 2000) reported that $\delta _{a} \equiv C^2 \approx 0.1613$, $\delta _{\theta } \approx 0.1711$ and $I\approx 400$, the value of I having been evaluated empirically, although more recent studies (Anazadehsayed et al. 2017, 2018) indicate a method to compute it.

The dimensionless form of the equation using length and time scales that are identified in Koehler et al. (1999, 2000)  ($x_{0}=(\gamma /\rho g)\delta _{\theta }^{1/2}/(2 L)$, $t_{0}=(\eta _{l} \delta _{\theta } I / (4 \delta _{a}L^{3})) \gamma /(\rho g)^{2}$) which necessarily differ from the analogous channel-dominated scales, as the parameters in the dimensional equations (2) and (5) also differ is

$$\begin{aligned} \frac{\partial \theta }{\partial \tau } - \frac{\partial }{\partial \xi }\cdot \frac{\partial \theta }{\partial \xi } + \frac{\partial \theta ^{3/2}}{\partial \xi } = 0. \end{aligned}$$

(6)

As noted for the channel-dominated case, it is possible to rewrite this in terms of a new variable $\varTheta \equiv \theta /\theta _{s}$  setting now ${{\hat{\xi }}} =\theta _{s}^{1/2}\,\xi$ and ${{\hat{\tau }}} = \theta _{s}\tau$ to obtain

$$\begin{aligned} \frac{\partial \varTheta }{\partial {\hat{\tau }}} - \frac{\partial }{\partial {\hat{\xi }}}\cdot \frac{\partial \varTheta }{\partial {\hat{\xi }}} + \frac{\partial \varTheta ^{3/2}}{\partial {\hat{\xi }}} = 0. \end{aligned}$$

(7)

Both Eqs. (4) and (7) are cast in the form of convection-diffusion equations where the convection term $\varTheta ^2$ or $\varTheta ^{3/2}$ and diffusivity term $\varTheta ^{1/2}$ or 1 is always unity when $\varTheta = 1$. This seems to enable a “fair” comparison between channel- and node-dominated theory (and compared to the formulation of Verbist and Weaire (1994); Verbist et al. (1996) is the reason for eliminating the factor of 2 that would otherwise appear in the channel-dominated FDE (4)).

2.3 Richards Equation

In Sects. 2.3.1–2.3.2, we review the derivation of Richards equation (RE) and its non-dimensionalisation, the dimensionless form facilitating comparison between RE and FDE.

2.3.1 Derivation of Richards Equation

Richards Eq. (Richards 1931) is deduced via the continuity and Darcy equations. Here,

$$\begin{aligned} {\partial {\theta }}/{\partial t} = {\partial q}/{\partial z} \end{aligned}$$

(8)

gives the continuity equation, where q is the Darcy flux which is measured downwards and distance z is measured upwards. This is a common sign convention for RE but differs from the usual convention for FDE.

Moreover, the vector flux is given as

$$\begin{aligned} {\mathbf {q}} = {K}({\hat{\mathbf {g}}} - \nabla {h}) \end{aligned}$$

(9)

where K is hydraulic conductivity, ${\hat{\mathbf {g}}}$ is the unit vector in the gravity direction and h is the capillary suction head (hereafter referred to as head). With our sign convention, $q = K((\partial h/\partial z) + 1)$. This is the Darcy equation for describing fluid flow in porous media and it is based on a potential or pressure differential. One key component is the capillary head factor which is analogous to a capillary suction term in the FDE. Substituting the expression for q into (8) gives the originally determined RE that describes fluid flow in porous media:

$$\begin{aligned} \frac{\partial \theta }{\partial t} = \frac{\partial }{\partial z} \left[ {K} \left( \frac{\partial {h}}{\partial z} + 1 \right) \right] . \end{aligned}$$

(10)

The capillary head is strictly speaking negative (i.e. it is a suction), but we can define a positive value via $h_+ = - h$. If h becomes more negative as we move downwards (i.e. as z decreases), then ${\partial h}/{\partial z}$ makes a positive contribution to the downward flux.

Richards equation as in Eq. (10) as determined above is given in a “mixed form” involving two variables, $\theta$ and h. It is more convenient to give the equation in terms of a single variable (Celia et al. 1990). Here, we opt for a so called “moisture-based form”

$$\begin{aligned} {\text{Moisture-based~ form}}\ &{\partial \theta }/{\partial t} - \nabla \cdot D(\theta )\nabla \theta - {\partial K(\theta )}/{\partial z} = 0 \end{aligned}$$

(11)

where $K(\theta )$ is hydraulic conductivity, while $D(\theta )$ is soil-water diffusivity satisfying $D\equiv K \, | \mathrm {d}h_{+}/\mathrm {d}\theta |$. As $\theta$ decreases, $h_{+}$ becomes increasingly positive, i.e. h becomes more negative. Hence, $\mathrm {d}h_{+}/\mathrm {d}\theta$ is negative, an absolute value being taken to return a positive quantity.

Equation (11) is in a form that is more readily solved than the mixed form. Regardless of the form in which the equation is written, the first term in RE represents accumulation. The second term represents capillary diffusivity, and the final term is the gravity conduction term.

Thus far, we have reviewed how to derive RE. We now consider its non-dimensionalisation.

2.3.2 Non-dimensionalisation of Richards’ Equation

The moisture-based form of the RE eliminates the head h which has units of length. There is a characteristic length scale in the soil-water retention curve relation (for h vs $\theta$) which is typically denoted as $\alpha ^{-1}$ for some parameter $\alpha$ (Van Genuchten 1980). Meanwhile, the hydraulic conductivity is $K \equiv K_s K_r(\theta )$. The multiplier $K_r(\theta )$ (relative hydraulic conductivity, hereafter RHC) is a dimensionless variable, whereas $K_s$ has the same units as K, namely $L T^{-1}$ in dimensional form and represents the hydraulic conductivity at saturation. We express $z=\xi /\alpha$ for the length scale ($\xi$ is measured positive upwards), $t=\tau /(\alpha K_s)$ for the time scale and $D(\theta )= K_s D_{r^*}(\theta )/\alpha$ for diffusivity. Here, $D_{r^*}(\theta )$ is relative diffusivity (hereafter RD).

Table 1 gives some physical properties of the soil types used in this work, helping to relate the dimensionless variables back to dimensional ones. As just mentioned, the length scale is given by $\alpha ^{-1}$, whereas $K_s$ gives a velocity scale and hence the time scale is $(\alpha K_s)^{-1}$.

Inserting these in Eq. (11), we obtain

$$\begin{aligned} \frac{\partial \theta }{\partial {\tau }} - \frac{\partial }{\partial {\xi }} \cdot {D_{r^*}(\theta )}\frac{\partial \theta }{\partial {\xi }} - \frac{\partial {K_r(\theta )}}{\partial {\xi }}= 0. \end{aligned}$$

(12)

Here relative diffusivity $D_r$ (hereafter RD) is defined such that $D = {K_s} {D_r} /\left( \alpha \left( {\theta_s} - {\theta_r} \right)\right), {\text{so}}\ {\text{that}} \ {D_r} =\left( {\theta_s}- {\theta_r }\right) {D_{r^*}}$. We may also recast this equation in the form

$$\begin{aligned} \frac{\partial \varTheta }{\partial {\hat{\tau }}}- \frac{\partial }{\partial {\hat{\xi }}} \cdot D_r(\varTheta )\frac{\partial \varTheta }{\partial {\hat{\xi }}} - \frac{\partial K_r (\varTheta )}{\partial {\hat{\xi }}} = 0, \end{aligned}$$

(13)

with $\varTheta$ indicating the rescaled moisture content

$$\begin{aligned} \varTheta = ({\theta - \theta _{r}})/({\theta _{s} - \theta _{r}}), \end{aligned}$$

where $\theta _{s}$ and $\theta _{r}$ are saturated and residual moisture contents. Note that $D_{r^*}(\theta) = K_r |{\text d}H_+/{\text d}\theta| \ {\text{where}}\ {\text{as}} \ D_r(\Theta) = K_r |{\text d}H_+/{\text d}\Theta|$. Moreover, rescaled distance and time are ${\hat{\tau }} = {\tau}(\theta _s-\theta _r)$ and ${\hat{\xi }}= \xi$.

The similarity between Eq. (13) and Eqs. (4) and (7) is clear, modulo a sign change in the final term due to measuring $\,{\hat{\xi }}$ in different directions. Comparing the above-mentioned equations, it is apparent that foam drainage analogue of diffusivity $D_r$ is $\varTheta ^{1/2}$ (channel-dominated FDE) or unity (node-dominated). The analogue of hydraulic conductivity $K_r$ is $\varTheta ^2$ (channel-dominated) or $\varTheta ^{3/2}$ (node-dominated).

Table 1 Soil-specific properties of the three example soils, with values as reported in Van Genuchten (1980). The parameter m will be discussed further in Sect. 3.1

3 Properties of Porous Media Within Richards Equation

To solve Richards Eq. (13), we need the functions $K_r(\varTheta )$ and $D_r(\varTheta )$ which are porous media properties. A number of these functions are available in literature (Brooks and Corey 1964; Van Genuchten 1980). In particular,  Van Genuchten (1980) derived hydraulic functions, namely conductivity and diffusivity from predictive conductivity models (hereafter PCM) given by Burdine (1953) and Mualem (1976) using a soil-water retention curve (which describes the suction head). Van Genuchten’s analytical expressions via the Mualem PCM (VGM) are generally considered to be the more accurate and more widely known versions (Van Genuchten 1980; Kim et al. 2008; Patel and Pradhan 2015). The PCM relates capillary suction head with water content and unsaturated hydraulic conductivity, eventually leading to capillary diffusivity. The equations that are derived (Van Genuchten 1980) depend on a parameter m which depends on the specific soil type. This model matches experimental data (Van Genuchten 1980), but the functional forms it predicts for RHC and RD turn out to be complicated, making it difficult to solve RE analytically. We deal with that issue in Sect. 4: first, however, we discuss the VGM porous media properties.

3.1 Head Function

The soil-water retention curve (SWRC) relates the rescaled water (moisture) content $\varTheta$ to the suction head h (in the soil-water retention curve) expressed below

$$\begin{aligned} \varTheta = \left[ {1 + H_+^n}\right] ^{-m} \,; \qquad H_+ = \big ( \varTheta ^{-1/m} - 1 \big )^{1/n}, \end{aligned}$$

(14)

where $H_+ = \alpha {h_+}$ is dimensionless head ($\alpha ^{-1}$ is a length scale and $h_+ = -h$), n and m are material parameters. Although fixed $m = 1$ with variable n has been used in literature (Van Genuchten 1980), usually m and n are related via $m = 1 - 1/n$ by Van Genuchten (1980) (using the predictive conductivity model (PCM) given by Mualem (1976)). Here, $n>1$, and hence, $0< m < 1$. Clayey soils are represented by comparatively low n values (m significantly less than 1), while higher n values (m close to 1) represent non-clayey soils (Stankovich and Lockington 1995). A plot of $H_+$ against $\varTheta$ for three soil samples with three different m values as determined by Van Genuchten (1980) is presented in Fig. 1.

Fig. 1

a The behaviour of dimensionless suction head ($H_+=\alpha h_{+} \equiv -\alpha {h}$) versus moisture content based on Eq. (14) for three soil types (Van Genuchten 1980), i.e. three different values of the parameter m. The inset gives a closer view of the behaviour near full saturation. b Plot of comparison of behaviour of head Eq. (14) and asymptotic behaviour (16) of head near saturation ($\varTheta \approx 1$)

3.1.1 Behaviour of Head Function in Small and Large $\varTheta$ Limits

From Fig. 1, we observe that changes in head values are very abrupt for small changes in $\varTheta$ in the limits of near residual and fully saturated moisture content limits. When head is significantly bigger than unity ($H_+ \gg 1$ or $\varTheta \ll 1$, i.e. in a very dry soil), we obtain analogously to the capillary suction head used by Brooks and Corey (1964)

$$\begin{aligned} \left. \varTheta \right| _{H_+ \gg 1} \approx H_+^{-n m}; \; \left. H_+ \right| _{\varTheta \ll 1} \approx \varTheta ^{-1/(n m)} \approx \; \varTheta ^{-(1-m)/m}. \end{aligned}$$

(15)

Thus, $H_+$ is large when $\varTheta$ becomes small, but the growth in $H_+$ is modest if $m \rightarrow 1$.

We observe also that for a near-saturated system ($\varTheta \approx 1$), the head would be much less than unity, i.e. $H_+ \ll 1$. This implies that Eq. (14) would approximate to

$$\begin{aligned} \left. \varTheta \right| _{H_+ \ll 1} \approx (1-m H_+^{1/{(1-m)}}) \,; \qquad \left. H_+ \right| _{\varTheta \approx 1} \approx \left( {(1-\varTheta )}/{m}\right) ^{1-m}. \end{aligned}$$

(16)

Figure  1b shows a comparison of original head (14) given by Van Genuchten (1980) and its asymptotic approximation (16) for $\varTheta \approx 1$. The two profiles match quite closely, and the approximation is more accurate as $m \rightarrow 1$. Note also that although $H_+ \rightarrow 0$ as $\varTheta \rightarrow 1$ as we expect, the approach is not smooth: in fact, the derivative of (14) goes to infinity in that limit

$$\begin{aligned} \left| {\mathrm {d}H_+}{/\mathrm {d}\varTheta } \right| \approx {(1-m)}{m^{-(1-m)}}(1-\varTheta )^{-m}. \end{aligned}$$

(17)

Note that for $m \rightarrow 1$, this predicts a large $| {\mathrm {d}H_+}/{\mathrm {d}\varTheta } |$ only at $\varTheta \rightarrow 1$ but much smaller $| {\mathrm {d}H_+/\mathrm {d}\varTheta } |$ elsewhere, i.e. a near step change in $H_+$ as Fig. 1b shows.

3.2 Relative Hydraulic Conductivity

Several relative hydraulic conductivity (RHC) functions have been proposed in literature (Brooks and Corey 1964; Van Genuchten 1980; Assouline and Tartakovsky 2001). Here, the more commonly used Van Genuchten–Mualem (VGM) function (Van Genuchten 1980) is employed, as has been done by Celia et al. (1990); Caputo and Stepanyants (2008); Patel and Pradhan (2015). It is given as

$$\begin{aligned} K_r(\varTheta )= \varTheta ^{1/2}[1-(1-\varTheta ^{1/m})^m]^2. \end{aligned}$$

(18)

This is based on the aforementioned predictive conductivity model (PCM) given by Mualem (1976). Another PCM (Burdine 1953) for which Van Genuchten (1980) derived another hydraulic conductivity equation is available albeit less commonly used. For low moisture contents (considered in Sects. 3.2.1 and 5.2.1), the models developed by Van Genuchten (1980) using both the Burdine and Mualem models become equivalent to RHC models given by Brooks and Corey (1964), although differences appear as $\varTheta$ grows. The Burdine case will not be discussed further here.

Fig. 2

a Plot of relative hydraulic conductivity (RHC) against moisture content for the two variants of the foam drainage equation (FDE) and three soil samples. The analogous RHC values for channel-dominated FDE (4) and node-dominated FDE (7) are $\varTheta ^2$ and $\varTheta ^{3/2}$, respectively. The dashed straight line plots the value of $\varTheta$, which helps to visualise the amount $K_r(\varTheta )$ falls below $\varTheta$. b Comparison of the true RHC Eq. (18) and its asymptotic form Eq. (20) in the $\varTheta \approx 1$ region for three soil samples. The dashed straight line at the top plots $\varTheta$, so the value of $\varTheta - K_r(\varTheta )$ is the difference between this line and the value on each curve

RHC values for different soils are plotted in Fig. 2a and compared against the FDE analogues read off from Eqs. (4) and (7).

Note that $K_r(\varTheta )=1$ when $\varTheta =1$ in all cases. However, RHC values for soils are less than those for foams, when $\varTheta < 1$. Also, the lower the m value, the lower is the RHC at any $\varTheta$.

3.2.1 Behaviour of $K_{r}(\varTheta )$ in Small and Large $\varTheta$ Limits

We examine the behaviour of the RHC function for small and large $\varTheta$ values. From Eq. (18), when $\varTheta \ll 1$, we determine

$$\begin{aligned} \left. K_r(\varTheta ) \,\right| _{\,\,\varTheta \ll 1} \approx m^2 \varTheta ^{(1/2+2/m)}. \end{aligned}$$

(19)

Consulting Eq. (19), matching power laws and  neglecting prefactors, we deduce that the node-dominated and channel-dominated FDE cases are akin to $m=2$ and $m=4/3$, respectively. Comparing these m values to those for soils (which have $0<m<1$), it leads to differences in drainage behaviour between foams and soils, at least in the small $\varTheta$ limit.

We also find the behaviour of Eq. (18) for large moisture content $\varTheta \approx 1$ values is

$$\begin{aligned} \left. K_r(\varTheta ) \,\right| _{\,\,\varTheta \approx 1} \approx \left( 1-{2}{m^{-m}}(1-\varTheta )^m \right) . \end{aligned}$$

(20)

The smaller the m value, the faster $K_r(\varTheta )$ falls as $\varTheta$ decreases.

Figure 2b compares true RHC Eq. (18) and its asymptotic form (20) for three different soils, i.e. three different m values (see Table 1). We observe that for all three m values, the asymptotic formula is a reasonable approximation. However, later on (see Sects. 4.3 and 5.2.4), we will be interested not only in values of $K_r(\varTheta )$, but also in values of $\varTheta - K_r(\varTheta )$. The values of $K_r(\varTheta )$ tend to be larger, and hence the values of $\varTheta - K_r(\varTheta )$ are smaller for the Hygiene Sandstone (and generally for cases with $m \approx 1$) than for the other soils, i.e. silt loam and Guelph loam. For $m \approx 1$ then, small errors in the $K_r(\varTheta )$ values can lead to large relative errors in $\varTheta - K_r(\varTheta )$. The implications of this will be explored further later (see Sect. 5.2.4 and the “Appendix”).

3.3 Relative Diffusivity

The general equation for determining relative diffusivity (RD) is based on the relationship for relative hydraulic conductivity (RHC) and head (soil-water retention curve, SWRC),

$$\begin{aligned} D_r(\varTheta ) = K_r(\varTheta ) \left| {\mathrm {d}H_+}/{\mathrm {d}\varTheta } \right| . \end{aligned}$$

(21)

We use Eq. (21) with Eqs. (14) and (18) to obtain the expression,

$$\begin{aligned} {D}_r(\varTheta ) = \dfrac{(1-m)}{m} \dfrac{(\varTheta ^{-1/m}-1)^{-m}}{\varTheta ^{(1/2 + 1/m)}}[1-(1-\varTheta ^{1/m})^m]^2. \end{aligned}$$

(22)

Equation (22) is plotted in Fig. 3 for three soil samples where it is compared with the analogous terms from the FDEs. We observe here that as $\varTheta \rightarrow 1$, $D_r(\varTheta ) \rightarrow \infty$ for the three soil samples (due to the behaviour of $| \mathrm {d}H_+/\mathrm {d} {\varTheta } |$ in that limit). By contrast, diffusivity is unity when the medium is fully saturated for the FDE variants (i.e. $D_r(\varTheta )=1$ at $\varTheta =1$).

Fig. 3

Profile of relative diffusivity against moisture content for the FDEs (4) and (7) and three soil samples based on Eq. (22). The relative (capillary) diffusivity terms in (4) and (7) are equivalent to $\sqrt{\varTheta }$ and 1, respectively. The functions for the soil samples are dependent on m for each soil sample

3.3.1 Behaviour of $D_r(\varTheta )$ in Small and Large $\varTheta$ Limits

We examine the behaviour of the diffusivity function for small and large $\varTheta$ values. We deduce from Eqs. (15) and (19) that when $\varTheta \ll 1$

$$\begin{aligned} \left. D_r(\varTheta ) \,\right| _{\,\,\varTheta \ll 1} \approx m\,(1-m)\,\varTheta \,^{(1/2+1/m)}. \end{aligned}$$

(23)

This Eq. (23) is plotted in Fig. 4a alongside the original Eq. (22). Hygiene Sandstone with higher m not only has a lower exponent, but also a lower prefactor. Thus, $D_r(\varTheta )$ for Hygiene Sandstone starts off larger than for Guelph loam or silt loam, but eventually it is overtaken by the two loams as $\varTheta$ increases. Note also that if we match the power law of Eq. (23) neglecting prefactors to the foam drainage relative diffusivities, the node- and channel-dominated cases now become akin to $m=-2$ (a value that would certainly not apply in the case of soils) and $m\rightarrow \infty$.

Fig. 4

Plot of relative diffusivity against moisture content for the three soil samples using Eq. (22) and asymptotic Eq. (23) in a, and (22) and asymptotic Eq. (24) in b. The solid lines are the true functions and the dashed lines are the asymptotic forms. In b, we use a zoomed view $0.99 \le \varTheta \le 1$ to facilitate comparison between (22) and (24)

We can also find a relationship for $\varTheta \approx 1$ given Eq. (22). The RHC function approaches unity in this limit. Hence, via Eq. (17)

$$\begin{aligned} \left. D_r(\varTheta )\right| _{\varTheta \approx 1} \approx \left| {\mathrm {d}H_+}/{\mathrm {d}\varTheta } \right| _{\varTheta \approx 1} \approx { (1-m) }{ m^{m-1}} (1-\varTheta )^{-m}. \end{aligned}$$

(24)

Here, $D_r(\varTheta ) \rightarrow \infty$ as $\varTheta \rightarrow 1$ as shown in Fig. 4b (a significant change from FDE for which $D_r(\varTheta ) \rightarrow 1$ as $\varTheta \rightarrow 1$). Here, Hygiene Sandstone has the larger (negative) exponent but the smaller prefactor. Thus, it has the largest $D_r(\varTheta )$ as $\varTheta \rightarrow 1$, but not as $\varTheta$ starts to fall. This implies (Fig. 4 ) that $D_r(\varTheta )$ is non-monotonic in m, depending on the $\varTheta$ value.

This completes our presentation of soil sample material properties.

Good estimates of the SWRC are essential (Van Genuchten 1980; Vogel and Cislerova 1988; Assouline et al. 1998), since the SWRC is used in the Mualem (1976) predictive model to obtain the RHC, and then both the retention curve and RHC are used to obtain the RD. In what follows, this information is used to obtain moisture content profiles.

4 Travelling Wave Solution

Various different types of solution have been studied for foam drainage, among them the so-called free drainage (in which an initially comparatively wet foam loses water from its bottom boundary) and forced drainage (in which water is added to the top of what is initially a comparatively dry foam) (Verbist et al. 1996; Koehler et al. 2000; Neethling et al. 2001).

Given the close analogy between foam drainage and Richards equation, similar scenarios are expected in soils: e.g. rainwater almost saturating soil and then draining down to the water table and/or irrigating an initially almost dry soil. Here we focus on travelling wave solutions that are relevant to the forced drainage case, as these are rather simpler to handle mathematically than free drainage is (Verbist et al. 1996; Cox et al. 2000).

Travelling wave solutions for RE have previously been proposed (Philip 1957a; Gilding 1991; Witelski 1997). In this section, we consider one such travelling wave solution for RE, analysing it specifically in the context of mathematical models for soil material properties (Burdine 1953; Mualem 1976; Van Genuchten 1980) described in Sect. 3. We consider a general expression in Sect. 4.1 for the travelling wave obtained from RE. We then deduce the velocity of the travelling wave in Sect. 4.2. Finally, we determine the equation that describes the shape of the travelling wave in Sect. 4.3.

4.1 General Form of Travelling Wave

From the general form of a travelling wave, supposing distance (depth) ${{\hat{\xi }}}$ is positive upwards and velocity $\upsilon$ is positive downwards, we search for a solution of the form

$$\begin{aligned} \varTheta = \varTheta \,({\hat{\xi }} + \upsilon {\hat{\tau }}). \end{aligned}$$

(25)

As length and time only appear in the combination of ${\hat{\xi }} + \upsilon {\hat{\tau }}$, we can express the solution of RE as a function of position rather than time, derivatives being related via $\partial /\partial {\hat{\tau }} = \upsilon \cdot \partial /\partial {\hat{\xi }}$. Integrating (13) with respect to ${\hat{\xi }}$ yields

$$\begin{aligned} \upsilon \varTheta - D_r(\varTheta ) \cdot {\partial \varTheta }/{\partial {\hat{\xi }}} - K_{r}(\varTheta ) = \mathrm {constant}. \end{aligned}$$

(26)

We impose boundary conditions that a long way upstream $\varTheta \rightarrow \varTheta _{1}$, and a long way downstream $\varTheta \rightarrow \varTheta _{2}$, with $\varTheta _{1} > \varTheta _{2}$ (i.e. high saturation at the top and low saturation at the bottom).

The wave velocity $\upsilon$ now obeys the Rankine–Hugoniot condition (Philip 1957b),

$$\begin{aligned} {\upsilon } = {(K_{r}(\varTheta _1) - K_{r}(\varTheta _2))}{/(\varTheta _1 - \varTheta _2)}. \end{aligned}$$

(27)

After some algebra, we obtain the constant on the right hand side of (26) as

$$\begin{aligned} \mathrm {constant} = {(K_{r}(\varTheta _1)\varTheta _2 - K_{r}(\varTheta _2)\varTheta _1)}/{(\varTheta _1 - \varTheta _2)}. \end{aligned}$$

(28)

Note that in the limit of a dry system at the bottom $\varTheta _{2} \rightarrow 0$, and Eq. (28) vanishes.

4.2 Evaluating the Travelling Wave Propagation Velocity

To study the travelling wave propagation velocity, we plot Eq. (27) for various moisture contents. We use the Van Genuchten–Mualem RHC relationships in Fig. 5 and plot $\upsilon$ vs $\varTheta _1$ at different initial moisture contents $\varTheta _2$ where $\varTheta _{2} < \varTheta _{1} \le 1$ ($\varTheta _{2} = 0, \,0.5$).

In the special case $\varTheta _{2} = 0$, Eq. (27) reduces to $K_r(\varTheta _1)/\varTheta _{1}$ (as in Fig. 5a).

Fig. 5

Travelling wave propagation velocity for foam drainage equations (FDE) and three soil samples from Eq. (27) for a case where a moisture content at the bottom $\varTheta _{2}=0$ and b moisture content at the bottom $\varTheta _{2}=0.5$, i.e. in an already wet system

We see in Fig. 5 that $\upsilon$ increases with increasing $\varTheta _1$ (and does so abruptly near $\varTheta _{1} = 1$ in the case of soils). Additionally, we observe that $\upsilon$ increases with increasing $\varTheta _{2}$. When $\varTheta _{2}$ is non-zero (i.e. in an already wet system with $\Theta_2=0.5$), the soil samples tend to have lower velocity than the foam drainage cases for most $\varTheta _1$ values, but it actually have higher velocity when $\varTheta _{1} \rightarrow 1$.

4.3 Profile of ${\hat{\xi }}$ Versus $\varTheta$

It turns out that travelling wave solutions satisfying Eqs. (26) and (27) only exist if $\varTheta _{1} > \varTheta _{2}$. Inserting Eqs. (27) and (28) in (26), and rearranging, we obtain

$$\begin{aligned} \frac{\mathrm {d} {\hat{\xi }}}{\mathrm {d}\varTheta } =\frac{D_r(\varTheta )(\varTheta _1 - \varTheta _2)}{K_{r}(\varTheta _1)(\varTheta - \varTheta _2) - K_{r}(\varTheta _2)(\varTheta - \varTheta _1)- K_{r}(\varTheta )(\varTheta _1 - \varTheta _2)}. \end{aligned}$$

(29)

In a special case $\varTheta _1 = 1$ and $\varTheta _2 = 0$ (the upper and lower limits possible for a travelling wave), Eq. (29) simplifies. In this case, $K_{r}(\varTheta _1)=1$, $K_{r}(\varTheta _2)=0$, $\upsilon = 1$, and the constant on the right hand side of Eq. (26) vanishes as presented by Parlange (1971); Witelski (1997). We deduce

$$\begin{aligned} {\mathrm {d} {\hat{\xi }}}/{\mathrm {d}\varTheta } = {D_r(\varTheta )}/ \left( {\varTheta - K_{r}(\varTheta )} \right) . \end{aligned}$$

(30)

Note $\varTheta - K_{r}(\varTheta )$ in the denominator: the relevance of this term was stated in Sect. 3.2.1.

The special case (30) is analysed further in Sect. 5. Equation (30) applies equally well in the case of RE and FDE. It is only the form of the RHC and RD functions that varies. For certain forms of RHC and RD, analytical solutions to (30) are obtained (Verbist et al. 1996; Koehler et al. 2000; Cox et al. 2000). In cases for which choices of RHC and RD do not allow analytical solution of (30), we compute ${\hat{\xi }}$ against $\varTheta$ via Simpson’s rule, implemented in MATLAB using a $\varTheta$ step size of 0.0001.

5 Results for Travelling Wave Profiles

In what follows, we focus on numerical and asymptotic analytical solutions obtained from Eq. (30). Specifically, we compare and contrast solutions to this equation using the parameters from the FDEs and also using soil material properties derived by Van Genuchten (1980) for porous soils in RE. In Sect. 5.1, we review travelling wave solutions to the two variants of the FDE (Verbist et al. 1996; Koehler et al. 1999, 2000) and their asymptotic behaviours in the limit of small and large moisture content ($\varTheta$). These FDE solutions are known from literature (Verbist et al. 1996; Koehler et al. 1999, 2000; Cox et al. 2000; Weaire and Hutzler 2001) allowing us to benchmark our approach. Subsequently, we consider the travelling wave solution to RE in Sect. 5.2, comparing back to the FDE. After that, Sect. 5.3 then considers an integrated quantity that we call “missing moisture”.

5.1 Foam Drainage Equations Solutions

In solving for FDE travelling wave solutions analogously to Eq. (30), we use the parameters describing hydraulic conductivity and capillary diffusivity for the two FDEs from Eqs. (4) and (7). The solutions correspond to travelling waves previously determined in literature (Verbist et al. 1996; Koehler et al. 2000), remembering however from Sect. 2.1 that there is a factor of 2 difference scaled out of Eqs. (3)–(4). Scaling out this factor gives a fairer comparison between the two FDEs as they have the same diffusivity and conductivity at full saturation.

The equation obtained in the case of channel-dominated FDE is given as

$$\begin{aligned} {\mathrm {d} {\hat{\xi }}} / {\mathrm {d}\varTheta } =1/ \big ({\sqrt{\varTheta } \;(1 - \varTheta )} \big ), \end{aligned}$$

(31)

the solution of which is given as

$$\begin{aligned} {\hat{\xi }} = 2 \, \, \mathrm {arctanh} \,\sqrt{\varTheta }\, + \,{c_{CD}} \end{aligned}$$

(32)

where $c_{CD}$ is an integration constant and can be set to zero. Likewise, we obtain for the node-dominated case

$$\begin{aligned} {\mathrm {d} {\hat{\xi }}} / {\mathrm {d}\varTheta } = 1/\big ({\varTheta \; (1 - \sqrt{\varTheta }) }\big ). \end{aligned}$$

(33)

Integrating (33), we obtain

$$\begin{aligned} {\hat{\xi }} = {2} \log \big ( {\sqrt{\varTheta }}/{(1 - \sqrt{\varTheta })} \big ) + {c_{ND}} \end{aligned}$$

(34)

where $c_{ND}$ is an integration constant that again we set equal to zero.

Equations (32) and (34) are plotted in Fig. 6a assuming $c_{CD} = c_{ND} = 0$. In the case of the node-dominated FDE, we observed solutions exhibiting inflection points with ${\hat{\xi }}$ extending to both positive and negative infinity. The channel-dominated solution also exhibits an inflection point, but the moisture content $\varTheta$ goes to zero at finite depth unlike the node-dominated case where $\varTheta$ only vanishes as ${\hat{\xi }} \rightarrow \infty$.

As $\varTheta \rightarrow 1$, Fig. 6a shows that ND FDE has larger ${\hat{\xi }}$ than CD FDE does. This follows because ${\mathrm {d} {\hat{\xi }}/\mathrm {d}\varTheta }$ is larger according to Eq. (30). Both ND FDE and CD FDE have a common RD value as $\varTheta \rightarrow 1$ but their $\varTheta - K_r(\varTheta )$ value in Eq. (30) can be approximated by $(1-\varTheta )({K_r}'(1)-1)$, where ${K_r}'$ denotes ${\mathrm {d} K_r/\mathrm {d} \varTheta }$. Figure 2a shows that ND FDE has smaller ${K_r}'(1)$ than CD FDE hence larger ${\mathrm {d} {\hat{\xi }}/\mathrm {d}\varTheta }$, leading in turn to larger ${\hat{\xi }}$ values.

5.1.1 Asymptotic Behaviour of the CD FDE

We analyse the FDEs in the limits of either $\varTheta \rightarrow 0$ or $\varTheta \rightarrow 1$ since that will form an interesting contrast from RE.  Considering Eq. (31) when $\varTheta \ll 1$, and integrating (following an approach of Witelski (1997)) we deduce

$$\begin{aligned} {\hat{\xi }} \approx 2 \sqrt{\varTheta } + {c_{CD}}, \end{aligned}$$

(35)

where we set the integration constant $c_{CD} =0$. We note that $\varTheta \rightarrow 0$ at finite ${\hat{\xi }}$ but nonetheless $\mathrm {d}{\hat{\xi }}/\mathrm {d}\varTheta$ is infinite as $\varTheta \rightarrow 0$ (so $\varTheta$ changes only slowly as ${\hat{\xi }}$ changes).

We can also consider the case when $\varTheta \approx 1$. Taking Eq. (31) in this limit and integrating (again using the approach of Witelski (1997)) we find

$$\begin{aligned} {\hat{\xi }} \approx \log \left( {1}/{(1 - \varTheta )} \right) + {c_{CD1}} \equiv \log \left( {1}/{(1 - \varTheta )} \right) + \log 4, \end{aligned}$$

(36)

where the integration constant $c_{CD1}$ is set to $c_{CD1}=\log 4$ to match with (32). This follows because $2 \, \mathrm {arctanh} \,(\sqrt{\varTheta })$ within (32) asymptotes to $\log (4/(1-\varTheta ))$ as $\varTheta$ approaches 1.

5.1.2 Asymptotic Behaviour of the ND FDE

We estimate the behaviour of the node-dominated FDE when the moisture content is very small ($\varTheta \ll 1$). Taking (33) in this limit and integrating gives

$$\begin{aligned} {\hat{\xi }} \approx \log (\varTheta ) + {c_{ND}} \end{aligned}$$

(37)

where we can set $c_{ND} = 0$. We can likewise obtain an expression for the node-dominated FDE in the limit $\varTheta \approx 1$, taking Eq. (33) for $\varTheta \approx 1$ and integrating gives

$$\begin{aligned} {\hat{\xi }} \approx 2 \log \left( {1}{/(1 - \varTheta )} \right) + {c_{ND1}} \equiv 2 \log \left( {1}/{(1-\varTheta )} \right) + \log 4 \end{aligned}$$

(38)

where the value of the constant $c_{ND1} = \log 4$ is set to match with (34). This follows because $2 \log (\sqrt{\varTheta } / (1-\sqrt{\varTheta }) )$ within (34) asymptotes to $2 \log (2/(1-\varTheta ))$ as $\varTheta$ approaches 1.

Comparing Eq. (38) with (36) confirms the deduction that compared to the channel-dominated system, the moisture content in the node-dominated system makes a slower approach as $\varTheta \rightarrow 1$, i.e. it requires a larger ${\hat{\xi }}$ to achieve a given $\varTheta$. Having now reviewed the FDE solutions from Verbist et al. (1996); Koehler et al. (2000), in what follows, we contrast them with solutions of RE.

Fig. 6

a Profile of travelling wave solution for node-dominated and channel-dominated solutions. The inset shows a zoomed view of the channel-dominated case. b Profile of travelling wave solution (using Simpson’s rule) for Richards Eq. (39) using the Van Genuchten soil material functions (18) and (22) for three soil types. c The same solution (as b) on semi-log axes

5.2 Van Genuchten Solution to Richards Equation

We seek travelling wave solutions for Richards equation (30) using soil functions derived by Van Genuchten (1980) (see Sect. 3). Substituting Eqs. (18) and (22) into (30) gives

$$\begin{aligned} \frac{\mathrm {d} {\hat{\xi }} }{ \mathrm {d} \varTheta } = \dfrac{(1-m) \cdot \varTheta ^{-1/m-3/2} [1 - (1-\varTheta ^{1/m})^{m}]^{2}}{m \cdot (\varTheta ^{-1/m}-1)^{m} \left( 1-\varTheta ^{-1/2}[1 - (1-\varTheta ^{1/m})^{m}]^{2}\right) }, \end{aligned}$$

(39)

which we integrate by Simpson’s rule. The solution for ${\hat{\xi }}$ is shown in Fig. 6b. For low $\varTheta$, there is a very sharp change in $\varTheta$ for a small change in ${\hat{\xi }}$ (a contrast from the FDE solutions). Figure 6c shows a semi-log profile of the solution to (39), making it easier to visualise behaviour at small ${\hat{\xi }}$ in particular. Looking towards larger $\varTheta$ values ($\varTheta \rightarrow 1$), however, there are very large values of ${\hat{\xi }}$. How these behaviours come about is considered below.

5.2.1 Asymptotic Behaviour of the Van Genuchten Solution

When $\varTheta \ll 1$, the dominant material function affecting the profile shape is diffusivity. Equation (30) reduces to

$$\begin{aligned} \left. {\mathrm {d} {\hat{\xi }}}/{ \mathrm {d}{\varTheta }} \right| _{\,\,\varTheta \ll 1} \approx {D_{r}(\varTheta )}/{\varTheta }, \end{aligned}$$

(40)

where the denominator is estimated as being approximately $\varTheta$, neglecting $K_r(\varTheta )$ relative to $\varTheta$ in small moisture content limits. Substituting from Eq. (23) and integrating gives (as found by Witelski (1997))

$$\begin{aligned} \left. {\hat{\xi }} \right| _{\,\,\varTheta \ll 1} \approx \frac{2m^2\,(1-m)}{(2+m)} \varTheta ^{(1/2 \, + \, 1/m)} + {c_{VGM}}, \end{aligned}$$

(41)

where the constant ${c_{VGM}}$ is set to zero to ensure ${\hat{\xi }}=0$ as $\varTheta \rightarrow 0$. This solution is shown in Fig. 7 via dashed lines for each soil type (m value). Noting the scale on the vertical axis, we observe that $\varTheta$ changes very abruptly with ${\hat{\xi }}$, or equivalently ${\hat{\xi }}$ changes only very slowly with $\varTheta$. Changes in ${\hat{\xi }}$ are slowest when m is small (owing to the power law in $\varTheta ^{(1/2+1/m)}$ in (41)), but are also surprisingly slow for Hygiene Sandstone ($m=0.9038$) compared to Guelph loam, owing to the prefactor $(1-m)$ in Eq. (41) which vanishes as $m \rightarrow 1$. As was the case in Fig. 4 for $D_r(\varTheta )$, the behaviour here for ${\hat{\xi }}$ is non-monotonic in m.

These slow changes in ${\hat{\xi }}$ (or equivalently abrupt changes in $\varTheta$) predicted by (41) differ from what is seen in the two FDEs in the $\varTheta \ll 1$ limit. The implication for a case in which a sensor had been set to stop irrigating a soil when $\varTheta$ reached some initial value, e.g. $\varTheta = 0.1$, is that substantially higher $\varTheta$ values $\varTheta = 0.2$, $\varTheta = 0.3$, etc., will follow just shortly thereafter.

5.2.2 Comparison of True and Asymptotic Solution ($\varTheta \ll 1$)

In Fig. 7, we compare solutions from the true (i.e. numerical) solution in Eq. (39) and its asymptotic version Eq. (41) to verify the accuracy of our analysis in Sect. 5.2.1. We observe that the asymptotic solution can barely be distinguished on the scale of the graph for silt loam ($m = 0.5146$) at least up to $\varTheta \approx 0.2$,  the ratio of the true to asymptotic solution at $\varTheta = 0.2$ being 1.0256. The accuracy of the asymptotic solution reduces as m increases. For $m=0.6377$ (Guelph loam) the ratio of the true to asymptotic solution at $\varTheta =0.2$ becomes 1.0510, whereas for $m=0.9038$ (Hygiene Sandstone) it is 1.1463. Nonetheless, the scale of the graph (up to less than ${\hat{\xi }} \approx 10^{-3}$) shows we are dealing with tiny ${\hat{\xi }}$ values throughout.

5.2.3 Asymptotic Behaviour when $\varTheta \approx 1$

When we consider the case $\varTheta \approx 1$, we can approximate the behaviour of Eq. (30) as

$$\begin{aligned} \left. \frac{ \mathrm {d} {\hat{\xi }} }{ \mathrm {d}\varTheta } \,\, \right| _{\,\,\varTheta \approx 1} \approx \,\,\frac{ |\mathrm {d}H_+ / \mathrm {d}\varTheta | }{1-K_r(\varTheta )} , \end{aligned}$$

(42)

where the numerator of (30) has been approximated by $|\mathrm {d}H_+ / \mathrm {d}\varTheta |$ (since $K_r(\varTheta ) \rightarrow 1$) and where the denominator $\varTheta - K_r(\varTheta )$ is expressed as $(1 - K_r(\varTheta )) - \,(1-\varTheta ) \,$, with $(1-\varTheta ) \ll (1 - K_r(\varTheta ))$ (see Eq. (20) with $m < 1$ here). It follows from (20) that $\mathrm {d}K_r(\varTheta )/\mathrm {d}\varTheta$ diverges at $\varTheta \rightarrow 1$, explaining the rapid decrease in $K_r(\varTheta )$ with decreasing $\varTheta$ for soils in Fig. 2a–b.  Inserting equations (16), (20) and (24) in (42), and integrating

$$\begin{aligned} \left. {\hat{\xi }} \,\, \right| _{\,\,\varTheta \approx 1} \approx \frac{ (1-m) m^{(2m-1)}}{2(2m-1)}(1-\varTheta )^{1-2m} + {c_{VGM1}}, \end{aligned}$$

(43)

where ${c_{VGM1}}$ can only be obtained by matching to the solution of Eq. (39) via Simpson’s rule. We have arbitrarily chosen to match at $\varTheta = 0.9$. The dashed lines in Fig. 7b are the asymptotic solutions given in (43). Out of the soils considered, Hygiene Sandstone (with the largest m) has the biggest ${\hat{\xi }}$ scaling as $(1-\varTheta )^{1-2m}$, and this is also bigger as $\varTheta \rightarrow 1$ than the $\log {(1/(1-\varTheta ))}$ terms obtained in the FDEs (see Eqs. (36) and (38)).

Fig. 7

Comparison of profiles for true or numerical function (39) and the asymptotic function (41) in a; and between (39) and the asymptotic function (43) in b

It is interesting to understand how these differences between FDE and RE come about. In Eq. (30), the denominator $\varTheta - K_r(\varTheta )$ vanishes in both the FDE case and the RE case but does so more slowly for RE. On the other hand, the numerator $D_r(\varTheta )$ is finite as $\varTheta \rightarrow 1$ for FDE but divergent for RE, and this more than compensates the more slowly vanishing numerator. The net result is that ${\hat{\xi }}$ grows more rapidly for soils than for foams.

The implication of (43) for an irrigation system is that if a sensor were to detect when the saturation at a certain location reached a certain critical level (e.g. say $\varTheta = 0.9$), it might be quite some time in the case of Hygiene Sandstone before substantially higher saturations, e.g. $\varTheta = 0.99$ or $\varTheta = 0.999$ arrive at the sensor. If the sensor reaching $\varTheta = 0.9$ is a signal that an irrigation system might eventually need to be switched off, and the higher level $\varTheta = 0.99$ or $\varTheta = 0.999$ is the level at which switch off is actually required, there could in the case of Hygiene Sandstone be a delay between the original signal and the eventual switch off.

5.2.4 Comparison of True and Asymptotic Solution ($\varTheta \approx 1$)

We again compare solutions from the numerical solution of RE and its asymptotic approximation but now for $\varTheta \rightarrow 1$. Specifically, Fig. 7b shows a comparison of the true solution from (39) (obtained via Simpson’s rule) and the asymptotic solution (43) (with a value matched to the numerical solution at $\varTheta = 0.9$). For each m, the two profiles (true profile vs asymptotic) show slight differences. For $m = 0.9038$ (Hygiene Sandstone), the asymptotic profile is an underestimate, whereas for $m = 0.5146$ (silt loam) and $m = 0.6377$ (Guelph loam) it is an overestimate. To highlight this, we take a ratio of ${\hat{\xi }}$ between the true (numerical) and asymptotic solutions as shown in Fig. 8.

 On the domain of $\varTheta$ plotted, the ratio reaches a maximum value of 1.1456 in the case of Hygiene Sandstone and reaches minimum values of 0.7981 & 0.9072 in the case of silt loam and Guelph loam, respectively.

Fig. 8

Profile of ratio of solution to equation for true ${\hat{\xi }}$ (Eq. (39)) to approximate ${\hat{\xi }}$ (Eq. (43))

We now examine the different behaviour in Fig. 8 for different m, identifying terms which lead to such behaviour by decomposing the governing equations (30) and (42). Differences could arise from the numerator ($D_r(\varTheta )$ itself broken down into two components, i.e. the product of $| \mathrm {d}H_+/\mathrm {d}\varTheta |$ and $K_r (\varTheta )$). Alternately, differences could arise from the denominators, respectively ($\varTheta - K_r(\varTheta )$) and ($1 - K_r(\varTheta )$) as shown in (30) and (42).

Fig. 9

Ratio of a true $\mathrm {d}H_+/\mathrm {d}\varTheta$ used in Eq. 30 (via Eq. (21)) to the approximate $\mathrm {d}H_+/\mathrm {d}\varTheta$ in equation (42) (via Eq. 17), and b the denominator in Eq. (30) to the asymptotic version in Eq. (42) (via Eq. 20)

Figure 9 shows considerable variation between the original and the asymptotic $| \mathrm {d}H_+/\mathrm {d}\varTheta |$ functions. The true $| \mathrm {d}H_+/\mathrm {d}\varTheta |$ obtained via (14) exceeds the asymptotic one given by (17), by up to $25\%$ at some points in the domain $0.9< \varTheta < 1$ in the case of silt loam. The ratio of true $K_r(\varTheta )$ in the original Eq. (18) to the asymptotic $K_r(\varTheta )$ (which can be taken equal to 1 when $\varTheta \rightarrow 1$) is the same as the solid lines shown in Fig. 2b. Whereas, the original $\mathrm {d}H_+/\mathrm {d}\varTheta$ is bigger than the asymptotic one, the original $K_r(\varTheta )$ is smaller to compensate, falling to less than half the asymptotic value in the case of silt loam at $\varTheta =0.9$. As $|\mathrm {d}H_+ / \mathrm {d}\varTheta |$ and $K_r(\varTheta )$ are multiplied together to form $D_r(\varTheta )$, the dominant effect is the original $K_r(\varTheta )$ being less than the asymptotic one (crudely approximated by unity here) and this is reflected in the ${\hat{\xi }}$ versus $\varTheta$ profiles for silt loam and Guelph loam in Fig. 7b and  also in the ratios in Fig. 8. The effect is less pertinent for Hygiene Sandstone since its $K_r(\varTheta )$ does not decrease as sharply with decreasing $\varTheta$ as happens for the other soils.

We now consider the denominators of (30) and (42). Figure 9b shows that using the asymptotic expression for $1-K_r(\varTheta )$ in lieu of $\varTheta - K_r(\varTheta )$ gives some deviation (about $30\%$) for the soils over the domain $0.9< \varTheta < 1$. It is accurate (to within $10\%$) at $\varTheta \approx 0.99$ for both loams. For Hygiene Sandstone however, the asymptotic $1 - K_r(\varTheta )$ exceeds the true $\varTheta - K_r(\varTheta )$ significantly even for $\varTheta \approx 0.99$. When m is close to 1, it appears that the assumption employed in (42) that $1 - \varTheta \ll 1 - K_r(\varTheta )$ is no longer applicable, for reasons that are explained in the “Appendix”.

The result is that for Hygiene Sandstone in particular, equation (30) has a smaller denominator than the asymptotic prediction even for $\varTheta$ really close to 1, hence the predicted ${\hat{\xi }}$ via Eq. (30) is larger as Figs. 7b and 8 show.

To summarise, numerical (i.e. Simpson’s rule) solutions generated for Richards equation have been examined. Analytical asymptotic approximations for $\varTheta \ll 1 \; \& \; \varTheta \approx 1$  captured the true solution behaviour reasonably well within $0<\varTheta <0.25$ and $0.9<\varTheta <1$, respectively, although the level of agreement is sensitive to a parameter m. This parameter m reflects the pore structure (capillary model) for each soil type (Van Genuchten 1980). Due to variability of pore structure and how it affects soil-water retention, the value of m can vary from one soil type to another (clayey soils have smaller m), and this is reflected in the moisture content profile.

5.3 Integral of ${\hat{\xi }}$

There is one final quantity we consider, namely $\int _{0}^{1} {\hat{\xi }}\, \mathrm {d}\varTheta$. For a front that has reached a given location, setting ${\hat{\xi }} = 0$ at the point where $\varTheta = 0$, the value $\int _{0}^{1} {\hat{\xi }} \, \mathrm {d}\varTheta$ is a measure of the amount of additional liquid that would be needed above the front to saturate fully the liquid above. In other words, it is the “missing” moisture that is required to saturate fully or flood the soil down to the current depth of the front. The integral is easily computed by dividing up the domain $0\le \varTheta \le 1$ into three subdomains: $\varTheta _{1}'\le \varTheta \le 1$ (for some $\varTheta _{1}'$ close to unity), $\varTheta _{2}'\le \varTheta \le \varTheta _{1}'$ (for some $\varTheta _{2}'\ll 1$), and $0\le \varTheta \le \varTheta _{2}'$. Integration is done analytically using asymptotic formulae for ${{\hat{\xi }}}$ in the first and third subdomain and numerically via Simpson’s rule in the second subdomain. We chose $\varTheta _{1}'=0.1$ and $\varTheta _{2}'=0.9$. The missing moisture evaluates to 0.2243 for Hygiene Sandstone ($m = 0.9038$), 0.1204 for Guelph loam ($m = 0.6377$) and 0.0808 for silt loam ($m = 0.5146$). Based on this, loam is more readily flooded than sandstone. By comparison for the channel-dominated FDE, the value of $\int _{0}^{1} {\hat{\xi }}\, \mathrm {d}\varTheta$ is 2. Hence, the amount of missing moisture is larger for foam than for the soil.

 Even though the ${\hat{\xi }}$ values for soils approach ${\hat{\xi }} \rightarrow \infty$ more rapidly as $\varTheta \rightarrow 1$ than the CD FDE does (i.e. a power law Eq. (43) rather than a logarithm Eq. (36)), the integral $\int _{0}^{1} {\hat{\xi }}\,\mathrm {d}\varTheta$ is greater for foam. This is because the channel-dominated case has significant ${\hat{\xi }}$ values even for smaller $\varTheta$, away from the neighbourhood of $\varTheta = 1$.

It is meaningless to compute $\int _{0}^{1} {\hat{\xi }}\, \mathrm {d}\varTheta$ for the node-dominated FDE as there is no finite ${\hat{\xi }}$ location at which $\varTheta$ falls to zero so there is a degree of arbitrariness in how much ${\hat{\xi }}$ is shifted up or down.

6 Discussion

In the previous section, we studied in detail the travelling wave solutions of Richards equation and foam drainage, focussing on asymptotic behaviour in various ranges of liquid saturation. We have not, however, discussed whether and if so when these solutions are relevant in real world infiltration problems.

The travelling wave solutions presented in this paper are long time solutions and thus are by themselves insufficient to describe real infiltration problems, since travelling waves take some time to develop under standard boundary conditions (Broadbridge and White 1988). In order to approach close to a travelling wave solution, depths of at least several dimensionless units have to be reached, or converting to dimensional lengths, several times $\alpha ^{-1}$ (see Table 1). Additionally, infiltration needs to proceed long enough for that to happen, i.e. several units of dimensionless time, which is easily related to dimensional time ($(\alpha K_s)^{-1}$ being the typical time scale in Table 1). An unsteady state simulation of liquid infiltration into soils could be safely ceased as soon as the simulated front shape matched the travelling wave, and it would continue to follow the travelling wave thereafter.

Whereas, long time solutions assume a given liquid saturation moving arbitrarily far up from the travelling wave front, in nature, what is set is the infiltration rate (rate at which liquid enters the soil (Parlange 1972; Broadbridge and White 1988), or analogously (for the FDE) enters the foam (Brito-Parada et al. 2013)) rather than setting liquid saturation at the top. In the long time limit, where the gradient of liquid saturation is negligible near the top, all infiltration at the top is due to gravity. Hence, setting the infiltration rate at the top defines the liquid saturation and vice versa (in the long time limit).

For shorter times, however, when there are still significant gradients of liquid saturation near the top, it is not necessary to have a high liquid saturation right at the top to achieve the same infiltration rate, since only part of the infiltration is contributed by gravity, the rest being contributed by  capillarity. It is possible to consider early-time solutions to Richards equation, in which liquid saturation at the top grows with time, rather than being fixed. There turns out to be a similarity solution (Broadbridge and White 1988; Witelski 1997; Caputo and Stepanyants 2008) in which liquid saturation is mostly uniform with height but shows a slight increase near the top localised over a small height.

The amount of the slight increase in liquid saturation gradually increases over time, as does the height over which it is registered. However, solutions at one time can be overlain with solutions at a different time by a suitable rescaling. In other words, the solutions at different times are self-similar provided we remain within the early-time limit. Note that because we are dealing with just slight increases over an initial liquid saturation in this limit, we should not need to know soil material properties over the entire range of liquid saturations, just the local behaviour near the initial saturation. We will not consider early-time similarity solutions any further in this paper having chosen to leave that aspect for future work, so we focussed instead on late-time travelling wave solutions for which porous media properties over the full range of saturations become pertinent.

We note however that early-time similarity solutions are found in broadly analogous multiphase flow systems (e.g. solid-liquid suspensions). Early on during gravity settling of a suspension, the suspension is mostly uniform, but the bottom of the suspension registers a small localised increase in solids fraction (Buscall and White 1987; Davis and Russel 1989): the amount increases, and the size of the region where it occurs grow over time. Likewise in pressure filtration of a suspension, early on, most of the suspension is uniform, with spatial change in solids fraction being confined to a thin layer, albeit the thickness of that layer grows in extent over time (Landman and White 1994).

7 Conclusion

We have considered solutions for Richards equation using relative hydraulic conductivity and relative diffusivity from Van Genuchten’s soil functions (Van Genuchten 1980), and as well as the solution to two foam drainage equation variants. Specifically, we considered travelling wave behaviour which tends to set in at long times for typical infiltration conditions (Broadbridge and White 1988). In general, for Richards equation these solutions (Philip 1957a; Parlange 1971; Broadbridge and White 1988) have been obtained via Simpson’s rule, although analytic approximations are available for very low moisture content ($\varTheta \ll 1$), and systems near saturation ($\varTheta \rightarrow 1$) (Witelski 1997). Analytical travelling wave solutions are available for the channel- and node-dominated foam drainage cases for general $\varTheta$ as originally derived by Koehler et al. (2000); Verbist and Weaire (1994); Verbist et al. (1996), with further simplifications in the $\varTheta \ll 1$ and $\varTheta \rightarrow 1$ limits.

The foam drainage travelling waves have been compared with those for Richards equation. It was found that travelling wave velocities $\upsilon$ tend to be lower (see e.g. Fig. 5) in soils compared to foams. The only exceptions are for a soil that is already comparatively wet and then substantially more liquid is added to bring it close to full saturation. This behaviour follows from the shape of the soil relative hydraulic conductivity which is comparatively small for most moisture contents $\varTheta$, but which grows substantially near $\varTheta =1$.

Profiles of $\varTheta$ versus spatial coordinate ${\hat{\xi }}$ indicate that rise in moisture content $\varTheta$ with ${\hat{\xi }}$ is very abrupt (in the limit of small $\varTheta$ values) in the case of soils (Richards equation). This contrasts with the FDE in which $\varTheta$ rises much more gradually. As soon as a small $\varTheta$ is detected at a certain depth, moisture contents much larger than before quickly follow. On the other hand, if $\varTheta$ at a certain depth is found to be relatively large, the rate at which full saturation is reached in soil moving up in ${{\hat{\xi }}}$ is surprisingly slow: much slower than in the case of foams, as indeed follows from the $\varTheta \rightarrow 1$ analytical asymptotic formulae that we present. These observations suggest how rapidly an irrigation system might need to be switched off once a predetermined saturation is achieved at a certain depth. Moreover, when one evaluates the total amount of water that has entered the soil behind a travelling wave that has penetrated to a certain depth, it turns out that (in relative terms) more water is added in the case of soil than foam. In other words, the “missing” moisture needed to attain full saturation is less for soil than for foam: this difference between soil and foam results from the abrupt changes in $\varTheta$ with position in soils in the small $\varTheta$ limit, notwithstanding the gradual changes seen for higher $\varTheta$. More generally, this shows that the solutions to Richards equation based on travelling waves are useful for understanding water transport in soils.

Appendix: Asymptotic Form of Relative Hydraulic Conductivity

The original relative hydraulic conductivity (hereafter RHC) used in this work is given by Eq. (18). It is observed in Sect. 3.2 (see e.g. Fig. 2b  for Hygiene Sandstone) that the asymptotic function (20) obtained for $\varTheta \approx 1$ overestimates RHC as $m \rightarrow 1$. Retaining an additional term in the asymptotic expansion of Eq. (18), we deduce another expression that also estimates the RHC as $\varTheta \rightarrow 1$

$$\begin{aligned} \left. K_r(\varTheta ) \right| _{\varTheta \approx 1} \approx \left( 1-\frac{2}{m^m}(1-\varTheta )^m - \frac{1}{2}(1-\varTheta ) \right) . \end{aligned}$$

(44)

Hence,

$$\begin{aligned} \left. \varTheta - K_r(\varTheta ) \right| _{\varTheta \approx 1} \approx \frac{2}{m^m}(1-\varTheta )^m - \frac{1}{2}(1-\varTheta ). \end{aligned}$$

(45)

For $m = 0.5146$ (silt loam) or $m=0.6377$ (Guelph loam), the first term on the right hand side of (44) is dominant and the approximation considered in this “Appendix” is not required. However, for $m \rightarrow 1$, e.g. $m=0.9038$ considered here for Hygiene Sandstone, the term $\frac{1}{2}(1-\varTheta )$ could be around one quarter of the value of $2m^{-m}(1-\varTheta )^m$. Dropping the $\frac{1}{2}(1-\varTheta )$ term, therefore causes $\varTheta - K_r(\varTheta )$ to be overestimated by a factor of 4/3. If we look at the ratio of the true $\varTheta - K_r(\varTheta )$ to an approximation based on just the first term $2m^{-m}(1-\varTheta )^m$, then the ratio could be around $\frac{3}{4}$ which is what we see for $m = 0.9038$ in Fig. 9b, where there is a range of $\varTheta$ values for which the data stay close to 0.75. On the other hand, if we look at a ratio between the true $\varTheta - K_r(\varTheta )$ and the approximate formula (45), values much closer to unity result (see Fig. 10). Overestimating $\varTheta -K_r(\varTheta )$ causes us to underestimate $\mathrm {d}{\hat{\xi }}/\mathrm {d}\varTheta$ and hence ${\hat{\xi }}$. Including an extra term in the estimate of $\varTheta -K_r(\varTheta )$ would improve the estimate of $\mathrm {d}{\hat{\xi }}/\mathrm {d}\varTheta$ via (30) but would also preclude us from obtaining an analytical formula for ${\hat{\xi }}$ vs $\varTheta$ close to $\varTheta \approx 1$, so we do not pursue it in this paper.

Fig. 10

Comparison of the ratio between the true $\varTheta - K_r(\varTheta )$ and the approximation from Eq. (20), namely $\varTheta - K_r(\varTheta ) \approx 1-K_r(\varTheta ) \approx 2m^{-m}(1-\varTheta )^m$ (dashed lines) and also (solid line) the ratio of the true $\varTheta - K_r(\varTheta )$ to the approximate $\varTheta - K_r(\varTheta )$ in (45)