Capillary rise in partially saturated rigid porous media

Abstract

We explore predictions of two models of one-dimensional capillary rise in rigid and partially saturated porous media. One is an existing one from the literature and the second is a free-boundary model based on Richards’ equation with two moving boundaries of the evolving partially saturated region. Both models involve the specification of saturation-dependent functions for local capillary pressure and permeability and connect to classical models for saturated porous media. Existing capillary-rise experiments show two notable regimes: (i) an early-time regime typically well-described by classical capillary-rise theory in a fully saturated porous media, and (ii) a long-time regime that has anomalous dynamics in which the capillary-rise height may scale with a non-classical power law in time or have more complicated dynamics. We demonstrate that the predictions of both models compare well with experimental capillary-rise data over early- and long-time regimes gathered from three independent studies in the literature. The model predictions also shed light on recent scaling laws that relate the capillary pressure and permeability of the partially saturated media to the capillary-rise height. We use these models to probe computationally observed permeability relationships to capillary-rise height. We demonstrate that a recently proposed permeability scaling for the anomalous capillary-rise regime is indeed realized and is particularly apparent in the lower portion of the partially saturated media. For our free-boundary model we also compute capillary pressure measures and show that these reveal the linear relation between the capillary pressure and capillary-rise height expected for a capillarity–gravity balance in the upper portion of the partially saturated porous media.

1. Introduction

A fundamental problem in porous media is capillary-driven flow into a dry porous material. Commonly observed in papers and soil, capillary phenomena are important in many industrial applications, as well as other areas such as biomechanics and environmental science. Experiments by Bell & Cameron (1906) well over a century ago revealed an empirical relationship of the form $h^n \sim t$ where $h$ is the distance the liquid has travelled through the porous media in time $t$ with the exponent $n$ ‘above 2 in most cases’. These experiments focused on early times and/or horizontal flows in which gravitational forces were negligible compared with capillary forces. Just over a century ago, Washburn (1921) obtained what is now a classical result characterizing a quasisteady one-dimensional fluid flow in a capillary tube driven by capillary pressure at the meniscus. In Washburn's model, the capillary-rise height follows a square-root in time behaviour at early times before transitioning to an equilibrium height (Jurin's height) in which capillary forces at the meniscus and cylinder walls balance the opposing gravitational force. Under the assumption that a porous material can be approximated as a bundle of capillary tubes, a corresponding result was then extended to capillary dynamics in a porous material.

Ninety years after Bell & Cameron's work, Delker, Pengra & Wong (1996) examined the problem of capillary rise of water in a porous media composed of glass beads. Their measurements of capillary rise spanned times scales of seconds to more than a day. The early capillary-rise dynamics were well described by Washburn's model. At later times, however, there was a clear deviation from the classical dynamics as the capillary-rise height continued to increase in what they called anomalous dynamics that appeared in some cases to be represented by a power law, $t^\beta$, and in other cases showed more complicated behaviour. Delker et al. (1996) posed a model that included interface pinning and an interface velocity that depended algebraically on the amount by which the driving force for motion exceeded a threshold value. Lago & Araujo (2001) followed this work with a similar set of experiments on capillary rise of water into porous media composed of glass beads and another set using Berea sandstone. For glass beads, their experiments showed a transition from the classical $t^{1/2}$ dynamics at early times to a different power law at long times, $t^\beta$, where $1/20 < \beta < 1/4$. Shikhmurzaev & Sprittles (2012) extended these ideas and argued that the introduction of a dynamic contact angle model incorporating different modes of meniscus motion related to dynamic wetting, a contact line pinning and interface depinning allows a good fit with these data. Lunowa et al. (2022) have also recently explored the influence of a dynamic contact angle, slip and inertia on capillary rise of fluids in cylindrical containers. Configurations with variable cross-section porous materials and/or radial (source type) flows (e.g. Reyssat et al. 2008; Xiao, Stone & Attinger 2012; Perez-Cruz, Stiharu & Dominguez-Gonzalez 2017) can also display departures from the classical $t^{1/2}$ capillary imbibition dynamics directly attributable to the macroscale geometry.

A valuable approach to develop predictive tools that are able to capture in a single framework both classical and anomalous dynamic regimes in capillary-rise problems is based on the general ideas of mixture theory and continuum descriptions of the porous media. Richards’ equation (Richards 1931), which accounts for partial saturation of the porous media, is based on conservation of mass coupled to a Darcy-type flow and requires the specification of closure relationships for the dependence of capillary pressure and hydraulic conductivity on the variable level of saturation through the media (e.g. see Pillai & Hooman 2013). With such closure relationships specified, this approach leads to very powerful models from which detailed spatial and temporal dynamics can be predicted computationally and/or analytically. Lockington & Parlange (2004) used Richards’ equation as the starting point for their model along with closure relationships in the form of exponential functions of capillary pressure for hydraulic conductivity and saturation. In their subsequent analysis based on a travelling-wave solution approximation they obtained closed-form expressions for capillary-rise height and time linked parametrically via a volume flux variable. They demonstrated that their model captured the early-time classical dynamics and that parameters related to their closure relations allowed for a range of anomalous dynamics in the long-time regime. Along with the predictions of our own related model, we shall revisit the Lockington & Parlange (2004) model and explore its predictions more deeply in the context of the Lago & Araujo (2001) and Delker et al. (1996) data and more recent developments.

Owing to the importance of porous media flows in a wide range of applications (perhaps most notably in soil science) as well as the great utility of theory like Richards’ equation, there exist a wealth of studies, measurements, data and models on the relationships between capillary pressure, hydraulic conductivity or permeability and media saturation. These are often expressed as soil water retention curves and relative hydraulic conductivity relationships and include well-known works by Leverett (1941), Brooks & Corey (1964), Mualem (1976) and van Genuchten (1980), as well as related ones adopted by Lockington & Parlange (2004) (but see also Bear (1972) and Hornung (1997)). More recent investigations include Kuang et al. (2020) who discuss a modification to the classical van Genuchten model for relative hydraulic conductivity and the work by Soldi, Guarracino & Jougnot (2017) which has a specific focus on hysteretic features. Another recent study by Johnson, Zyvoloski & Stauffer (2019) incorporates an additional porosity dependence of these functions to be used in problems in which the porosity evolves in time (e.g. when the porous matrix itself freezes, melts/dissolves or otherwise deforms). In general the pressure versus media saturation depends on whether imbibition or draining (e.g. irrigation and fluid recovery in soil science terminology) occurs. This is known as capillary hysteresis. Our focus is exclusively on imbibition (no hysteresis) and the key features of the capillary pressure versus saturation curves in this case are (i) a sharp rise in the capillary pressure at the low end of saturation (at a residual saturation which may or may not be zero), and (ii) below a fixed pressure the porous media remains completely saturated. We choose our closure relations from Brooks & Corey (1964), who document such relationships for a wide variety of porous materials. The dependence of capillary pressure on saturation is indeed an idealized version of a more general relationship that accounts for non-equilibrium pore-scale dynamics. In particular, our model posed below reflects an assumption that the time scales for equilibration of the local interfacial and volume dynamics are rapid compared with the time scales associated with capillary rise (e.g. see Gray & Miller 2014, § 11.7, pp. 459–460).

This capillary-rise problem, along with related investigations of capillary rise in soft porous materials has seen much attention recently. Mirzajanzadeh, Deshpande & Fleck (2019), for example, investigated the capillary-rise dynamics in cellulose sponges and developed a model in which the fluid motion was diffusively controlled beyond the Jurin height. They justified this model by pointing out an insensitivity of the dynamics to the direction of gravity and that the water continued to rise after the water reservoir had been removed. Mirzajanzadeh, Deshpande & Fleck (2020) examined the two-dimensional deformation that accompanies capillary rise into an initially compressed sponge. Other investigations include the work of Kim, Moon & Kim (2016) on hemiwicking, as well as on capillary rise in cellulose sponges by Kim, Ha & Kim (2017) and Ha et al. (2018). These and related studies on one-dimensional capillary rise in soft porous materials have been reviewed recently by Ha & Kim (2020).

In their review, Ha & Kim (2020) outlined scaling arguments that capture the fundamental physics associated with various observed power law, $t^{\beta }$, dynamics. As we are interested in exploring these relationships in our computational models, here we outline two of their arguments that capture the essence of the $t^{1/2}$ and $t^{1/4}$ power laws for partially saturated capillary rise in a porous media (see also Kim et al. 2017). The essential ingredients in their arguments are Darcy's equation and careful interpretation of quantities such as capillary pressure and permeability in partially saturated porous media. Darcy's law expresses that the rate of change of capillary-rise height $h(t)$ times liquid fraction is related to the flux so that

(1.1)\begin{equation} \frac{{\rm d}h}{{\rm d}t} \sim{-} \frac{k}{\mu} \left( \frac{\partial p}{\partial z} + \rho g \right), \end{equation}

where $k$ is the permeability of the material, $\mu$ is the dynamic viscosity, $p$ is the pressure, $\rho$ is the fluid density, $g$ is the gravitational acceleration and $z$ is the vertical coordinate.

Early-time dynamics: Kim et al. (2017) and Ha & Kim (2020) argued that at early times the fluid fills the void space on the macroscale and the capillary pressure is inversely proportional to the macroscale pore radius, $r$, so that $p_c \sim \sigma / r$, where $\sigma$ is the surface tension. Here, the permeability is dominated by the macroscale features and has $k \sim r^2$. The pressure gradient then scales with $-p_c/h$ and, for sufficiently small $h$, dominates gravitational effects. Then,

(1.2)\begin{equation} \frac{{\rm d}h}{{\rm d}t} \approx \frac{k}{\mu} \frac{p_c}{h} \sim \frac{\sigma r}{\mu h}, \end{equation}

which leads to the scaling $h \sim (\sigma r t/\mu )^{1/2}$ consistent with the classical Washburn $t^{1/2}$ dynamics.

Long-time dynamics: Kim et al. (2017) and Ha & Kim (2020) argued that at later times, as the fluid rises in the sponge, it fails to fill the void space on the macroscale. Darcy's equation (1.1) still applies, but gravity now plays a critical role, the driving capillary pressure reflects imbibition processes in microscale pores, and the permeability scale must reflect the partially saturated nature of the flow through the macropores. In particular, beyond Jurin's height, there is a balance between capillarity and gravity in partially filled macropores in which $\sigma /r_g \sim \rho g h$ where $r_g$ is a representative length scale for the meniscus in the macropore (Kim et al. (2017) refer to this length, $\lambda$ in their notation, as the radius of the meniscus curvature in the partially filled macropores – see their figure 4). This length scale also sets the scale for the permeability so that $k \sim r_g^2$. These two relationships also indicate that $k \sim (\sigma /(\rho g h))^2$ so the permeability, apparently, is inversely proportional to the capillary-rise height squared. That the capillary pressure and permeability, both space- and time-dependent quantities in the partially saturated porous media, have a particular scaling with respect to capillary-rise height, a purely time-dependent quantity, is a subtle point that we shall explore more deeply in our present work. Continuing with Kim, Ha and coworker's arguments, further imbibition is driven by the capillary pressure operating on the microscale pore size so that $p_c \sim \sigma /r_m$ where $r_m$ is a microscale pore where $r_m \ll r_g \ll r$. These scalings together give

(1.3)\begin{equation} \frac{{\rm d}h}{{\rm d}t} \approx \frac{k}{\mu} \frac{p_c}{h} \sim \frac{r_g^2}{\mu} \frac{\sigma}{r_m h} \sim \frac{\sigma^3}{\mu r_m \rho^2 g^2}\frac{1}{h^3}, \end{equation}

from which the $t^{1/4}$ scaling follows:

(1.4)\begin{equation} h \sim \left( \frac{\sigma^3 t}{\mu r_m \rho^2 g^2} \right)^{1/4}. \end{equation}

Recognizing these three disparate length scales in the scaling argument above indicates that there are also three disparate scales for capillary pressure in which the capillary pressure on the micropore-scale is much larger than that associated with the fully or partially filled macropores. Ha et al. (2018) showed that a related scaling argument for cases in which swelling occurs has a velocity-dependent micropore length scale $r_m$, and reveals a $t^{1/5}$ scaling that appears consistent with experimental observations on capillary rise of water in cellulose sponges (Siddique, Anderson & Bondarev 2009; Kim et al. 2017; Ha et al. 2018; Ha & Kim 2020).

In the present work, we focus on two models – the one by Lockington & Parlange (2004) and one we derive here using a slightly more general free-boundary formulation – both of which contain essentially the same physics outlined in the scaling arguments above. Our first goal will be to explore the predictions of these models in direct comparison with the capillary rise data of Lago & Araujo (2001) and Delker et al. (1996). As such, we shall focus only on capillary imbibition and not address capillary-pressure hysteresis (e.g. see Mitra et al. 2020). In addition to the capillary-rise experimental data, we also use these model predictions to directly explore the clever but subtle scaling arguments proposed by Kim et al. (2017) and Ha & Kim (2020) that connect the permeability to capillary rise dynamics. Since both of these models have capillary pressure and permeability that vary in both space and time through the partially saturated porous media, we introduce different time-dependent measures (effectively integrals over space) for the capillary pressure and permeability that, a posteriori, can be compared with the capillary-rise dynamics. We believe that this is the first time such a comparison with these scaling predictions, specifically with respect to the capillary pressure and permeability of the partially saturated media, has been made. In making this comparison we shed light on deeper interpretations of these scaling laws and their connection to the observed anomalous capillary-rise dynamics.

Our paper is organized as follows. In the next section we derive our model. Included here are the specific capillary pressure and permeability relationships from Brooks & Corey (1964) that we employ. Section 3 provides a review of the Lockington & Parlange (2004) model and a comparison of their related capillary pressure and permeability relationships to those of Brooks & Corey (1964). In § 4 we compare both the Lockington–Parlange model and our model predictions to the capillary-rise data of Lago & Araujo (2001) and Delker et al. (1996). This section includes parameter estimations for both of these models as well as connections to the classical Washburn model predictions. In § 5 we turn to the permeability and capillary pressure measures and compare these predictions to recently proposed scaling relations associated with the anomalous capillary-rise dynamics. This section also includes a brief comparison with one of the capillary-rise experiments reported in Kim et al. (2017). Section 6 contains the conclusions. A few additional technical details are given in Appendix A for a similarity solution and in Appendix B for an analysis of permeability measures.

2. Mathematical modelling

The model formulation we use here is based on continuum mechanics and mixture theory where at each point in space and time in the porous material one can define phase volume fractions, along with field variables such as velocity and pressure. The model invokes averaging in the sense that a single point in space represents a small elemental volume in which all three phases – solid, liquid and gas – may be present in some proportion. This approach has been used extensively, especially in the context of flows in deformable porous materials, and builds on the pioneering work of Biot (1941a,b,c) with applications in biomechanics (e.g. Holmes 1983,  1984,  1985; Holmes & Mow 1990; Lai, Hou & Mow 1991; Barry & Aldis 1993,  1997) infiltration (e.g. Preziosi, Joseph & Beavers 1996; Sommer & Mortensen 1996; Michaud, Sommer & Mortensen 1999; Ambrosi & Preziosi 2000; Billi & Farina 2000) printing (e.g. Chen & Scriven 1990; Fitt et al. 2002; Anderson 2005) and capillary-driven imbibition in various contexts (e.g. Siddique et al. 2009; Siddique & Anderson 2011; Anderson & Siddique 2013; Mirnyy et al. 2013).

The configuration of interest is one-dimensional capillary rise in a rigid porous material as outlined below. This is a common experimental configuration with available data and theory with which we can compare our work (e.g. Delker et al. 1996; Lago & Araujo 2001; Kim et al. 2017). In the configuration of interest shown in figure 1, a reservoir supplies fluid at the bottom of the porous material, $z=0$ at time $t=0$, one interface position, $z=h_{\ell }^S(t)$, marks the top of the fully saturated region which occupies $0 < z < h_{\ell }^S(t)$, and a second interface position, $h_{\ell }(t)$, defines the top of the wet region of the porous material, so that the partially saturated porous region occupies $h_{\ell }^S(t) < z < h_{\ell }(t)$. This formulation leads to Richards’ equation as a description for the flow in the partially saturated region (Richards 1931). Although this type of model has been studied extensively in other configurations and contexts (e.g. Witelski 2003; Lockington & Parlange 2004; Pillai & Hooman 2013; Tafreshi & Bucher 2013; Perez-Cruz et al. 2017) we briefly outline its derivation and the coupling to the saturated region here.

Figure 1. This figure shows the configuration under consideration for capillary rise in partially saturated porous media.

In the partially saturated region $h_{\ell }^S(t) < z < h_{\ell }(t)$ the mass and momentum balances for a one-dimensional configuration are

(2.1)$$\begin{gather} \frac{\partial \phi_{\ell}}{\partial t} + \frac{\partial}{\partial z} \left( w_{\ell} \phi_{\ell} \right) = 0, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \phi_{g}}{\partial t} + \frac{\partial}{\partial z} \left( w_{g} \phi_{g} \right) = 0, \end{gather}$$

(2.3)$$\begin{gather}- \phi_{\ell} \frac{\partial p_{\ell}}{\partial z} - \rho_{\ell}^T \phi_{\ell} g - K_{s \ell} w_{\ell} + K_{g\ell} ( w_g - w_{\ell} ) = 0, \end{gather}$$

(2.4)$$\begin{gather}- \phi_{g} \frac{\partial p_g}{\partial z} - \rho_{g}^T \phi_{g} g - K_{s g} w_g + K_{g\ell} ( w_{\ell} - w_g ) = 0, \end{gather}$$

where $\phi _{\ell }$ and $\phi _g$ are the liquid and gas volume fractions, $w_{\ell }$ and $w_g$ are the vertical components of the liquid and gas velocities, $p_{\ell }$ and $p_g$ are the liquid and gas pressures and $\rho _{\ell }^T$ and $\rho _g^T$ are assumed constant bulk liquid and gas densities. We assume a known, fixed, solid volume fraction $\phi _s$ so that the condition $\phi _{\ell } + \phi _g + \phi _s =1$ effectively establishes $\phi _g$ in terms of $\phi _{\ell }$ and $\phi _s$. The friction coefficients, $K_{ij}$, are associated with relative motion between phase $i$ and phase $j$. We shall also assume that the gas phase density is negligible and from here onward set $\rho _{g}^T =0$. Further, we assume the gas phase viscosity is sufficiently small to generate negligible friction with either liquid or solid phases – that is, the $K_{g\ell }$ and $K_{s g}$ terms are assumed negligible. Therefore, the only relative velocity of consequence is that between the liquid and (motionless) solid phase. This leads to a hydrostatic pressure for the gas, assuming $\phi _g \neq 0$, so that (2.4) is replaced with $\partial p_g /\partial z =0$ or $p_g = p_A$ where $p_A$ is atmospheric pressure at $z=h_{\ell }$. Then, the liquid momentum equation (2.3) is

(2.5)\begin{equation} K_{s\ell}(\phi_{\ell}) w_{\ell} = \phi_{\ell} \frac{\partial p_c}{\partial z} - \phi_{\ell} \rho_{\ell}^T g, \end{equation}

where we have defined the capillary pressure $p_c$ as

(2.6)\begin{equation} p_{\ell}- p_g ={-}p_c(\phi_{\ell}). \end{equation}

Note that both $K_{s\ell }(\phi _{\ell })$ and capillary pressure $p_c(\phi _{\ell })$ depend on the liquid fraction $\phi _{\ell }$. Equation (2.6) expresses that because the pore space is partially saturated the pore-scale liquid pressure and the gas pressure are different and in general vary with saturation (e.g. Leverett 1941; Brooks & Corey 1964; Bear 1972; Mualem 1976; van Genuchten 1980; Hornung 1997; Gray & Miller 2014; Soldi et al. 2017; Kuang et al. 2020). We outline more details for the function $p_c$ below.

Inserting $w_{\ell }$ into the mass-conservation equation (2.1) gives an evolution equation for $\phi _{\ell }$,

(2.7)\begin{equation} \frac{\partial \phi_{\ell}}{\partial t} ={-} \frac{\partial}{\partial z} \left[ \frac{\phi_{\ell}^2}{K_{s\ell}(\phi_{\ell})} \left( \frac{\partial p_c}{\partial z} - \rho_{\ell}^T g\right) \right], \end{equation}

on $h_{\ell }^S(t) < z < h_{\ell }(t)$. Equation (2.2) determines the evolution of $w_g$ but is otherwise decoupled from the dynamics. Identifying $\phi _{\ell }^2/K_{s\ell }(\phi _{\ell }) = (k_{sat}/ \mu ) k_{r\ell }(\phi _{\ell })$ where $k_{sat}$ is the permeability at full saturation, $\mu$ is the viscosity and $k_{r\ell }(\phi _{\ell })$ is the relative permeability of the (wetting) liquid phase (whose form will be specified below), and keeping in mind that $p_c$ is a function of $\phi _{\ell }$ we can identify this as Richards’ equation,

(2.8)\begin{equation} \frac{\partial \phi_{\ell}}{\partial t} ={-} \frac{\partial}{\partial z} \left[ \frac{k_{sat} k_{r\ell}(\phi_{\ell}) }{\mu} \left( \frac{\partial p_c }{\partial z} - \rho_{\ell}^T g \right) \right]. \end{equation}

One boundary condition at $z=h_{\ell }(t)$ prescribes that the capillary pressure there is known; $p_c = p_c^{1}$. In view of (2.6) note that a non-zero value of $p_c^1$ corresponds to a jump in liquid and gas pressures at this boundary owing to capillary effects. In our formulation the capillary pressure in the partially saturated media is assumed to be a monotonic function of liquid fraction (or, in (2.23a,b) saturation). Therefore, we can express this condition on the capillary pressure as a boundary condition that fixes the liquid fraction at $z=h_{\ell }(t)$; that is,

(2.9)\begin{equation} \phi_{\ell} = \phi_{\ell}^1 \quad\text{at }z=h_{\ell}(t), \end{equation}

where $p_c(\phi _{\ell }^1) = p_c^1$. A second boundary condition is $w_{\ell }(z=h_{\ell }(t)) = {{\rm d}h_{\ell }}/{{\rm d}t}$ which gives the evolution of $h_{\ell }(t)$,

(2.10)\begin{equation} \frac{{\rm d} h_{\ell}}{{\rm d}t} = \frac{k_{sat} k_{r\ell}(\phi_{\ell}^1)}{\mu \phi_{\ell}^1} \left. \left( \frac{\partial p_c }{\partial z} - \rho_{\ell}^T g \right) \right|_{z=h_{\ell}}. \end{equation}

Finally, at the bottom of the partially saturated region, $z=h^S_{\ell }(t)$, we assume that the capillary pressure is that associated with full saturation, $p_c(\phi _{\ell } = 1-\phi _s)$. Again, in terms of a boundary condition on the liquid fraction we express this as $\phi _{\ell } = 1- \phi _s$ at $z=h^S_{\ell }(t)$.

In the fully saturated region, where the porosity $\phi _{\ell } = 1 - \phi _s$ is constant, we have

(2.11)$$\begin{gather} w_{\ell} ={-} \frac{k_{sat}}{\mu ( 1 - \phi_s ) } \left( \frac{\partial p_{\ell} }{\partial z} + \rho_{\ell}^T g \right), \end{gather}$$

(2.12)$$\begin{gather}\frac{\partial w_{\ell}}{\partial z} = 0. \end{gather}$$

Boundary conditions correspond to

(2.13)$$\begin{gather} p_{\ell} = p_A + \rho_{\ell}^T g z_R, \quad \text{at }z=0, \end{gather}$$

(2.14)$$\begin{gather}p_{\ell} = p_A - p_c(\phi_{\ell}=1-\phi_s),\quad \text{at }z=h_{\ell}^S(t), \end{gather}$$

(2.15)$$\begin{gather}w_{\ell}(z=h_{\ell}^S(t)) = \frac{{\rm d} h_{\ell}^S(t)}{{\rm d}t},\quad \text{at }z=h_{\ell}^S(t), \end{gather}$$

where (2.13) reflects that the fluid reservoir is maintained at height $z=z_R$ above the bottom of the porous region and (2.14) reflects that fully saturated conditions apply at $z=h_{\ell }^S$. The third condition (2.15) prescribes that the interface $z=h_{\ell }^S(t)$ moves with the fully saturated fluid velocity there. In this fully saturated region, we can solve for pressure by noting that (2.11) and (2.12) imply $\partial ^2 p_{\ell }/\partial z^2 =0$; that is, $p_{\ell }$ is a linear function of $z$. Applying boundary conditions (2.13) and (2.14) on $p_{\ell }$ then gives

(2.16)\begin{equation} p_{\ell} = p_A + \rho_{\ell}^T g z_R + \frac{z}{h_{\ell}^S(t)} \left( - p_c(\phi_{\ell}=1-\phi_s) - \rho_{\ell}^T g z_R \right). \end{equation}

Substituting this result for pressure into (2.11) and inserting the resulting expression for $w_{\ell }$ into (2.15) gives an ordinary differential equation for $h_{\ell }^S(t)$, the dimensionless version of which is listed in § 2.2, (2.21).

2.1. Equilibrium conditions

The model given above has a long-time, equilibrium solution, with interface positions

(2.17)$$\begin{gather} h_{\ell}(t \rightarrow \infty) = h_{\ell}^{S}(t \rightarrow \infty) + \frac{p_c(\phi_{\ell} = \phi_{\ell}^1) - p_c(\phi_{\ell}=1-\phi_s)}{\rho_{\ell}^T g} = z_R + \frac{p_c^1}{\rho_{\ell}^T g}, \end{gather}$$

(2.18)$$\begin{gather}h_{\ell}^{S}(t \rightarrow \infty) = z_R + \frac{p_c(\phi_{\ell}=1-\phi_s)}{\rho_{\ell}^T g} = z_R + h_e, \end{gather}$$

where we have introduced the hydraulic equilibrium height defined by $h_e = p_c(\phi _{\ell }=1-\phi _s)/(\rho _{\ell }^T g)$ which marks the height of the fully saturated region in equilibrium (see also Lago & Araujo 2001). Note that (2.17) follows from setting $w_{\ell }=0$ in (2.5), cancelling a common factor $\phi _{\ell }$ and integrating from $z=h_{\ell }^{S}$ to $z=h_{\ell }$. Similarly, (2.18) follows from setting $w_{\ell }=0$ in (2.11), integrating from $z=0$ to $z=h_{\ell }^{S}$ and applying boundary conditions (2.13) and (2.14). Note that if $p_c(\phi _{\ell } = \phi _{\ell }^1) = p_c^1$ is finite then $h_{\ell }(t \rightarrow \infty )$ is finite (this is the case in our model as explained further below). In the mathematical limit $p_c^1 \rightarrow \infty$, then $h_{\ell }(t \rightarrow \infty ) \rightarrow \infty$ (this is the limiting case associated with $\phi _{\ell }^1 \rightarrow 0$ for the capillary-pressure in the Lockington–Parlange model as explained further below).

2.2. Non-dimensionalization

We summarize our model in terms of dimensionless variables, denoted by bars, and introduce the saturation variable, $S = \phi _{\ell }/(\phi _g + \phi _{\ell })$, representing the proportion of the pore space occupied by liquid. Let $z = H_0 \bar {z}$, $h_{\ell } = H_0 \bar {h}_{\ell }$, $h_{\ell }^S = H_0 \bar {h}_{\ell }^S$, $t = T_0 \bar {t}$ and $p_c = P_0 \bar {p}_c$ where $H_0$ is a length scale, $T_0$ is a time scale, and $P_0$ is a pressure scale. Inserting these into the governing equations and choosing $H_0 = z_R + h_e$, $T_0 = ((1-\phi _s)\mu H_0^2)/(k_0 P_0)$ and $P_0 = p_c(S=1) = \rho _{\ell }^T g h_e$ (i.e. capillary pressure value at full saturation) gives

(2.19)$$\begin{gather} \frac{\partial S}{\partial \bar{t}} ={-} \frac{\partial}{\partial \bar{z}} \left[ \bar{k}_{r\ell} \left( \frac{\partial \bar{p}_c}{\partial \bar{z} } - \frac{1}{1-\bar{z}_R} \right) \right], \quad \text{on } \bar{h}_{\ell}^S < \bar{z} < \bar{h}_{\ell}, \end{gather}$$

(2.20)$$\begin{gather}\frac{{\rm d}\bar{h}_{\ell}}{{\rm d}\bar{t}} = \left. \left[ \frac{\bar{k}_{r\ell} }{S} \left( \frac{\partial \bar{p}_c}{\partial \bar{z}} - \frac{1}{1-\bar{z}_R} \right) \right] \right|_{\bar{z}=\bar{h}_{\ell}}, \end{gather}$$

(2.21)$$\begin{gather}\frac{{\rm d} \bar{h}_{\ell}^S}{{\rm d}\bar{t}} = \frac{1}{1-\bar{z}_R} \left( \frac{1}{\bar{h}_{\ell}^S} - 1 \right), \end{gather}$$

where $\bar {z}_R = z_R/(z_R + h_e)$. Since typical capillary-rise experiments have $z_R \ge 0$ and $h_e >0$ we expect $0 \le \bar {z}_R < 1$. These equations are subject to boundary conditions $S=1$ at $\bar {z} = \bar {h}_{\ell }^S$ and $S=S_1$ at $\bar {z} = \bar {h}_{\ell }$, where $S_1 = \phi _{\ell }^1/(1-\phi _s)$. The initial conditions for these equations are $\bar {h}_{\ell }^S( \bar {t}=0) = 0$ and $\bar {h}_{\ell }( \bar {t}=0) = 0$. We establish the initial evolution for $\bar {t} \in [0,\bar {t}^*]$ (where $\bar {t}^*$ is some small but non-zero value of time) via a similarity solution which describes the formation of the partially saturated medium. Then, this similarity solution is used as a starting saturation profile $S(\bar {z},\bar {t}^*)$ for $\bar {z} \in [\bar {h}_{\ell }^S(\bar {t}^*),\bar {h}_{\ell }(\bar {t}^*)]$, where $\bar {h}_{\ell }(\bar {t}^*) > \bar {h}_{\ell }^S(\bar {t}^*) > 0$ to continue integration for $\bar {t} > \bar {t}^*$. Details are given below after we specify the capillary pressure and permeability functions, $\bar {p}_c(S)$ and $\bar {k}_{r\ell }(S)$. The dynamics of $\bar {h}_{\ell }^S$ are the classical Washburn dynamics whereas the dynamics represented by $\bar {h}_{\ell }$ account for partial saturation.

To connect these results to others of Lago & Araujo (2001) and Lockington & Parlange (2004) we note that the time scale $T_0$ can be expressed as

(2.22)\begin{equation} T_0 = \frac{(1-\phi_s) \mu H_0^2}{k_{sat} P_0} = \frac{z_R + h_e}{v_g} \frac{1}{1-\bar{z}_R}, \end{equation}

where we introduce the velocity scale $v_g = K_{sat} / (1-\phi _s)$ where $K_{sat} = k_{sat} \rho _{\ell }^T g/\mu$ is the saturated hydraulic conductivity (see also Lago & Araujo 2001, equation (3); Lockington & Parlange 2004).

2.3. Brooks–Corey expressions for capillary pressure and permeability

We use capillary pressure and relative permeability expressed as functions of saturation based on Brooks & Corey (1964). Specifically, these are characterized by

(2.23a,b)\begin{equation} \bar{p}_c(S) = \left( \frac{1 - S_r}{S - S_r} \right)^{1/\lambda_{BC}}, \quad \bar{k}_{r\ell}(S) = \left( \frac{S - S_r}{1 - S_r} \right)^{3+ 2/\lambda_{BC}}, \end{equation}

where $S_r$ is a residual saturation and $\lambda _{BC}$ is a constant. We have chosen to use these forms of capillary pressure and permeability as these have been established for a wide range of porous media including the case of glass beads used in the experiments of Lago & Araujo (2001) and Delker et al. (1996). We further demonstrate below that these forms match closely to those used in the Lockington–Parlange model and this allows us to directly compare the predictions of these two models.

If we define

(2.24)\begin{equation} \bar{M}(S) ={-} \frac{\partial \bar{p}_c}{\partial S} = \frac{1}{\lambda_{BC}} \frac{1}{1-S_r} \left( \frac{1-S_r}{S-S_r} \right)^{1 + 1/\lambda_{BC}}, \end{equation}

the governing equations can be written as

(2.25)$$\begin{gather} \frac{\partial S}{\partial \bar{t}} = \frac{\partial}{\partial \bar{z}} \left[ \bar{k}_{r\ell}(S) \left( \bar{M}(S) \frac{\partial S}{\partial \bar{z} } + \frac{1}{1-\bar{z}_R} \right) \right], \quad \text{on }\bar{h}_{\ell}^S < \bar{z} < \bar{h}_{\ell}, \end{gather}$$

(2.26)$$\begin{gather}\frac{{\rm d}\bar{h}_{\ell}}{{\rm d}t} ={-}\left. \left[ \frac{\bar{k}_{r\ell}(S) }{S} \left( \bar{M}(S) \frac{\partial S}{\partial \bar{z} } + \frac{1}{1-\bar{z}_R} \right) \right] \right|_{\bar{z}=\bar{h}_{\ell}}, \end{gather}$$

(2.27)$$\begin{gather}\frac{{\rm d} \bar{h}_{\ell}^S}{{\rm d}\bar{t}} = \frac{1}{1-\bar{z}_R} \left( \frac{1}{\bar{h}_{\ell}^S} - 1 \right), \end{gather}$$

subject to boundary conditions $S=1$ at $\bar {z} = \bar {h}_{\ell }^S$ and $S=S_1$ at $\bar {z} = \bar {h}_{\ell }$. The initial conditions use a similarity solution formulation shown in the next subsection. In all of our calculations we shall assume that $S_r=0$.

2.4. Similarity solution

When $\bar {h}_{\ell } \ll 1$ and $\bar {h}_{\ell }^S \ll 1$, the governing equations can be solved by introducing a similarity variable $\eta = \bar {z} /(2 \sqrt {\bar {t}})$ along with $\bar {h}_{\ell } = 2 \lambda _{\ell } \sqrt { \bar {t}}$ and $\bar {h}_{\ell }^S = 2 \lambda _{\ell }^S \sqrt { \bar {t}}$ and seeking a solution of the form $S(\bar {z},\bar {t}) = S(\eta )$ with constants $\lambda _{\ell }$ and $\lambda _{\ell }^S$ to be determined. In particular, in this limit the terms $1/(1-\bar {z}_R)$ in (2.25) and (2.26) along with the term $-1$ in parentheses on the right-hand side of (2.27) are absent (these are the original gravity terms). Then, we find that (2.25) can be rewritten in terms of the similarity variable $\eta$ as

(2.28)\begin{equation} {-}2 \eta \frac{{\rm d} S}{{\rm d}\eta} = \frac{{\rm d}}{{\rm d}\eta} \left[ \bar{k}_{r\ell}(S) \bar{M}(S) \frac{{\rm d} S}{{\rm d}\eta} \right], \end{equation}

for $\lambda _{\ell }^S < \eta < \lambda _{\ell }$ while (2.26) and (2.27) can be written as equations governing $\lambda _{\ell }$ and $\lambda _{\ell }^S$

(2.29)\begin{equation} 2 \lambda_{\ell} ={-} \left. \frac{\bar{k}_{r\ell}(S_1)}{S_1} \bar{M}(S_1) \frac{{\rm d} S}{{\rm d}\eta} \right|_{\eta = \lambda_{\ell}}, \quad \lambda_{\ell}^S = \sqrt{ \frac{1}{2 (1- \bar{z}_R) } }. \end{equation}

Note that in order for the partially saturated boundary to advance faster than the saturated boundary we need $\lambda _{\ell } > \lambda _{\ell }^S$. This condition is confirmed for the Brooks–Corey functions if $\lambda _{BC}>0$, as outlined in a solution of this problem in Appendix A.

3. Related Lockington–Parlange model

Before we present solutions for the free-boundary model outlined above, we give a brief summary of the related model by Lockington & Parlange (2004) and its connection to ours and the capillary pressure and permeability assumptions from Brooks & Corey (1964). The Lockington–Parlange model begins with Richards’ equation and extracts a parametric closed form prediction for the capillary-rise height based on a travelling-wave solution approximation along with specific choices for the functional dependence of capillary pressure and permeability (hydraulic conductivity) with respect to saturation.

In their notation, the Lockington–Parlange model (equations (34) and (35) in Lockington & Parlange (2004)) is given by

(3.1)$$\begin{gather} Z_f = \frac{1}{A} \ln \left( \frac{Q_0 + 1}{Q_0} \right) + \frac{1}{Q_0 + 1}, \end{gather}$$

(3.2)$$\begin{gather}T = \ln \left( \frac{Q_0 + 1}{Q_0} \right) - \frac{1}{Q_0 + 1} + \frac{1}{A \sqrt{1 + 3 \epsilon} } \left[ \frac{1}{Q_0} - \ln \left( \frac{Q_0 + 1}{Q_0} \right) \right], \end{gather}$$

where $Z_f$ and $T$ represent the dimensionless capillary-rise height and time, respectively, defined parametrically in terms of the dimensionless flux $Q_0$ at the bottom of the porous medium. These quantities are related to their dimensional counterparts, $z_f$, $t$, and $q_0$, by the expressions

(3.3a–c)\begin{equation} Z_f = \frac{z_f}{z_R + h_e},\quad T = \frac{t v_g}{z_R + h_e},\quad Q_0 = \frac{q_0}{K_{sat}}. \end{equation}

Two key dimensionless parameters, $A$ and $\epsilon$, appearing in this model are

(3.4a,b)\begin{equation} A = \alpha (z_R + h_e), \quad \epsilon = \frac{\beta - \alpha}{\alpha}. \end{equation}

The classic Washburn model is represented in the Lockington–Parlange model by the limit $A \rightarrow \infty$ (in which case the parameter $\epsilon$ drops out of the problem). The Lockington–Parlange model assumes $\epsilon \ll 1$. The parameters $\alpha$ and $\beta$ are constants in relationships between hydraulic conductivity, saturation and capillary pressure.

In an effort to avoid confusion about which capillary-rise height prediction we discuss, we retain the use of $Z_f$ and $T$ to indicate the predictions for the Lockington–Parlange model but note these have interpretations as $\bar {h}_{\ell }$ and $\bar {t}$ in our model (and similarly for dimensional variables). The parameter values such as $v_g$, $z_R$ and $h_e$ are the same across the two models although we note that Lockington & Parlange (2004) refer to $h_e$ as the air-entry pressure head ($h_a$ in their notation). The parameters $A$ and $\epsilon$ are unique to the Lockington–Parlange model. Lockington & Parlange (2004) introduce $\hat {h} \le 0$ (in their notation this is $h$) to denote their ‘matric potential head’ variable which can be interpreted as a measure of the negative of our capillary pressure via $\hat {h} = - p_c/(\rho _{\ell }^T g)$. With this notation in mind, their expressions for hydraulic conductivity and saturation are

(3.5)\begin{equation} K(\hat{h}) = K_{sat} \left\{ \begin{array}{@{}ll@{}} \exp \left[ \beta (\hat{h} + h_e) \right] & \hat{h} \le - h_e, \\ 1 & \hat{h} >{-} h_e, \end{array} \right. \end{equation}

and

(3.6)\begin{equation} S(\hat{h}) = \left\{ \begin{array}{@{}ll@{}} \exp \left[ (\beta - \alpha) (\hat{h} + h_e) \right] & \hat{h} \le - h_e, \\ 1 & \hat{h} >{-} h_e. \end{array} \right. \end{equation}

We note that for $\hat {h} \le - h_e$,

(3.7)\begin{equation} \frac{K}{K_{sat}} = \exp \left[ \beta (\hat{h} + h_e) \right] = \left( \exp \left[ (\beta- \alpha) (\hat{h} + h_e) \right] \right)^{\beta/(\beta- \alpha)} = S^{\beta/(\beta- \alpha)}, \end{equation}

indicating a power law relation between permeability and saturation. Also for $\hat {h} \le - h_e$ the hydraulic diffusivity ($K \,{\rm d}\hat {h}/{\rm d}S$) in Lockington & Parlange (2004) is

(3.8)\begin{equation} K \frac{{\rm d}\hat{h}}{{\rm d}S} = \frac{K_{sat}}{(\beta - \alpha)} S^{\alpha/(\beta- \alpha)} , \end{equation}

(see their equation (11) which is expressed in terms of volumetric water content), which is also a power law in terms of saturation.

To put the Lockington–Parlange model and ours on similar footing, below we connect these forms for permeability and hydraulic diffusivity to those using the Brooks & Corey (1964) expressions in (2.23a,b). Matching the Lockington & Parlange (2004) result for permeability in (3.7) with the Brooks & Corey (1964) result $\bar {k}_{r\ell }(S)$ in (2.23a,b), requires $\lambda _{BC} = \lambda _{BC}^{P}$ where $\lambda _{BC}^{P}$ is defined by the relation

(3.9)\begin{equation} 3 + 2 / \lambda_{BC}^{P} = \beta / ( \beta - \alpha) = \frac{1}{\epsilon} + 1. \end{equation}

That is, $\epsilon = \lambda _{BC}^{P}/(2 + 2 \lambda _{BC}^{P})$, or $\lambda _{BC}^{P} = 2 \epsilon /(1 - 2 \epsilon )$. On the other hand, note that the Brooks & Corey (1964) expression for hydraulic conductivity, $\bar {k}_{rl} \bar {M}$, is

(3.10)\begin{align} \bar{k}_{r\ell}(S) \bar{M}(S) & = \frac{1}{\lambda_{BC} (1-S_r)} \left( \frac{S - S_r}{1 - S_r} \right)^{3+ 2/\lambda_{BC}} \left( \frac{1 - S_r}{S - S_r} \right)^{1/\lambda_{BC} + 1} \nonumber\\ & = \frac{1}{\lambda_{BC} (1-S_r)} \left( \frac{S - S_r}{1 - S_r} \right)^{2+ 1/\lambda_{BC}}. \end{align}

To match this exponent with the Lockington & Parlange (2004) result (3.8) requires that $\lambda _{BC} = \lambda _{BC}^{D}$ where $\lambda _{BC}^{D}$ is defined by the relation

(3.11)\begin{equation} 2+ 1/\lambda_{BC}^{D} = \alpha/(\beta- \alpha) = \frac{1}{\epsilon}. \end{equation}

That is, $\epsilon = \lambda _{BC}^{D}/(1 + 2\lambda _{BC}^{D})$, or $\lambda _{BC}^{D} = \epsilon /(1 - 2\epsilon )$, which is a relationship that differs from that required to match the permeabilities (in particular $\lambda _{BC}^{P} = 2 \lambda _{BC}^{D}$).

So, the Brooks & Corey (1964) forms differ slightly from those of Lockington & Parlange (2004); technically either the permeability or hydraulic diffusivity, but not both simultaneously, can be matched. That said, the forms are effectively very similar as can be seen in a graphical comparison of the capillary pressure and permeability functions in figure 2. These comparisons capture the general features that (i) the capillary pressure builds from a unit reference pressure as the saturation drops away from one, and has a sharp increase at low saturation, and (ii) the permeability grows monotonically with saturation.

Figure 2. Panel (a) shows the dimensionless capillary pressure versus saturation curves for the Brooks–Corey model with $\lambda _{BC}=7.3$ and $S_r=0$ (black curve) and the Lockington–Parlange model with $\epsilon = \lambda _{BC}/(2 + 2 \lambda _{BC})$ and $\alpha = 14$ (red dashed curves and circles) and the Lockington–Parlange model with $\epsilon = \lambda _{BC}/(1 + 2\lambda _{BC})$ and $\alpha = 14$ (blue dotted curves and circles). Panel (b) shows the dimensionless permeability versus saturation curves for the Brooks–Corey model and Lockington–Parlange model for the same parameters and colour scheme as in (a). Note that the parameters for the red curves/circles were chosen to exactly match the Brooks–Corey and Lockington–Parlange permeability functions while the parameters for the blue curves/circles were chosen to match exponents of the Brooks–Corey and Lockington–Parlange hydraulic diffusivity functions. The key observation is that the Brooks–Corey functions match very closely with those employed in the Lockington–Parlange model.

Several other features of the Lockington–Parlange model are worth noting here. The early-time asymptotic behaviour of the Lockington–Parlange model in (3.1) and (3.2) corresponds to $Q_0 \gg 1$. We find that

(3.12)$$\begin{gather} Z_f = \left( 1 + \frac{1}{A} \right) \frac{1}{Q_0} - \left( 1 + \frac{1}{2A} \right) \frac{1}{Q_0^2} + O(Q_0^{{-}3}), \end{gather}$$

(3.13)$$\begin{gather}T = \frac{1}{2} \left( 1+\frac{1}{ A \sqrt{1 + 3\epsilon} } \right) \frac{1}{Q_0^2} + {O }(Q_0^{{-}3}). \end{gather}$$

This indicates that the early-time relationship between $Z_f$ and $T$ is

(3.14)\begin{equation} Z_f \sim \left( 1 + \frac{1}{A} \right) \sqrt{\frac{2 A \sqrt{1 + 3\epsilon} \; T }{1+ A \sqrt{1 + 3\epsilon} } } . \end{equation}

As noted by Lockington & Parlange (2004), with $A \rightarrow \infty$ the Washburn limit is recovered. However, outside of this limit ($1/A \neq 0$) both $A$ and $\epsilon$ influence the early-time dynamics and in general the details of the early-time Lockington & Parlange (2004) solution differ from those of the Washburn solution, although both follow the square-root in time scaling. The long-time asymptotic behaviour of the Lockington–Parlange model in (3.1) and (3.2) corresponds to $Q_0 \ll 1$. These two expressions can be expanded in this limit to reveal that

(3.15)$$\begin{gather} Z_f = \frac{1}{A} \ln \left( \frac{1}{Q_0} \right) + 1 + O(Q_0), \end{gather}$$

(3.16)$$\begin{gather}T = \frac{1}{A \sqrt{1 + 3\epsilon}} \frac{1}{Q_0} + \left(1 - \frac{1}{A \sqrt{1 + 3\epsilon}} \right) \ln \left( \frac{1}{Q_0} \right) - 1 + O(Q_0). \end{gather}$$

With $A \rightarrow \infty$ these give the expected Washburn limit $Z_f \rightarrow 1$. For finite $A$, these expressions show that $Z_{f}$ grows logarithmically in the long-time limit.

Our goals for revisiting the Lockington–Parlange model are three-fold. First, the Lockington–Parlange model is a relatively simple model involving closed form expressions and as such serves as a useful predictive tool and comparison for our free-boundary model. Second, to our knowledge, no quantitative comparison of the Lockington–Parlange model has been made to experimental data to the point where one can identify appropriate parameters $A$ and $\epsilon$ (Lockington & Parlange (2004) did demonstrate in their paper that their model recovers qualitatively the features of the Lago & Araujo (2001) experiments). Third, with reasonable comparison with both classical and anomalous dynamics in hand, solutions to the Lockington–Parlange model (as well as our model solutions) will be used to reveal detailed information about permeability dynamics. In particular, we shall directly explore the recently proposed permeability versus capillary-rise dynamic scaling proposed by Kim et al. (2017) and Ha & Kim (2020).

4. Comparison with experimental data

Lago & Araujo (2001) and Delker et al. (1996) report capillary rise of water in a vertical column of glass beads connected to a reservoir of fluid. We obtained the Lago & Araujo (2001) capillary rise versus time experimental data from their figure 10(a) using MATLAB's grabit.m software. We collected points from each of the curves corresponding to $150$–$180\ \mathrm {\mu }{\rm m}$, $180$–$212\ \mathrm {\mu }{\rm m}$, $212$–$250\ \mathrm {\mu }{\rm m}$ and $250$–$300\ \mathrm {\mu }{\rm m}$ glass bead sizes, and for simplicity refer to these below in tables and graphs by labels $150\ \mathrm {\mu }{\rm m}$, $180\ \mathrm {\mu }{\rm m}$, $212\ \mathrm {\mu }{\rm m}$ and $250\ \mathrm {\mu }{\rm m}$, respectively. These represent time and height data denoted by $(t_{exp}^i,h_{exp}^i)$ for $i=1, \ldots, N_{total}$ where $N_{total}$ is the number of data points. We used a similar procedure to obtain the Delker et al. (1996) capillary-rise data from their figure 2(a). Four data sets were obtained corresponding to their experiments for $180\ \mathrm{\mu }{\rm m}$, $253\ \mathrm {\mu }{\rm m}$, $359\ \mathrm {\mu }{\rm m}$ and $510\ \mathrm {\mu }{\rm m}$ bead diameters.

A technical but important detail about the Lago & Araujo (2001) capillary-rise data is that at their time zero the fluid height is at the reservoir height, $z_R$ (see figure 1). Figure 2(a) of Delker et al. (1996) shows capillary-rise height relative to $z=0$ with their starting time also corresponding to when the capillary-rise height is equal to $z_R$. To put the Delker et al. (1996) data in the form of the Lago & Araujo (2001) data, we shift the Delker et al. (1996) height values by $z_R$ (in their notation this is $Z - z_R$). In contrast, our model (and that of the Lockington–Parlange model) has time $t=0$ corresponding to $h_{\ell } = 0$. Therefore, in order to compare our predictions with the experiments we plot $h_{\ell }(t) - z_R$ versus $t - t_R$ where $t_R$ is defined to be the time for which $h_{\ell }(t_R) = z_R$ (in the Lockington–Parlange model we similarly plot $z_f - z_R$ versus $t - t_R$). Thus, ‘early time’ in the experimental plots means $t - t_R \ll 1$ in our model notation. A consequence of this is that the plotted capillary-rise height $h_{\ell }(t) - z_R$ for early time is linear in $t- t_R$ whenever $z_R \neq 0$. Further details of this are explained below.

In addition to comparing the Lago & Araujo (2001) and Delker et al. (1996) data with our model and the Lockington–Parlange model we briefly revisit the Washburn model. From a parameter estimation point of view, the Washburn model, as outlined below, involves three-dimensional parameters $h_e$, $z_R$ and $v_g$. For the Lockington–Parlange model five parameters appear: the dimensional ones $h_e$, $z_R$ and $v_g$, along with the dimensionless parameters $A$ and $\epsilon$. For our free-boundary model there are five as well: $h_e$, $z_R$, $v_g$ and the two dimensionless parameters $S_1$ and $\lambda _{BC}$.

Our optimization procedure to identify parameters involves minimizing the mean square error of capillary-rise height at the available experimental time points. We have used absolute errors in all cases described here but note that relative errors and/or other weighting schemes could be implemented. Our numerical scheme uses MATLAB's fmincon to search the parameter space. We have explored both interior-point and sqp algorithms with a variety of initial guesses to help find the best solution within our search space. In each case we limit our parameter search space with upper and lower bounds on each parameter. Some of these bounds are based on previous estimates in Lago & Araujo (2001) and/or Delker et al. (1996), especially with respect to the parameters $h_e$, $v_g$ and $z_R$, and are explained in detail below. The models of primary interest – the Lockington–Parlange model and our free-boundary model – have five-dimensional search spaces and even with bounds on these parameters we cannot guarantee that our reported solutions are global optimizers. Rather, we seek parameter values that are consistent with prior information and also reasonably fit the capillary-rise versus time data. With reasonable fits to the capillary-rise dynamics in both early- and long-time regimes we then further explore our model predictions on capillary pressure and permeability dynamics, which are not as easily accessible experimentally, and help shed further light on the capillary-rise phenomena.

The first step of our parameter estimation procedure involves fitting a linear function $t - t_R = P_1 (h_{\ell } - z_R)$ to early time data of Lago & Araujo (2001) and Delker et al. (1996) where $P_1$ has units of seconds per centimetre. That this early-time data has this linear structure, rather than the celebrated $t^{1/2}$ scaling, relates to the space and time shift of the data, and will be demonstrated in the context of the models below. We report numerical values of $P_1$ in table 1 for each Lago & Araujo (2001) and Delker et al. (1996) data set and note that these values of $P_1$ are independent of the theoretical model (Washburn, Lockington–Parlange, $\ldots$). The value of $N$ indicates the number of data points used in the fits. We see that the value of $P_1$ is relatively insensitive at small values of $N$ but necessarily drifts to larger values as more time points are included and the dynamics depart from the early-time behaviour. We have listed estimates and standard deviations in each case on the lines marked by $*$ in table 1.

Table 1. Fitted leading coefficient, $P_1$, from the model independent relation $t - t_R = P_1 (h_{\ell } - z_R)$ using capillary-rise data of Lago & Araujo (2001) and Delker et al. (1996). Here $N$ is the number of early time points used in the fit. Note that the case for Delker et al. (1996) 510 $\mathrm {\mu }$m beads only has six data points that appear to fall in the early-time regime. The mean and standard deviation for $P_1$ using the first four points ($N=6,7,8,9$ for Lago & Araujo (2001) and $N=3,4,5,6$ for Delker et al. (1996)) are listed on the lines marked by the $*$.

As shown below, in the Washburn model we can identify $P_1 = P_1^{Wash} = z_R/(v_g h_e)$ which fixes a relationship between the three parameters $h_e$, $v_g$ and $z_R$. In the Lockington–Parlange model $P_1 = P_1^{LP}$ where $P_1^{LP}$ depends on $z_R$, $h_e$ and $v_g$ as well as on $A$ and $\epsilon$. Similarly, for our free-boundary model $P_1 = P_1^{SA}$ where $P_1^{SA}$ depends on $z_R$, $h_e$ and $v_g$ as well as on $S_1$ and $\lambda _{BC}$. With these relationships, we choose to replace our search variable $v_g$ with the search variable $P_1$. The reported mean plus-or-minus the standard deviation for each data set in table 1 then provides upper and lower bounds for the $P_1$ search space.

Before reporting further details on the parameter estimation for the Washburn, Lockington–Parlange and free-boundary models, we give an example solution from our free-boundary model in figure 3 with dimensional heights of the wetting front (solid curve), the saturated/unsaturated front (dashed curve – , i.e. Washburn solution) and the early time linear relationship (dashed–dotted curve), along with the Lago & Araujo (2001) $150\ \mathrm {\mu }{\rm m}$ data (open circles). This plot shows the relationship of the two interface positions, $h_{\ell }^S$ and $h_{\ell }$, along with the early-time linearization. The boundary $h_{\ell }^S$ (which is equivalent to the Washburn model) reaches an equilibrium height while $h_{\ell }$ continues to increase in good approximation to the experimental observations.

Figure 3. An example of our free-boundary model predictions: Lago & Araujo (2001) data for $150\ \mathrm {\mu }{\rm m}$ spheres (open circles); the linear approximation $h_{\ell } - z_R \sim (1/P_1^{SA}) (t - t_R)$ (dash–dotted line); the prediction $h_{\ell } - z_R$ versus $t- t_R$ (solid curve); $h_{\ell }^S - h_{\ell }^S(t_R)$ versus $t-t_R$ (dotted line). Note that $h_{\ell }^S$ is equivalent to the Washburn solution.

4.1. Washburn model predictions

The dimensional form of the classical Washburn model (e.g. using the Lockington–Parlange equations (3.1) and (3.2) with $A \rightarrow \infty$), is

(4.1)\begin{equation} t_{Wash} = \frac{h_e + z_R}{v_g} \ln \left( \frac{z_R + h_e}{z_R + h_e - z_f} \right) - \frac{z_f}{v_g}, \end{equation}

where $t_{Wash}$ is time and $z_f$ is the capillary-rise height relative to fixed position $z=0$. This form has $z_f=0$ when $t_{Wash}=0$ and $z_f \rightarrow z_R + h_e$ as $t_{Wash} \rightarrow \infty$. It is worth noting that for early time $z_f$ and $t_{Wash}$ have the classical relationship

(4.2)\begin{equation} z_f \sim \sqrt{2 v_g (z_R + h_e) t_{Wash}} \sim t_{Wash}^{1/2}. \end{equation}

As previously noted, in order to compare equation (4.1) with the experiments of Lago & Araujo (2001) we need a time- and space-shifted version. The capillary-rise height and time variables used in Lago & Araujo (2001) (see their figure 10a) are the shifted versions $z_f - z_R$ and $t_{Wash} - t_R$, where $t_R$ is defined as the time for which $z_f = z_R$. That is,

(4.3)\begin{equation} t_R= \frac{h_e + z_R}{v_g} \ln \left( \frac{z_R + h_e}{h_e} \right) - \frac{z_R}{v_g}, \end{equation}

and

(4.4)\begin{equation} t_{Wash} - t_R = \frac{h_e + z_R}{v_g} \ln \left( \frac{h_e}{h_e - (z_f - z_R)} \right) - \frac{z_f - z_R}{v_g}. \end{equation}

In this view the early-time dynamics are characterized by the relationship

(4.5)\begin{equation} t_{Wash} - t_R = \frac{z_R}{v_g h_e} (z_f - z_R) + \frac{z_R + h_e}{2 v_g h_e^2} (z_f - z_R)^2 + \cdots. \end{equation}

That is, unless the reservoir height $z_R=0$, the capillary-rise height in the Lago & Araujo (2001) interpretation scales linearly with early times. The dynamics are still classical in the sense that they are characterized by the Washburn equation but the capillary-rise height does not display the commonly referred to $t^{1/2}$ power law (e.g. see their figure 10a and/or our figure 4). Here we define $P_1^{Wash} = z_R/(v_g h_e)$.

Figure 4. Panel (a) shows capillary-rise height versus time predictions with the Lockington–Parlange model using the parameter values listed in table 3 and the Lago & Araujo (2001) data. Panel (b) shows capillary-rise height versus time predictions with the Lockington–Parlange model using the parameter values listed in table 3 and the Delker et al. (1996) data. The triangles indicate power-law slopes of $1$ (left-hand triangle), $1/4$ (upper right-hand triangle) and $1/20$ (lower right-hand triangle).

We note that if given capillary-rise data in the form $(t_{Wash}, z_f)$ the model in (4.1) would allow at most the determination of two independent parameters, $v_g$ and the combination $h_e + z_R$ (that is, one cannot get $z_R$ and $h_e$ individually without additional information). On the other hand, given sufficient data in the form $(t_{Wash}-t_R, z_f- z_R)$ the model in (4.4) does appear to allow the determination of three independent parameters, $v_g$, $h_e$ and $z_R$.

We comment here on upper and lower bounds for the parameters $P_1$, $z_R$ and $h_e$. First, as already mentioned, the search space for the parameter $P_1$ was limited for each data set by the range listed in table 1. Second, bounds for the parameter $z_R$ were identified in two different ways for the Lago & Araujo (2001) and Delker et al. (1996) data sets. Delker et al. (1996) reported in their experiments that $z_R$ was approximately $4$ cm and so for these cases we have simply allowed a small departure ($\pm 0.1$ cm) from this value. In most cases (with the Washburn model, as well as the Lockington–Parlange model, and our free-boundary model) we found that the upper bound $z_R = 4.1$ cm was an active constraint. In general, while extending the range of $z_R$ could allow small numerical reductions in our objective function there were not significant graphical improvements that would seem to justify allowing further departures from $4$ cm estimate reported by Delker et al. (1996). Lago & Araujo (2001) did not report values for their $z_R$. However, they did report ranges for $h_e$ and for $v_g$ (see their figure 11a,b) and since they also report fits based on the Washburn model we identified bounds for their $z_R$ using the Washburn relation $P_1 = z_R/(h_e v_g)$ and the $P_1$ values reported in table 1. The $z_R$ values we found for the Lago & Araujo (2001) data were in all cases within these bounds with no active constraints. Third, for the parameter $h_e$ we chose generous bounds that included the full range of $h_e$ reported for the Lago & Araujo (2001) data sets (see their figure 11a). We expect $h_e \sim \sigma /(\rho g D_b)$ (i.e. Jurin height) where $\sigma$ is surface tension, $\rho$ is fluid density, $g$ is gravity and $D_b$ is the bead diameter. Using $\sigma = 72\,{\rm mN}\,{\rm m}^{-1}$, $\rho = 10^3\,{\rm kg}\,{\rm m}^{-3}$ and $g=9.8\,{\rm m}\,{\rm s}^{-2}$ gives $h_e \sim [4.9, 4.1, 3.5, 2.9]\,{\rm cm}$ for the Lago & Araujo (2001) data and $h_e \sim [4.1, 2.9, 2.0, 1.4]\,{\rm cm}$ for the Delker et al. (1996) data.

The bounds for $P_1$, $z_R$ and $h_e$ are shown in table 2; these are also later used for the Lockington–Parlange model and our free-boundary model. Our Washburn model estimates for $h_e$, $z_R$, $P_1$ and $v_g$ are listed in table 2. It is important to note that for the Washburn model one must also decide how much of the capillary rise data to fit (i.e. make a choice for $N$) as the long-time behaviour of the data is not described by the Washburn model. Recall that the Washburn model predictions have the form of the dotted curve shown in figure 3. Making $N$ too large and trying to fit the capillary-rise data in the anomalous regime will artificially drive up the value of $h_e$. Making $N$ too small and not including enough information about the transition away from the early-time dynamics will not allow identification of $h_e$ and $z_R$. That said, the parameter values in table 2 are generally consistent with various earlier estimates reported by Lago & Araujo (2001) and Delker et al. (1996). Our primary purpose in generating estimates for these parameters from the Washburn model is to have a point of reference when we obtain estimates for these same parameter values for the Lockington–Parlange model and our free-boundary model, both of which are better suited to characterize the long-time dynamics.

Table 2. Fitted Washburn model parameters $h_e$ and $z_R$ for capillary-rise data of Lago & Araujo (2001) and Delker et al. (1996) using values of $P_1$ in the range specified in table 1, and the resultant value $v_g = z_R/(P_1 h_e)$. The range indicated by the brackets are those imposed as constraints during the optimization procedure. The number of points, $N$, used for the fits is listed in parentheses and has been chosen in a range where the estimates for $h_e$ and $z_R$ are relatively insensitive to $N$.

4.2. Lockington–Parlange model predictions

The Lockington–Parlange model with dimensionless capillary-rise height $Z_f$ and time $T$ is given by (3.1) and (3.2). We note that $Z_f$ reaches the dimensionless reservoir height, $Z_R = z_R/(z_R + h_e)$, at a time $T_R$ for which $Q_0=Q_R$. That is,

(4.6)$$\begin{gather} Z_R = \frac{1}{A} \ln \left( \frac{Q_R + 1}{Q_R} \right) + \frac{1}{Q_R + 1}, \end{gather}$$

(4.7)$$\begin{gather}T_R = \ln \left( \frac{Q_R + 1}{Q_R} \right) - \frac{1}{Q_R + 1} + \frac{1}{A \sqrt{1 + 3 \epsilon} } \left[ \frac{1}{Q_R} - \ln \left( \frac{Q_R + 1}{Q_R} \right) \right], \end{gather}$$

where (4.6) determines the value of $Q_R$ for a given reservoir height $Z_R$ and the corresponding time is given by (4.7).

The early time, $0 \le T-T_R \ll 1$, form of $Z_f - Z_R$ is

(4.8)\begin{equation} Z_f - Z_R = Q_R \frac{1 + \dfrac{1}{A} \left( 1 + \dfrac{1}{Q_R} \right)}{1 + \dfrac{1}{A \sqrt{1 + 3 \epsilon} } \left( 1 + \dfrac{1}{Q_R} \right)} \left( T-T_R \right) + O\left( T-T_R \right)^2. \end{equation}

The corresponding dimensional quantities are $z_f - z_R = (Z_f - Z_R) (z_R + h_e)$ and $t-t_R = (T - T_R) (z_R + h_e)/v_g$. Then, the dimensional leading-term in (4.8) is

(4.9)\begin{equation} z_f - z_R \sim v_g \left[ Q_R \frac{1 + \dfrac{1}{A} \left( 1 + \dfrac{1}{Q_R} \right)}{1 + \dfrac{1}{A \sqrt{1 + 3 \epsilon} } \left( 1 + \dfrac{1}{Q_R} \right)}\right] \left( t-t_R \right). \end{equation}

That is, the early-time data is fit by taking $P_1^{LP} = P_1$ where

(4.10)\begin{equation} P_1^{LP} = \frac{1}{v_g Q_R} \frac{1 + \dfrac{1}{A \sqrt{1 + 3 \epsilon} } \left( 1 + \dfrac{1}{Q_R} \right)}{1 + \dfrac{1}{A} \left( 1 + \dfrac{1}{Q_R} \right)}, \end{equation}

which is a non-trivial function of the ratio $h_e/z_R$ (through $Z_R$ and $Q_R$), $A$ and $\epsilon$. When $A \rightarrow \infty$, we find that $Q_R \rightarrow h_e/z_R$ and $P_1^{LP} \rightarrow z_R/(v_g h_e)$ as in the leading behaviour of the Washburn model analogue (4.5).

We fit the full Lockington–Parlange model $z_f - z_R$ versus $t - t_R$ for all times available in the Lago & Araujo (2001) and Delker et al. (1996) data sets. Our search space is over the parameters $h_e$, $z_R$, $A$, $\epsilon$ and $P_1$ with ranges as listed in brackets in table 3. The bounds on $P_1$, $z_R$, and $h_e$ are the same as those already identified in the discussion of the Washburn model. The search space for the parameter $A$ was effectively left unconstrained although for technical reasons we typically set bounds $[0 \ldots 100]$. We used $[0 \ldots 1]$ as the search space for the parameter $\epsilon$ with it in mind that the Lockington–Parlange model assumed $\epsilon \ll 1$. The value of $v_g$ was obtained through the condition $P_1^{LP} = P_1$.

Table 3. Lockington–Parlange model parameters $A$, $\epsilon$, $h_e$ and $z_R$ for capillary-rise data of Lago & Araujo (2001) and Delker et al. (1996) using $P_1$ in the range indicated in table 1 with the resultant value $v_g$ obtained from the relationship (4.10). The range indicated in the $A$, $\epsilon$, $h_e$, $z_R$ and $P_1$ columns correspond to the allowed search space.

The parameter value estimates for the Lockington–Parlange model are listed in table 3. In some of these cases the constraints on the parameter are active, which generally indicates that one can find solutions with numerically smaller objective function values if these bounds were extended. Two examples of note are the solutions for the Delker et al. (1996) 359 $\mathrm {\mu }$m and 510 $\mathrm {\mu }$m cases which have active lower bounds on $h_e$. In these two cases the Lockington–Parlange model predictions appear to fit the data only moderately well (see figure 4). These graphical comparisons could be improved slightly by allowing significantly smaller values of $h_e$. For example, for the 359 $\mathrm {\mu }$m Delker et al. (1996) data the parameter set [$A=9.32$, $\epsilon = 0.0036$, $h_e = 0.39$ cm, $z_R = 4.1$ cm and $P_1 = 11.9$ s cm$^{-1}$] has $F=3.6056$ which is a graphical and numerical improvement over the case listed in table 3 and shown in figure 4. However, the value $h_e = 0.39$ cm appears to be inconsistent (approximately a factor of 10 too small) with what would be expected using the Jurin scaling and the estimate for $h_e$ for the 180 $\mathrm {\mu }$m bead diameter; that is, $\approx 6.7/2 = 3.35$ cm.

It appears from figure 4 that for the majority of these cases, except for largest two bead diameter results of Delker et al. (1996), the Lockington–Parlange model is able to capture essential features of the dynamics from the early-time, classical dynamics all the way through the long-time, anomalous dynamics, evolution. Lago & Araujo (2001) have already reported an estimated range of scaling exponents for their long-time power laws between ${1/20}$ and ${1/4}$. For comparison purposes we have indicated these particular scaling estimates on the graph. For both the Lago & Araujo (2001) data and the Delker et al. (1996) data the anomalous scaling in the long-time dynamics seems to lie somewhere in between these values. The speed and ease of use of the Lockington–Parlange model are advantages over the free-boundary model discussed next, although we shall see that the ability to characterize the Delker et al. (1996) large bead data, appears to be an advantage of the free-boundary model.

4.3. Free-boundary model with Brooks–Corey function predictions

Here we compare our model with Brooks–Corey functions for capillary pressure and permeability to the capillary-rise experiments of Lago & Araujo (2001) and Delker et al. (1996) These solutions were obtained computationally as follows. An initial condition for the partial differential equation for saturation on $\bar {h}_{\ell }^S < \bar {z} < \bar {h}_{\ell }$ was obtained using the similarity solution. In particular, the similarity solution starts with initial conditions $\bar {h}_{\ell }^S(\bar {t}=0)=\bar {h}_{\ell }(\bar {t}=0)=0$ and is used to advance forward in time to a small but non-zero value $\bar {t}^*$ where $\bar {h}_{\ell }^S(\bar {t}^*)$ and $\bar {h}_{\ell }(\bar {t}^*)$ are non-zero and a saturation profile $S(\bar {z},\bar {t}^*)$ from the similarity solution is also obtained. An example saturation curve from the similarity solution used to start the numerical simulations is shown in figure 5(a). Note that this example has dimensional time $t^* = 0.0081$ s and has $h_{\ell }(t^*) > h_{\ell }^S(t^*) >0$. The partial differential equation was discretized in space over the partially saturated region with $N_{pde}+1$ evenly spaced points. We used a range of values of $N_{pde}$ up to $N_{pde}=2048$. Several example solutions for $S(z,t)$ are shown in figure 5(b). We observe that the saturation decreases monotonically at higher positions in the partially saturated media with the variation more pronounced at later times. This example corresponds to a base-case set of parameters from which we explore the solution dependence on the parameters $S_1$ and $\lambda _{BC}$ further below. Parameter estimation for this model was performed in a similar manner to that done for the Washburn model and the Lockington–Parlange model. Details specific to this model involving the interpretation of parameter $P_1$ and the search space used for parameters $S_1$ and $\lambda _{BC}$ are outlined below.

Figure 5. Panel (a) shows the similarity solution $S(z,t)$ plotted against $z- z_R$ at the time value $t= t^* = 0.0081$ s. This function $S(z,t^*)$, along with the values $h_{\ell }(t^*)$ and $h_{\ell }^S(t^*)$, from the similarity solution are used as the initial conditions (at $t=t^*$) for the numerical integration of the free-boundary problem defined in (2.25), (2.26) and (2.27). Note that the red ‘o’ in (a) at $S=1$ corresponds to $h_{\ell }^S(t^*) - z_R$ (i.e. $z=h_{\ell }^S(t^*)$) while the red ‘x’ in this same plot near $S=0$ corresponds to $h_{\ell }(t^*) - z_R$ (i.e. $z=h_{\ell }(t^*)$). Panel (b) shows the full numerical solution $S(z,t)$ – with initial condition indicated in panel (a) – at various times plotted against the shifted coordinate $z - z_R$. In this plot, the ‘o’ marks at $S=1$ correspond to $h_{\ell }^S(t) - z_R$ at the indicated times while the ‘x’ marks near $S=0$ correspond to $h_{\ell }(t) - z_R$ at those times. The parameter values used for this example are $\lambda _{BC}=3.32$, $S_1=0.01$, $h_e = 9.8$ cm, $P_1 = 18.6$ s cm$^{-1}$ and $z_R = 7.5$ cm, which correspond to the values used for the solid curves in figure 6(a,b).

As a first step in the parameter identification for this free-boundary model we identify a parameter that we can associate with the experimentally identified values of $P_1$ (see table 1) which characterize capillary rise dynamics near the time $t = t_R$ for which $h_{\ell } = z_R$. The local expansion of our dimensionless solution is

(4.11)\begin{equation} \bar{h}_{\ell} - \bar{z}_R \sim \left. \frac{{\rm d} \bar{h}_{\ell}}{{\rm d}\bar{t}} \right|_{\bar{t} = \bar{t}_R} ( \bar{t} - \bar{t}_R), \end{equation}

where from (2.26) we have that

(4.12)\begin{equation} {\mathcal{D}} \equiv \left. \frac{{\rm d} \bar{h}_{\ell}}{{\rm d}\bar{t}} \right|_{\bar{t} = \bar{t}_R} ={-} \frac{1}{\lambda_{BC}} S_1^{1 + 1/\lambda_{BC}} \left. \frac{\partial S}{\partial \bar{z}} \right|_{\bar{z} = \bar{h}_{\ell}(\bar{t}_R)} - \frac{1}{1-\bar{z}_R} S_1^{2 + 2/\lambda_{BC}}. \end{equation}

We do not have an analytical expression for the saturation derivative appearing here but it is computed as part of our solution. In dimensional form we have $t - t_R \sim P_1^{SA} (h_{\ell } - z_R)$ where

(4.13)\begin{equation} P_1^{SA} = \frac{h_e + z_R}{h_e v_g} \frac{1}{ {\mathcal{D}} }. \end{equation}

Before discussing additional details of our model comparison with the Lago & Araujo (2001) and Delker et al. (1996) data we explore how the two model parameters $S_1$ and $\lambda _{BC}$ influence the dynamics. Figure 6 shows capillary-rise height $h_{\ell }(t) - z_R$ versus $t-t_R$ as the two parameters $S_1$ and $\lambda _{BC}$ are individually varied. The long-time dynamics are sensitive to both of these parameters. Note that using the equilibrium expressions in (2.17) and (2.18) with the Brooks–Corey expression for capillary pressure in (2.23a,b) gives

(4.14)\begin{equation} h_{\ell}(t \rightarrow \infty) = z_R + h_e \left( \frac{1}{S_1} \right)^{{1}/{\lambda_{BC}}} . \end{equation}

This shows that $h_{\ell }(t \rightarrow \infty ) \rightarrow \infty$ as $S_1 \rightarrow 0$ while $S_1 = 1$ recovers the Washburn equilibrium. These trends can be observed in figure 6. The smaller $S_1$ leads to a higher capillary-rise height. All of our numerical simulations have used non-zero $S_1$ and so all of these predictions eventually reach an equilibrium height. In this particular example in figure 6 the data of Lago & Araujo (2001) for the 150 $\mathrm {\mu }$m bead diameter are shown and their final experimental data point has a capillary-rise height of around $26$ cm. For sufficiently small $S_1$ the theoretical equilibrium height $h_{\ell }(t \rightarrow \infty )$ is well above this value. In general, the Lago & Araujo (2001) data and the Delker et al. data provide capillary-rise data out to times around 24 hours. We are not aware of experimental data for these systems that could test the model predictions, such as the possibility of an eventual plateau in capillary-rise height, beyond these times.

Figure 6. Panel (a) shows four capillary-rise solutions as $S_1$ varies with $\lambda _{BC}=3.32$, $h_e = 9.8$ cm, $P_1 = 18.6$ s cm$^{-1}$ and $z_R = 7.5$ cm. Panel (b) shows four capillary-rise solutions as $\lambda _{BC}$ varies with $S_1 = 0.01$, $h_e = 9.8$ cm, $P_1 = 18.6$ s cm$^{-1}$ and $z_R = 7.5$ cm. For reference we have also included the Lago & Araujo (2001) $150\ \mathrm {\mu }$m data.

Our best parameter estimates in comparison with the Lago & Araujo (2001) and Delker et al. (1996) data are listed in table 4 with the corresponding graphical results shown in figure 7. The search space in each case is indicated and for parameters $h_e$, $z_R$ and $P_1$ these bounds are the same as used previously for the Washburn and Lockington–Parlange models. The upper and lower bounds for the parameters $S_1$ and $\lambda _{BC}$ are also listed and were chosen to effectively keep these parameters unconstrained; none of these constraints are active. Again we do not guarantee that the reported solutions are globally optimal in the prescribed feasible set but in all of these cases the graphical comparison with the data is excellent, especially in the long-time, anomalous, regime.

Table 4. Parameter values for $S_1$, $\lambda _{BC}$, $h_e$ and $z_R$ for our free-boundary model with Brooks–Corey functions for capillary rise data of Lago & Araujo (2001) and Delker et al. (1996) using the values of $P_1$ indicated by the range given in table 1 with the resultant value $v_g$ obtained from the expression (4.13). The range indicated in the $S_1$, $\lambda _{BC}$, $h_e$, $z_R$ and $P_1$ columns correspond to ranges allowed for these parameters in our parameter search.

Figure 7. Panel (a) shows capillary-rise height versus time predictions with our model using the parameter values listed in table 4 and the Lago & Araujo (2001) data. Panel (b) shows capillary-rise height versus time predictions with our model using the parameter values listed in table 4 and the Delker et al. (1996) data.

The solutions reported in tables 3 and 4 along with the graphical counterparts in figures 4 and 7 show that both the Lockington–Parlange and free-boundary models can capture many features of the capillary-rise dynamics for both early- and long-time regimes. We note that there is some variation across models in terms of the predicted values of $h_e$, $v_g$ and $z_R$. As we have shown, the early-time fit establishes a value of the parameter $P_1$ for a given data set, but depending on the model (Washburn, Lockington–Parlange or the free-boundary formulation) the relationships that $h_e$, $v_g$ and $z_R$ have with the parameter $P_1$ differ. That is, other parameters ($A$ and $\epsilon$ for Lockington–Parlange and $\lambda _{BC}$ and $S_1$ for the free-boundary model) also are involved in the corresponding expressions for $P_1$ and, consequently, the fit values of $h_e$, $v_g$ and $z_R$ (if not measured independently in an experiment) are at least slightly model dependent. We have chosen to allow a range for each of these parameters, as noted in the tables. The Lockington–Parlange model perhaps does slightly better numerically for the Lago & Araujo (2001) data while the free-boundary model appears to do better for the Delker et al. (1996) data. We note that the two cases corresponding to the largest bead diameter for the Delker et al. (1996) data show that the capillary-rise dynamics has positive concavity. The prediction of the free-boundary model also has this characteristic for these two cases. The 359 $\mathrm {\mu }$m and $510\ \mathrm {\mu }{\rm m}$ predictions suggest that the larger bead diameter prediction reaches a larger capillary-rise height (compare $S_1=0.0007$ for the 359 $\mathrm {\mu }$m case with $S_1 = 0.001$ for the 510 $\mathrm {\mu }$m case) for times beyond where the experimental data is available. We presume this particular observation may change were further experimental data available for these cases beyond the reported final time.

5. Permeability and capillary pressure measures and their dynamic scaling

In the introduction we reviewed scaling laws initially developed for capillary-rise phenomena in partially saturated porous media by Kim et al. (2017) and Ha & Kim (2020). Their arguments suggested that at early times both the capillary pressure and the permeability scale independently of capillary-rise height. At later times there is a capillarity–gravity balance in which capillary pressure scales linearly with height and the permeability scales inversely proportional to capillary-rise height squared. We explore these ideas further here.

In both of the capillary-rise models presented in this work – the Lockington–Parlange model and our free-boundary model with Brooks–Corey functions – the capillary pressure and permeability are defined throughout the partially saturated region (as well as the completely saturated region) and therefore have both spatial and temporal variation. The question then arises that if, for example, the scaling law $k \sim 1/h^2$ holds, where specifically in the partially saturated region does it hold? Throughout the entire partially saturated region, or in some specific region? By examining the computed permeability function through the identification of representative permeability measures over different regions in the partially saturated media we can address this question. For our free-boundary model we discuss analogue capillary pressure measures as well.

We begin with the Lockington–Parlange model where integrals of the permeability can be expressed in simple terms. Since those results were not part of the original work presented in Lockington & Parlange (2004) we provide details in Appendix B. The four permeability measures we use here correspond to mean permeabilities taken from the bottom quarter to the top quarter of the partially saturated region. In particular, we define

(5.1)$$\begin{gather} \bar{K}^{eff}_{I}(t) ={-}\frac{4}{\Delta Z} \int_{0}^{{-}u_1^*} \frac{\exp [A (1+\epsilon) u]}{1+Q_0 \exp [{-}A u] } \,{\rm d}u, \end{gather}$$

(5.2)$$\begin{gather}\bar{K}^{eff}_{II}(t) ={-}\frac{4}{\Delta Z} \int_{{-}u_1^*}^{{-}u_2^*} \frac{\exp [A (1+\epsilon) u]}{1+Q_0 \exp [{-}A u] } \,{\rm d}u, \end{gather}$$

(5.3)$$\begin{gather}\bar{K}^{eff}_{III}(t) ={-}\frac{4}{\Delta Z} \int_{{-}u_2^*}^{{-}u_3^*} \frac{\exp [A (1+\epsilon) u]}{1+Q_0 \exp [{-}A u] } \,{\rm d}u, \end{gather}$$

(5.4)$$\begin{gather}\bar{K}^{eff}_{IV}(t) ={-}\frac{4}{\Delta Z} \int_{{-}u_3^*}^{-\infty} \frac{\exp [A (1+\epsilon) u]}{1+Q_0 \exp [{-}A u] } \,{\rm d}u, \end{gather}$$

where $\Delta Z = Z_f - Z_e$, $Z_e = 1/(Q_0 + 1)$ and $Z_f$ and $Q_0$ are as in (3.1) and (3.2). For $i= 1,2,3$ we have

(5.5)\begin{equation} {-}u_i = \frac{1}{A} \ln \left[{-}Q_0 + (1+Q_0) {\rm e}^{\omega_i} \right],\quad \omega_i ={-} \frac{i}{4} A \Delta Z, \end{equation}

which define the integration limits. We also define an overall effective permeability for the whole partially saturated region which has $\bar {K}^{eff}(t) = ( \bar {K}^{eff}_{I}(t) + \bar {K}^{eff}_{II}(t) + \bar {K}^{eff}_{III}(t) + \bar {K}^{eff}_{IV}(t))/4$. We discuss these results below in the context of similar permeability measures from the free-boundary model.

For the free-boundary model we use the same choices for effective permeability. In this approach, the saturation variable is computed numerically and the corresponding permeability is defined by the Brooks–Corey model in (2.23a,b). In particular, dividing the partially saturated porous material into four quarters leads to

(5.6)$$\begin{gather} \bar{K}^{eff}_{I}(t) = \frac{4}{\bar{h}_{\ell}- \bar{h}_{\ell}^S} \int_{\bar{h}_{\ell}^S}^{\bar{h}_{\ell}^S + (\bar{h}_{\ell}- \bar{h}_{\ell}^S)/4} \bar{k}_{r\ell} (S(z,t)) \,{\rm d}z, \end{gather}$$

(5.7)$$\begin{gather}\bar{K}^{eff}_{II}(t) = \frac{4}{\bar{h}_{\ell}- \bar{h}_{\ell}^S} \int_{\bar{h}_{\ell}^S + (\bar{h}_{\ell}- \bar{h}_{\ell}^S)/4}^{\bar{h}_{\ell}^S + (\bar{h}_{\ell}- \bar{h}_{\ell}^S)/2} \bar{k}_{r\ell} (S(z,t)) \,{\rm d}z, \end{gather}$$

(5.8)$$\begin{gather}\bar{K}^{eff}_{III}(t) = \frac{4}{\bar{h}_{\ell}- \bar{h}_{\ell}^S} \int_{\bar{h}_{\ell}^S + (\bar{h}_{\ell}- \bar{h}_{\ell}^S)/2}^{\bar{h}_{\ell}^S + 3(\bar{h}_{\ell}- \bar{h}_{\ell}^S)/4} \bar{k}_{r\ell} (S(z,t)) \,{\rm d}z, \end{gather}$$

(5.9)$$\begin{gather}\bar{K}^{eff}_{IV}(t) = \frac{4}{\bar{h}_{\ell}- \bar{h}_{\ell}^S} \int_{\bar{h}_{\ell}^S + 3(\bar{h}_{\ell}- \bar{h}_{\ell}^S)/4}^{\bar{h}_{\ell}} \bar{k}_{r\ell} (S(z,t)) \,{\rm d}z. \end{gather}$$

We also examine $\bar {K}^{eff}(t) = ( \bar {K}^{eff}_{I}(t) + \bar {K}^{eff}_{II}(t) + \bar {K}^{eff}_{III}(t) + \bar {K}^{eff}_{IV}(t))/4$ which is the mean over the whole partially saturated region. Analogue measures for capillary pressure with $\bar {k}_{r\ell }$ replaced by $\bar {p}_c$ will also later be discussed.

Figure 8 shows two plots of effective permeability versus capillary-rise height. Figure 8(a) corresponds to the Lockington–Parlange model for the example based on the Lago & Araujo (2001) 150 $\mathrm {\mu }$m data. Figure 8(b) corresponds to the free-boundary model for the example based on the Delker et al. (1996) 359 $\mathrm {\mu }$m data. In both examples from top to bottom are the permeability estimates $\bar {K}^{eff}_{I}$ (cyan), $\bar {K}^{eff}_{II}$ (blue), $\bar {K}^{eff}_{III}$ (green) and $\bar {K}^{eff}_{IV}$ (red). The black dashed line shows the mean permeability, $\bar {K}^{eff}$, over the entire partially saturated region. The left-hand vertical dashed line marks approximately the boundary between the early-time classical dynamics and the transition region while the right-hand vertical dashed line marks approximately the boundary between the transition region and the long-time anomalous dynamics. Compare the Lago & Araujo (2001) 150 $\mathrm {\mu }$m capillary-rise data and theoretical solution in figure 4 for the Lockington–Parlange model and the Delker et al. (1996) 359 $\mathrm {\mu }$m data and theoretical solution in figure 7 for the free-boundary model. A first observation in figure 8 is that for early times (capillary-rise height to the left of the left-hand dashed line) each of the permeability measures is approximately independent of capillary-rise height. This is consistent with the early-time scaling predictions of Kim, Ha and coworkers. A second observation from figure 8 is that for long times (capillary-rise height to the right of the right-hand vertical dashed, black, line) there is a significant variation of each of the permeability measures with respect to height. The overall mean permeability, $\bar {K}^{eff}$ (black dashed curve), clearly varies strongly with $h$. For the Lockington–Parlange model this overall mean appears to scale less strongly with $h$ than the Kim et al. (2017) scaling $k \sim 1/h^2$. For free-boundary model example in figure 8(b) the $k \sim 1/h^2$ appears close to that observed in the lower quarter of the partially saturated regime and also for overall mean permeability in the high range of capillary-rise height. We see that this the overall mean is dominated by the behaviour of the permeability in the lowest region of the partially saturated porous media (cyan curve, $\bar {K}^{eff}_{I}$). Higher up in the partially saturated region, as measured by $\bar {K}^{eff}_{II}$, $\bar {K}^{eff}_{III}$ and $\bar {K}^{eff}_{IV}$, the permeability varies more strongly with capillary-rise height. That is, here $k \sim 1/h^\nu$ with $\nu > 2$. Qualitatively similar results are found for the other Lago & Araujo (2001) data as well as for the Delker et al. (1996) data. Additionally, we note that the choice of sampling by quarter-regions is representative of other choices for partitioning the partially saturated region, including pointwise sampling in space. These analyses show that the permeability, as suggested in the scaling arguments of Kim et al. (2017) and Ha & Kim (2020) does indeed change considerably from early time to late time and seems to be an essential difference distinguishing the early-time classical capillary rise dynamics from the long-time anomalous dynamics.

Figure 8. Panel (a) shows different permeability measures for the Lockington–Parlange model versus capillary-rise height for Lago & Araujo (2001) $150\ \mathrm {\mu }$m data. Panel (b) shows different permeability measures for the free-boundary model versus capillary-rise height for Delker et al. (1996) $359\ \mathrm {\mu }$m data. In each plot, the thin vertical dashed lines mark the approximate transition heights between the early-, intermediate- and long-time regimes.

5.1. Comparison with Kim et al. (2017) data

As a final comparison of the model predictions we examine data from Kim et al. (2017) in their figure 2(c) showing capillary rise dynamics of turpentine in a cellulose sponge. This particular data set shows very clearly both the classical $t^{1/2}$ power law at early times as well as the $t^{1/4}$ power law at later times. As such, it offers a good opportunity to test the model predictions in terms of capillary-rise dynamics as well as scaling arguments related to permeability and capillary pressure posed by Kim et al. (2017) and Ha & Kim (2020).

In figure 9(a) we show this Kim et al. (2017) data along with the predictions of the Lockington–Parlange model (red curve) and our free boundary model (blue curve). For the Lockington–Parlange model we used $A=1.2918$, $\epsilon = 0.001425$, $h_e = 0.9015$ cm, $v_g = 0.4309$ cm s$^{-1}$ and $z_R=0.01031$ cm. We note that the value of $h_e$ falls approximately where Kim et al. (2017) estimate the Jurin height marking the distinction between early stages (solid squares) and late stages (open squares) while $z_R$ is approximately zero. The free-boundary model uses $S_1 = 0.015$, $\lambda _{BC}=1.5$, $h_e=1.6$ cm, $v_g = 0.5625$ cm s$^{-1}$ and $z_R=0.01$ cm. Both models reasonably capture both early- and long-time dynamics.

Figure 9. Panel (a) shows two capillary-rise height predictions (Lockington–Parlange model, red; free-boundary problem, blue) plotted against the data for capillary rise of turpentine in a cellulose sponge from Kim et al. (2017) (open and solid squares; see also their figure 2c). Panels (b) and (c) show, for the free-boundary model, various effective capillary pressure measures and permeability measures, respectively, versus capillary-rise height taken over the same four quarters of the partially saturated media.

We recall here two key components of the scaling arguments of Kim, Ha and coworkers that lead to the long-time $t^{1/4}$ power law prediction for capillary-rise height. First, in the partially saturated media there is a balance between capillarity and gravity in partially filled macropores so that $p \sim \sigma /r_g \sim \rho g h$ where $r_g$ is a length scale for the meniscus at the macropore-scale. In particular the capillary pressure is predicted to scale linearly with capillary-rise height, $h$. Second, they argue that the permeability scales with this macropore-scale which together with the capillary-pressure balance gives the prediction that the permeability scales like $k \sim r_g^2 \sim (\sigma / (\rho g h))^2$. We can test these predictions with our solutions.

For the free-boundary model solution shown in figure 9(a) we compute capillary pressure measures of the same form as those in (5.6)–(5.9) except with $\bar {k}_{rl}(S)$ replaced with $\bar {p}_c(S)$ (both using (2.23a,b)). Figure 9(b) shows mean capillary pressures over the lower quarter (cyan), second quarter (blue), third quarter (green), top quarter (red) and entire partially saturated region (black, dashed) versus capillary-rise height. In the early-time regime all of the capillary pressure measures appear to be effectively independent of capillary-rise height. In the long-time regime the overall mean capillary pressure as well as any of the means in the upper three-quarters of the partially saturated region all appear to closely follow a linear scaling with capillary-rise height, confirming a balance between capillarity and gravity in the partially filled macropores. Both early- and long-time observations are consistent with the scaling arguments put forward by Kim et al. (2017) and Ha & Kim (2020). We note that in the lower quarter of the partially saturated media the effective capillary pressure appears to scale more weakly with $h$ (cyan curve); this reflects the fact that at the very bottom of the partially saturated regime the pressure limits to the constant reference pressure $\rho g h_e$.

Also for the free-boundary model solution shown in figure 9(a) we compute mean permeabilities as defined in (5.6)–(5.9): $\bar {K}^{eff}_{I}$ (cyan), $\bar {K}^{eff}_{II}$ (blue), $\bar {K}^{eff}_{III}$ (green), $\bar {K}^{eff}_{IV}$ (red) and $\bar {K}^{eff}$ (black, dashed). These are shown in figure 9(c). We observe that for early times (on the left-hand side of the vertical dashed blue line) all permeability measures are roughly independent of capillary-rise height, consistent with the scaling arguments of Kim et al. (2017) and Ha & Kim (2020). We observe that for the long-time regime (on the right-hand side of the vertical dashed black line) the permeability throughout the partially saturated region decreases strongly with capillary-rise height. For reference, scaling trends $k \sim 1/h^2$ and $k \sim 1/h$ are indicated by the upper and lower triangles. The overall mean permeability, $\bar {K}^{eff}$, (black dashed curve) appears to be dominated by the trends in the lower portion of the partially saturated regime (cyan curve), and both appear to scale with $1/h^\nu$ where $1 < \nu < 2$. The mean permeabilities in the upper portion of the partially saturated regime scale more strongly with height suggesting that the Kim et al. (2017) and Ha & Kim (2020) permeability scaling $k \sim 1/h^2$ is achieved over some intermediate region in the lower portion of the partially saturated media.

Taken together, these observed relationships for the mean capillary pressure and mean permeability with capillary-rise height in the partially saturated porous media, which we compute directly from solutions to our model, appear to support the scaling arguments of Kim et al. (2017) and Ha & Kim (2020).

6. Conclusions

We have explored predictions of two models of capillary-rise of a fluid in partially saturated porous media. One model is the Lockington & Parlange (2004) model that is based on a travelling-wave approximation in Richards’ equation. The second model, which we refer to as the free-boundary model, also arrives at Richards’ equation with two moving boundaries corresponding to the interface between fully saturated and partially saturated regions in the porous material and the wet/dry interface at the top of the partially saturated region. Both models use saturation-dependent functions for capillary pressure and permeability. We compare the model predictions with three different sets of experimental data from Delker et al. (1996), Lago & Araujo (2001) and Kim et al. (2017). These experiments have documented both the classical, early-time, dynamics of capillary rise as well as the long-time, anomalous, dynamics. We give parameter estimates for each model in comparison with these experiments and find that these generally give good graphical comparison with the data over the full range of dynamics.

We build on these results to give predictions for the relationship of effective permeability and capillary pressure measures to the capillary-rise height. Such relationships have recently been identified in scaling arguments by Kim et al. (2017) and Ha & Kim (2020). In particular, for the anomalous regime, those authors predict a capillarity–gravity balance in which the capillary pressure scales linearly with capillary-rise height and the permeability scales inversely proportional to the square of the capillary-rise height. Together these give anomalous capillary-rise behaviour of the form $h \sim t^{1/4}$. In both of the computational models we explore in this work, we have been able to numerically compute different permeability measures in the form of mean permeabilities taken over different regions of the partially saturated media. Additionally, for our free-boundary model we have computed analogous capillary pressure measures. The capillary–gravity balance is evident in our observed linear relation between capillary-pressure and capillary-rise height throughout the upper portion of the partially saturated media. We observe that the permeability scaling inversely proportional to capillary-rise height squared appears to be achieved in the lower portion of the partially saturated porous media. All of these observations support the scaling arguments of Kim et al. (2017) and Ha & Kim (2020).