2.1. An averaged continuum multi-phase model for isothermal drying

We model the two-phase flow and phase change of the liquid and a mixture of vapour and inert gas (e.g. air) through a uniform-porosity porous medium using an averaged continuum model, as derived through volume averaging arguments by Whitaker (Reference Whitaker1977). Our starting point is the volume-averaged equations summarised in section VI of Whitaker (Reference Whitaker1977). The assumptions made in their derivation are also listed in that section. For our setting, the most important of these assumptions are that the liquid phase density is constant, the porous medium is homogeneous, the gas flow does not affect the liquid momentum balance except via the capillary pressure relation and that this capillary pressure is well defined. In particular, we emphasise that both the liquid and gas phases are assumed to be continuous. We further simplify this model of Whitaker (Reference Whitaker1977) by assuming that the porosity is constant, the drying process is isothermal (so that there is no temperature variation), the vapour–gas mixture has a constant intrinsic average density, and that the jump between the intrinsic averaged gas mixture and liquid pressures can be described by a capillary pressure. Under these assumptions, the dependent variables in our model are the liquid saturation, $S$ (the volume fraction of liquid within the pore space, so that $1-S$ is the volume fraction of the vapour–gas mixture within the pore space), the liquid and gas-mixture pressures, $p^L$ and $p^G$ , respectively, the liquid and gas-mixture Darcy fluxes, $\boldsymbol{U}^L$ and $\boldsymbol{U}^G$ , respectively, and the vapour mass density, $\rho ^v$ , within the gas phase. Specifically, $\rho _v$ , $p^G$ and $p^L$ are the intrinsic average vapour density and fluid pressures (i.e. averaged only across the relevant phase). The model consists of Darcy’s law and a conservation-of-mass constraint for each phase, incorporating the phase change with rate $\mathcal{M}$ , an advection–diffusion equation for the transport of vapour and the assumed capillary pressure relation. The model is

(2.1a) \begin{align}\boldsymbol{U}^L&=-\frac {k_0}{\mu ^L}k^L(S)\big({\boldsymbol\nabla} p^L+\rho ^Lg\boldsymbol{e}_z\big), \end{align}

(2.1b) \begin{align}\boldsymbol{U}^G&=-\frac {k_0}{\mu ^G}k^G(S)\big({\boldsymbol\nabla} p^G +\rho ^Gg\boldsymbol{e}_z\big), \end{align}

(2.1c) \begin{align}\phi \rho ^L S_t+\rho ^L{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \boldsymbol{U}^L&=-\mathcal{M}(\rho ^v,S), \end{align}

(2.1d) \begin{align}-\phi \rho ^G S_t+\rho ^G{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \boldsymbol{U}^G&=\mathcal{M}(\rho ^v,S), \end{align}

(2.1e) \begin{align}\phi (1-S)\rho ^v_t+\boldsymbol{U}^G\,{\boldsymbol\cdot}\, {\boldsymbol\nabla} \rho ^v&=D{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, (D^{\textit{eff}}(S){\boldsymbol\nabla} \rho ^v)+\left (1-\frac {\rho ^v}{\rho _G}\right )\mathcal{M}(\rho ^v,S), \end{align}

(2.1f) \begin{align}p^L&=p^G-p_c(S), \end{align}

with subscript $t$ denoting the partial derivative with respect to time $t$ . Here, $\phi$ is the porosity, $k_0$ is the absolute permeability of the medium, $\rho ^L$ and $\rho ^G$ are the liquid and gas-mixture densities, $D$ is the diffusivity of vapour in air through the dry porous medium, and $\mu ^L$ and $\mu ^G$ are the liquid and gas-mixture dynamic viscosities, all of which we assume to be constant. The liquid and gas-mixture phase pressures differ by the capillary pressure, $p_c(S)$ , which we assume to be a known function of $S$ . The relative permeabilities, $k^L(S)$ and $k^G(S)$ , and the effective diffusivity through the gas-occupied pore space, $D^{\textit{eff}}(S)$ , are also assumed to be known functions of $S$ only, while the mass-transfer rate between the phases, $\mathcal{M}(\rho ^v,S)$ , is expected to depend on both the local vapour density and liquid saturation. Of course, all of these constitutive laws, $p_c$ , $k^L$ , $k^G$ , $D^{\textit{eff}}$ and $\mathcal{M}$ , additionally depend on the pore-scale geometry (sometimes encoded simply by a dependence on the porosity $\phi$ ) and on the properties of the solid–liquid–gas combination, such as the contact angle of the contact lines. We assume the porous material to be homogeneous and the solid structure rigid and unchanging during the drying process, so that the dependence of the constitutive laws on these properties does not change during the course of the drying process. We explicitly include the constitutive laws’ dependence on $S$ and $\rho ^v$ since these vary spatially and temporally during the drying, as discussed in the following section.

2.2. Constitutive relations

As discussed previously, the form of the constitutive laws (or effective parameters) $p_c(S)$ , $k^L(S)$ , $k^G(S)$ , $D^{\textit{eff}}(S)$ and $\mathcal{M}(\rho ^v,S)$ depend on the properties of the solid porous matrix, as well as the liquid and inert-gas pair considered. Expressions for some of these constitutive laws are given by Whitaker (Reference Whitaker1977), although these are difficult to evaluate explicitly as they are given in terms of the (unknown) microscale flow and transport. Others, such as $p_c(S)$ and $D^{\textit{eff}}(S)$ , do not follow from the volume averaging and must be determined empirically. Alternatively, upscaling via homogenisation methods typically results in more tractable methods by which the constitutive relations may be computed for a given pore geometry, although, as discussed in § 1, this is still only likely to be possible under strong restrictions on the pore geometry and phase distribution within the pore space, as in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2022). In practice, constitutive relations are often fit empirically or assumed based on heuristic arguments to solve models such as (2.1) numerically. In our analysis in § 5, we will leave the constitutive relations unfixed as far as is possible to maintain generality and investigate how the qualitative behaviour of the drying system depends on the properties of their functional forms. Indeed, we will later show that the choice of these functional forms can strongly affect the type of drying behaviour observed. In this section, we discuss the typical properties of the constitutive laws and some simple examples.

Generally, we require that $p_c$ is a monotonic decreasing function of $S$ , with $p_c=0$ at $S=1$ and $p_c\rightarrow \infty$ as $S$ becomes small. The relative permeabilities of each phase should be an increasing function of the volume fraction of that phase, so that $k^L(S)$ is increasing and $k^G(S)$ is decreasing, with $k^L(1)=k^G(0)=1$ and $k^L(0)=k^G(1)=0$ . In reality, there may be significant hysteresis between wetting and drying regimes, due to the history dependence of the pore-scale flow (Blunt Reference Blunt2017). Since we are interested only in drying, we will assume a single relation for each of $p_c(S)$ , $k^L(S)$ , $k^G(S)$ . Commonly used forms with these properties include

(2.2) \begin{align} k^L(S)=S^a, \qquad k^G(S)=(1-S)^b, \qquad p_c(S)=S^{-c}-S^{d}, \end{align}

for some powers $a,\,b,\,c,\,d\gt 0$ (Fowler Reference Fowler2011; Blunt Reference Blunt2017). The model of Van Genuchten (Reference Van Genuchten1980) is also popular in simulating drying (Ouedraogo et al. Reference Ouedraogo, Cherblanc, Naon and Bénet2013; Li et al. Reference Li, Vanderborght and Smits2019). We will keep the constitutive laws general as much as possible through our analysis, but for examples and simulations, will use (2.2) with $a=b=3$ and $c=d=1$ (Fowler Reference Fowler2011) unless otherwise stated.

The evaporation rate $\mathcal{M}$ is also assumed to be some known function of $\rho ^v$ and $S$ . From the volume-averaging argument, we might expect to be able to decompose

(2.3) \begin{equation} \mathcal{M}=\mathcal{A}(S)E(\rho ^v,S), \end{equation}

where $\mathcal{A}(S)$ is a measure of the surface area to volume ratio of the liquid–gas interface on the microscale and $E(\rho ^v,S)$ is the local evaporative flux rate on the microscale. This $\mathcal{A}(S)$ is sometimes referred to as the specific interfacial area of the surface in the drying literature (Ahmad et al. Reference Ahmad, Talbi, Prat, Tsotsas and Kharaghani2020; Wu & Prat Reference Wu and Prat2022). Ahmad et al. (Reference Ahmad, Talbi, Prat, Tsotsas and Kharaghani2020) used simulations of a pore-network model to approximate $\mathcal{A}(S)$ , and found for their simulations that $\mathcal{A}\rightarrow 0$ as $S\rightarrow 0$ and $\mathcal{A}$ was non-monotone in $S$ . The local evaporation rate is often taken to be linear in the vapour density $\rho ^v$ (or equivalently, since we assume an isothermal system, in the vapour partial pressure). An example of such a linear model for $E$ is a Hertz–Knudsen law (Hertz Reference Hertz1882)

(2.4) \begin{equation} E(\rho ^v)=C(\rho ^v_*-\rho ^v), \end{equation}

for constant $C$ (usually inversely proportional to the temperature), and where $\rho ^v_*$ is the saturation vapour density. As well as the local vapour density, the local rate $E$ might also depend on the liquid–gas interface curvature and, therefore, the pore-scale geometry of the porous medium. It is for this reason that $E$ may also depend on the saturation $S$ , as a proxy for this pore-scale geometry dependence. In general, $E$ might depend on surface adsorption or dual-porosity effects, which could also be modelled via a dependence on $S$ . We will assume throughout that the porous medium is simple and non-reactive with the liquid – for instance, an array of glass beads – and therefore that $E$ is independent of $S$ .

The effective diffusivity $D^{\textit{eff}}(S)$ of vapour through the porous medium is expected to decrease with $S$ as the pore-volume occupied by the gas mixture is reduced. We require $D^{\textit{eff}}(1)=0$ and $D^{\textit{eff}}(0)=1$ (the effective diffusivity of vapour through the porous medium with no liquid present has been incorporated into the dimensional constant $D$ ). In general, the form of $D^{\textit{eff}}(S)$ will also depend on the pore-scale geometry (Bruna & Chapman Reference Bruna and Chapman2015) but for simplicity, we will take $D^{\textit{eff}}(S)=(1-S)^{1/2}$ in numerical simulations, unless otherwise stated.

2.3. Problem of interest and appropriate boundary conditions

In this paper, we consider drying from an effectively one-dimensional porous medium with an impermeable base at $z=0$ and open to the surrounding air at the top, $z=L$ , so that $L$ is the depth of the porous medium. (If there were no gravitational effect, this geometry also captures the symmetric drying of a porous medium between $z=(-L,L)$ , with a line of symmetry at $z=0$ .) While the porous medium is in reality three-dimensional, by imposing initial and boundary conditions that are independent of $x$ and $y$ , we expect the dependent variables to only depend on $z$ spatially. Furthermore, the drying process, which we will see is limited by diffusion of vapour out of the material, is expected to be stable, and thus we do not anticipate higher dimensional effects becoming important.

By eliminating $U^L$ , $U^G$ (the $z$ -directional scalar Darcy fluxes) and $p^L$ using (2.1a ), (2.1b ) and (2.1f ), our model may be viewed as three, spatially second-order, (diffusive) equations for $S$ , $p^G$ and $\rho ^v$ . We therefore require one boundary condition on each of $S$ , $p^G$ and $\rho ^v$ at each edge of the domain, $z=0$ and $z=L$ , as well as a consistent initial condition for each of these three variables.

At the impermeable base, we impose no flux of liquid, gas or vapour, so that

(2.5) \begin{equation} U^L=U^G=DD^{\textit{eff}}\rho ^v_z=0 \qquad \text{ at }z=0. \end{equation}

Here, we use subscript $z$ to denote derivative with respect to $z$ . These conditions (2.5) may be interpreted as corresponding to one condition on each of $S$ , $p^G$ and $\rho ^v$ , respectively.

At the surface of the material, there can be no flux of liquid out of the material and the gas pressure is atmospheric, so that

(2.6) \begin{equation} U^L=0, \quad p^G=p^G_{\textit{atm}} \qquad \text{ at }z=L. \end{equation}

These may be viewed as one boundary condition on each of $S$ and $p^G$ . The final boundary condition for $\rho _v$ on $z=L$ depends on how vapour is transported through the surrounding air. Most generally, we might couple the drying porous medium with a model for the external vapour transport. For simplicity, we consider two simpler alternatives. First, we might impose that the vapour density is at a constant value, an atmospheric humidity, at the top of the porous medium, so that

(2.7) \begin{equation} \rho ^v=\rho ^v_\infty \qquad \text{ at }z=L, \end{equation}

which is valid if the vapour is very quickly removed from the surface of the porous medium, for instance, if there is a fast flow of dry air across it. Alternatively, we might impose a mass flux or Newton-cooling law, with a mass transfer coefficient $m$ , so that the flux of vapour out of the top of the porous medium is given by

(2.8) \begin{equation} \rho ^vU^G-D D^{\textit{eff}}\rho ^v_z=m(\rho ^v-\rho ^v_\infty ) \qquad \text{ at }z=L. \end{equation}

In the limit of fast mass transfer through the surrounding air, $m\rightarrow \infty$ , this mass-flux boundary condition (2.8) reduces to the Dirichlet condition (2.7). Either of (2.7) or (2.8) completes the set of boundary conditions on the model.

Initially, we will assume known saturation, vapour density and gas pressure profiles. Typically, we might assume these are uniform with $S=S_0\lt 1$ , $p^G=p^G_{\textit{atm}}$ and $\rho ^v=\rho ^v_\infty$ .

The correct boundary conditions to impose at the surface of a drying porous medium are widely debated in the drying literature, especially in relation to simplified continuum models for which the liquid and vapour are assumed to be in local equilibrium (so $\rho ^v=\rho ^v_*$ throughout), or the conservation of mass equations are combined into a single diffusive equation for the saturation (Li et al. Reference Li, Vanderborght and Smits2019; Talbi & Prat Reference Talbi and Prat2021; Wu & Prat Reference Wu and Prat2022). The conditions stated previously, (2.6) and either of (2.7) or (2.8), are physically motivated and by systematically reducing our model in § 5, we derive the corresponding boundary conditions for our simplified models. It will be clear what the correct boundary conditions should be for our reduced models because they are the appropriate reduction of the boundary conditions for the above mentioned full model.

5.1. Stage 1: the ‘constant rate period’

Since $\textit{Pe}\gg 1$ , we see from (3.3e ) that the evaporation rate is zero and the vapour density must take its saturation value $\rho =1$ everywhere except in a thin boundary layer, of width $1/\sqrt {\textit{Pe}}\ll 1$ , at the surface $z=1$ of the porous medium (this boundary-layer width is determined by a balance in (3.3e )). Since evaporation only occurs in this thin boundary-layer region, we see from (3.3c )–(3.3d ) that there can be at most an $O(1/\sqrt {\textit{Pe}})$ flow generated in either fluid. The time scale of the drying is determined by the time scale of changing saturation over the entire domain. Balancing the two terms on the left of (3.3c ) over the $O(1)$ length scale, we see that the time scale is longer, $O(\sqrt {\textit{Pe}})$ , since the flux of liquid is smaller. This is intuitive, since evaporation only takes place in a thin layer, so the time taken to evaporate all of the liquid must be proportionately longer. A schematic diagram illustrating the structure of the solution in stage 1 is shown in figure 3.

Figure 3. Schematic illustrating the solution structure during stage 1.

Motivated by these scaling choices, to study stage 1, we make the change of variables

(5.1) \begin{align} t=\sqrt {\textit{Pe}}\, \tau , \qquad U^G=\frac {1}{\sqrt {\textit{Pe}}}V^G, \qquad U^L=\frac {1}{\sqrt {\textit{Pe}}}V^L. \end{align}

5.1.1. Outer region

On an $O(1)$ length scale and with the rescalings (5.1), the model (3.3) then becomes

(5.2a) \begin{align}\textit{Ca}\textit{Pe}^{-1/2}V^L&=-k^L(S)\big ( p^L_z+ \textit{Bo}\big )\kern-2pt, \end{align}

(5.2b) \begin{align}\xi ^{-1}\textit{Ca}\textit{Pe}^{-1/2}V^G&=-k^G(S)\big ( p^G_z+\delta \textit{Bo}\big )\kern-2pt, \end{align}

(5.2c) \begin{align}\textit{Pe}^{-1/2}\big (S_\tau +V^L_z\big )&=-\mathcal{A}(S)E(\rho ), \end{align}

(5.2d) \begin{align}\textit{Pe}^{-1/2}\big (-\delta S_\tau +V^G_z\big )&=\mathcal{A}(S)E(\rho ), \end{align}

(5.2e) \begin{align}\textit{Pe}^{-1/2}\big (\delta \nu (1-S)\rho _\tau +\nu V^G\rho _z\big )&=\textit{Pe}^{-1}\big (D^{\textit{eff}}(S)\rho _z\big )_z+(1-\nu \rho )\mathcal{A}(S)E(\rho ). \end{align}

We see from (5.2e ) that

(5.3) \begin{equation} O\big(\textit{Pe}^{-1/2}\big)=\left (1-\nu \rho \right )\mathcal{A}(S)E(\rho ). \end{equation}

Thus, since $\rho \leqslant1$ , $\nu \lt 1$ and $S\gt 0$ , we must have that the vapour density is at its saturation value $\rho =1$ , so that $E=0$ to leading order in $\textit{Pe}^{-1/2}$ . Indeed, looking for an $O(\textit{Pe}^{-1/2})$ correction to $\rho$ in (5.2e ) results in that correction being zero (since the leading-order $\rho =1$ has no derivatives) and, arguing inductively in this way, over an order 1 length scale, we see that

(5.4) \begin{equation} \rho =1-\text{ exponentially small correction}. \end{equation}

Since the evaporation rate is therefore exponentially small, from (5.2d ), along with the assumption that $\delta \ll 1$ , we must have a uniform gas flux $V^G_z=0$ and so, using the boundary condition $V^G=0$ at $z=0$ , we see that

(5.5) \begin{equation} V^G=0, \end{equation}

while from (5.2c ), we find that

(5.6) \begin{equation} S_\tau +V^L_z=0. \end{equation}

Since $V^G=0$ , so long as $\delta \textit{Bo}\ll 1$ , then from (5.2b ), we must have a uniform gas pressure $p^G_z=0$ through the outer region.

Furthermore, if $\textit{Ca}\ll\textit{Pe}^{1/2}$ and $\textit{Bo}\ll 1$ , then from the liquid Darcy law (5.2a ), we see that, since the liquid flow is slow, the liquid pressure $p^L=p^G-p_c(S)$ is also uniform to leading order. Thus, we must have the expected uniform (although time-varying) saturation, $S$ , throughout this outer region.

This uniform saturation profile simplifies the problem greatly. By integrating (5.6) over the outer region, we see that the rate of change of saturation,

(5.7) \begin{equation} S_\tau =-\big[V^L\big]^1_0=-V^{L,BL}, \end{equation}

is precisely the flow of liquid $V^{L,BL}$ into the boundary layer at the surface. Furthermore, since $S_\tau$ is independent of $z$ , the liquid flux,

(5.8) \begin{equation} V^L=-S_\tau z=V^{L,BL}z, \end{equation}

is linear in $z$ .

In summary, at leading order through the outer region, there is no evaporation and no gas flow, the saturation, vapour density and gas pressure are all uniform (with the vapour density at its saturation point 1), and the liquid flux is linear in $z$ , determined by the rate at which liquid is drawn into the boundary layer. To understand the flow rate $V^{L,BL}$ into the boundary layer, we must examine the boundary-layer processes in more detail.

5.1.2. Boundary layer at $z=1$

In addition to the scalings (5.1), within the boundary layer at $z=1$ , we set

(5.9) \begin{equation} z=1+\frac {1}{\sqrt {\textit{Pe}}}X, \end{equation}

noting that the top of the porous medium at $z=1$ is at $X=0$ , and as $X\rightarrow -\infty$ , we leave the boundary layer and enter the outer region. In the boundary layer, (3.3) become

(5.10a) \begin{align}\textit{Ca}\textit{Pe}^{-1}V^L&=-k^L(S)\left ( p^L_X+ \textit{Pe}^{-1/2} \textit{Bo} \right ), \end{align}

(5.10b) \begin{align}\xi ^{-1}\textit{Ca}\textit{Pe}^{-1}V^G&=-k^G(S)\left ( p^G_X+\delta\textit{Pe}^{-1/2} \textit{Bo} \right ), \end{align}

(5.10c) \begin{align}\textit{Pe}^{-1/2}S_\tau +V^L_X&=-\mathcal{A}(S)E(\rho ), \end{align}

(5.10d) \begin{align}-\delta\textit{Pe}^{-1/2} S_\tau +V^G_X&=\mathcal{A}(S)E(\rho ), \end{align}

(5.10e) \begin{align}\delta \nu\textit{Pe}^{-1/2}(1-S)\rho _\tau +\nu V^G\rho _X&=\big(D^{\textit{eff}}(S)\rho _X\big)_X+\left (1-\nu \rho \right )\mathcal{A}(S)E(\rho ). \end{align}

We see that (5.10c )–(5.10e ) are all quasi-steady to leading order in $\textit{Pe}^{-1/2}\ll 1$ , with no time derivatives.

Furthermore, since we assume that both $\textit{Ca}\ll \xi\textit{Pe}$ and $\textit{Ca}\ll \sqrt {\textit{Pe}}\ll\textit{Pe}$ , we see from (5.10a )–(5.10b ) that both $p^L_X=0$ and $p^G_X=0$ to leading order. Thus, the capillary pressure $p_c(S)=p^G-p^L$ is spatially uniform and so the saturation $S$ is spatially uniform to leading order (although of course it varies with time $\tau$ ). Since we also found $S$ to be independent of $z$ in the outer region, the saturation is spatially uniform throughout the porous medium and evolves according to the temporal ODE (5.7).

Within the boundary layer, the model reduces at leading order to the ODE system

(5.11a) \begin{align}V^L_X&=-\mathcal{A}(S)E(\rho ), \end{align}

(5.11b) \begin{align}V^G_X&=\mathcal{A}(S)E(\rho ), \end{align}

(5.11c) \begin{align}\nu V^G\rho _X&=\big(D^{\textit{eff}}(S)\rho _X\big)_X+\left (1-\nu \rho \right )\mathcal{A}(S)E(\rho ), \end{align}

which must be solved at each time $\tau$ or, equivalently, for each saturation value $S(\tau )$ . We note that (5.11b )–(5.11c ) form a closed system for $\rho$ and $V^G$ , upon solution of which $V^L$ is given by (5.11a ).

By combining (5.11a )–(5.11b ) to obtain an expression of total conservation of mass, we see that $V^L+V^G$ is independent of $X$ and, by both matching with the outer region as $X\rightarrow -\infty$ and applying the boundary conditions at $X=0$ , we find that

(5.12) \begin{equation} V^L+V^G=V^{L,BL}=V^{G,\textit{out}}, \end{equation}

where $V^{G,\textit{out}}$ is the gas flux out of the top of the porous medium, at $X=0$ or $z=1$ . The ODE (5.7) for $S$ may therefore equivalently be written as

(5.13) \begin{equation} S_\tau =-V^{G,\textit{out}}. \end{equation}

Similarly, we note from (5.11b ) that the net evaporation rate (4.4) is

(5.14) \begin{equation} \text{ net evaporation rate }=\int _{z=0}^1\mathcal{A}(S)E(\rho )\,\mathrm{d}z=\frac {1}{\sqrt {\textit{Pe}}}\int _{X=-\infty }^0V^G_X\,\mathrm{d}X=\frac {V^{G,\textit{out}}}{\sqrt {\textit{Pe}}}. \end{equation}

To solve (5.11b )–(5.11c ) for $\rho$ and $V^G$ (and so find $V^{G,\textit{out}}$ ), we must specify a functional form for $E(\rho )$ . In the following subsection, we use a linear dependence and find an explicit solution in that case. For general $E(\rho )$ , we may combine (5.11b )–(5.11c ) to find a first integral

(5.15) \begin{equation} V^G\kern-2pt\left (1-\nu \rho \right )+D^{\textit{eff}}\rho _X=0, \end{equation}

where we have fixed the integration constant to be zero by matching with the outer region, in which $V^G$ and $\rho _X$ are both zero. Thus, our third-order problem (5.11b )–(5.11c ) reduces to the autonomous, first-order system

(5.16a) \begin{align} D^{\textit{eff}}\rho _X=-\left (1-\nu \rho \right )V^G, \quad \quad V^G_X =\mathcal{A}E(\rho ), \end{align}

with

(5.16b) \begin{equation} \rho \rightarrow 1, \quad V^G\rightarrow 0, \qquad \text{ as }X\rightarrow -\infty . \end{equation}

Furthermore, with the change of variables

(5.17) \begin{equation} X=\sqrt {\frac {D^{\textit{eff}}(S)}{\mathcal{A}(S)}}Y, \qquad V^G=\sqrt {\mathcal{A}(S)D^{\textit{eff}}(S)}\hat {V}^G, \end{equation}

the $S$ dependence is removed from (5.16), leaving the system

(5.18a) \begin{align} \rho _Y&=-\left (1-\nu \rho \right )\hat {V}^G, \end{align}

(5.18b) \begin{align} \hat {V}^G_Y&=E(\rho ), \end{align}

with

(5.18c) \begin{equation} \rho \rightarrow 1, \quad \hat {V}^G\rightarrow 0, \qquad \text{ as }Y\rightarrow -\infty . \end{equation}

Figure 4. Sketched phase plane portrait for the autonomous system (5.18) in the boundary layer during stage 1. For the required matching with the outer region, the solution trajectory must originate at the critical point $\rho =1$ , $\hat {V}^G=0$ , shown in red.

We can understand (5.18) to some extent for generic $E$ with a phase-plane analysis. The phase plane is sketched in figure 4. In particular, the critical point of the system is at $\rho =1$ , $\hat {V}^G=0$ , which we classify as a saddle, so long as $E(\rho )$ is decreasing in $\rho$ at the saturation point $\rho =1$ . Physical solutions must satisfy $\rho \in [0,1]$ and $\hat {V}^G\gt 0$ . In this range, we see from (5.18a ) that $\rho$ is monotone decreasing, while $\hat {V}^G$ is monotone increasing, since $E(\rho )\geqslant0$ . In particular, at $\rho =0$ , we must have $\rho _Y$ non-zero for $\hat {V}^G\gt 0$ , and so we expect solution trajectories to intersect the $\rho =0$ axis at finite $Y$ and finite $\hat {V}^G\gt 0$ . For the matching conditions (5.18c ) as $Y\rightarrow -\infty$ , the solution trajectory of interest is the unique trajectory that approaches the critical point as $Y\rightarrow -\infty$ , marked in red in figure 4.

Over time, we expect $S$ to decrease. The effective vapour diffusivity $D^{\textit{eff}}$ is likely to be decreasing in $S$ . The interfacial area fraction $\mathcal{A}$ might have non-montone variation with $S$ , but it is reasonable to assume this is increasing with $S$ , especially as $S$ becomes small. Thus, over time, we might expect variation in the width $\sqrt {D^{\textit{eff}}/\mathcal{A}}$ of the boundary layer (from (5.17)) and, especially as $S$ begins to become small, we expect the boundary-layer width to increase.

5.1.3. Solution when $E$ is linear in $\rho$

As discussed in § 2.2, many models for the evaporation rate $E$ are linear in the vapour density $\rho$ . In this section, we specify the Hertz–Knudsen relation (2.4) so that, in our dimensionless variables,

(5.19) \begin{equation} E(\rho )=1-\rho . \end{equation}

With this (commonly used) form for $E$ , we can make analytical progress in solving for boundary-layer processes and, thus, for the evolution of the saturation, $S$ , via (5.13).

With the linear form (5.19), the system (5.18) becomes

(5.20a) \begin{align} \rho _Y &= -\kern-0.2pt \left (1-\nu \rho \right )\hat {V}^G, \end{align}

(5.20b) \begin{align} \hat {V}^G_Y&=1-\rho . \end{align}

Combining these, we find the single, separable equation

(5.21) \begin{equation} \frac {\mathrm{d}\rho }{\mathrm{d}\hat {V}^G}=\frac {\rho _Y}{\hat {V}^G_Y}=-\frac {1-\nu \rho }{1-\rho }\hat {V}^G, \end{equation}

with solution

(5.22) \begin{equation} \left(\hat {V}^G\right)^2=\frac {2}{\nu }\left (1-\rho +\left (\frac {1}{\nu }-1\right )\log \left (\frac {1-\nu }{1-\nu \rho }\right )\right ), \end{equation}

where we have applied the matching conditions with the outer problem to fix the integration constant, specifying that $\rho \rightarrow 1$ and $\hat {V}^G\rightarrow 0$ as $Y\rightarrow -\infty$ .

We can now use this relation (5.22) between $\hat {V}^G$ and $\rho$ to find a single, first-order ODE for $\rho$ , by eliminating $\hat {V}^G$ in (5.20a ), say. Indeed, doing so results in the implicit solution for $\rho$

(5.23) \begin{equation} \int _{\rho _\infty }^\rho \frac {1}{(1-\nu \tilde {\rho })\sqrt {1-\tilde {\rho }+\left (1/\nu -1\right )\log \left ((1-\nu )/(1-\nu \tilde {\rho })\right )}}\,\mathrm{d}\tilde {\rho }=- \sqrt {\frac {2}{\nu }}Y, \end{equation}

if we impose the boundary condition $\rho =\rho _\infty$ at $Y=0$ .

However, we already gain significant information from this first integral (5.22). Applying the boundary condition for $\rho$ at the top of the porous medium ( $z=1$ or $Y=0$ ) in (5.22) gives us precisely the value $V^{G,\textit{out}}=V^G|_{X=0}=\sqrt {\mathcal{A}D^{\textit{eff}}}\hat {V}^G|_{Y=0}$ needed to understand the overall evaporation rate $S_\tau =-V^{G,\textit{out}}$ as in (5.13). We can do this for either choice of boundary condition discussed in § 2.3, and now look at each in turn.

Firstly, if we impose the Dirichlet condition for the vapour density $\rho =\rho _\infty$ at the top $Y=0$ , then immediately from (5.22), we see that

(5.24) \begin{align} S_\tau =-V^{G,\textit{out}}&=-\sqrt {\frac {2\mathcal{A}(S)D^{\textit{eff}}(S)}{\nu }}\sqrt {1-\rho _\infty +\left (\frac {1}{\nu }-1\right )\log \left (\frac {1-\nu }{1-\nu \rho _\infty }\right )}. \end{align}

Although we have worked with $\nu$ throughout as an order 1 parameter, it is helpful to note at this point that $\nu =O(10^{-2})$ is actually fairly small. In this case, (5.24) becomes the simpler expression

(5.25) \begin{align} S_\tau &=-V^{G,\textit{out}}=-\sqrt {\mathcal{A}(S)D^{\textit{eff}}(S)}\left (1-\rho _\infty +O(\nu )\right ) \quad \text{ if }\nu \ll 1. \end{align}

In the case where $\nu \ll 1$ , we may also integrate (5.23) to find an explicit solution for $\rho$ and, therefore, $V^G=\sqrt {\mathcal{A}D^{\textit{eff}}}\hat {V}^G$ via (5.22) in the boundary layer, which is simply $V^G=\sqrt {\mathcal{A}D^{\textit{eff}}}(1-\rho )+O(\sqrt {\nu })$ in this limit. We find

(5.26a) \begin{align}& \rho =1-(1-\rho _\infty )\exp \left (\sqrt {\frac {\mathcal{A}}{D^{\textit{eff}}}}X\right ), \qquad V^G=\sqrt {\mathcal{A}D^{\textit{eff}}}(1-\rho _\infty )\exp \left (\sqrt {\frac {\mathcal{A}}{D^{\textit{eff}}}}X\right ), \end{align}

(5.26b) \begin{align} &\qquad\qquad\qquad V^L=\sqrt {\mathcal{A}D^{\textit{eff}}}(1-\rho _\infty )\left (1-\exp \left (\sqrt {\frac {\mathcal{A}}{D^{\textit{eff}}}}X\right )\right ), \end{align}

where we have used (5.12) to fix $V^L$ and returned to spatial variable $X=\sqrt {(D^{\textit{eff}}/\mathcal{A})}\,Y$ .

We saw in § 5.1.1 that $V^L$ is linear in $z$ in the outer region, while $\rho =1$ and $V^G=0$ are constant. We form the leading-order composite solution expansions, valid over the entire domain (outer and boundary layer):

(5.27a) \begin{align}\rho &=1-(1-\rho _\infty )\exp \left (-\sqrt {\frac {\mathcal{A}\textit{Pe}}{D^{\textit{eff}}}}(1-z)\right ), \end{align}

(5.27b) \begin{align}U^G&=\sqrt {\frac {\mathcal{A}D^{\textit{eff}}}{\textit{Pe}}}(1-\rho _\infty )\exp \left (-\sqrt {\frac {\mathcal{A}\textit{Pe}}{D^{\textit{eff}}}}(1-z)\right ), \end{align}

(5.27c) \begin{align}U^L&=\sqrt {\frac {\mathcal{A}D^{\textit{eff}}}{\textit{Pe}}}(1-\rho _\infty )\left (z-\exp \left (-\sqrt {\frac {\mathcal{A}\textit{Pe}}{D^{\textit{eff}}}}(1-z)\right )\right ). \end{align}

This completes our analysis of stage-1 drying with the Dirichlet boundary condition; we will consider the flux condition shortly. First, however, we note from (5.24) that if the product $\mathcal{A}(S)D^{\textit{eff}}(S)$ is constant (independent of $S$ ), then the net evaporation rate is constant, as expected in this ‘constant rate period’ of the drying. Various experimental measurements show that the net evaporation rate is usually not strictly constant during stage 1, but could either rise or fall (Wu & Prat Reference Wu and Prat2022). We would expect $D^{\textit{eff}}(S)$ to be a decreasing function of $S$ : with less liquid in the pore space, the effective rate of the diffusion of vapour is increased. The dependence of the liquid–gas surface area $\mathcal{A}(S)$ will depend heavily on the pore geometry, but we might expect that $\mathcal{A}$ is an increasing function of $S$ , at least for fairly small $S$ . Our analysis suggests that it is the interplay between these effects that determine whether the evaporation rate, given by (5.24) and proportional to $\sqrt {\mathcal{A}D^{\textit{eff}}}$ , rises or falls over stage 1.

In figure 5, we show the time variation of the net evaporation rate for three different pairs of $\mathcal{A}$ and $D^{\textit{eff}}$ , showing a rising, falling and truly constant rate. The drying rate during stage 1 is very well predicted by (5.24) (dashed green curves) for all cases. The time at which stage 1 ends is not very well captured by (5.24), although this improves for larger $\textit{Pe}$ and smaller $\textit{Ca}$ (this is discussed further in § 5.2).

Figure 5. Stage-1 drying. Numerical solutions (solid curves) of (3.3) compared with the large- $\textit{Pe}$ approximation (dashed green) where $S$ is the solution of the ODE (5.24), and the analytical expressions for $\rho$ , $U^L$ and $U^G$ are (5.27). (a) Comparison of numerical (solid colour) and asymptotic (green dashed) net evaporation rates (4.4) for various $\mathcal{A}$ and $D^{\textit{eff}}$ (with $\textit{Ca}=1$ ). Crosses mark the point where $S$ computed from (5.24) reaches $\textit{Ca}\textit{Pe}^{-1/2}$ . (b–d) Numerical and approximate profiles of $S$ , $\rho$ , $U^L$ and $U^G$ at uniform time steps through stage 1, with $\mathcal{A}=0.5$ , $D^{\textit{eff}}=1$ constant (with $\textit{Ca}=0.1$ ). In panel (d), the arrow shows direction of increasing time for the numerical $U^L$ profile. Throughout, we take $\textit{Pe}=250$ , $\nu =0.2$ , $\xi =0.5$ and all other parameters as in table 2.

In figure 5(b–d), we compare the numerically computed profiles of $S$ , $\rho$ , $U^L$ and $U^G$ against the solution of the asymptotic reduction (5.24) and the corresponding composites (5.27). We take $\mathcal{A}$ and $D^{\textit{eff}}$ constant, so that the approximate $\rho$ , $U^L$ and $U^G$ profiles (5.27) are independent of $S$ , and so of time. The vapour density $\rho$ and gas velocity $U^G$ are very well approximated by (5.27), as are $S$ and $U^L$ for the majority of stage 1. However, as the saturation $S$ becomes smaller, it is no longer spatially uniform and the liquid flux $U^L$ is reduced. This is due to the viscous drag force, neglected in our stage-1 approximation, becoming important and hindering the capillary flow. We discuss this further in § 5.2.

The above mentioned analysis is for the case of a Dirichlet boundary condition for $\rho$ at the surface. In the alternative case, the boundary condition for $\rho$ at the surface of the porous medium is the Newton-cooling mass flux condition

(5.28) \begin{equation} \nu U^G\rho -\textit{Pe}^{-1}D^{\textit{eff}}\rho _z=m(\rho -\rho _\infty ) \qquad \text{ at }z=1. \end{equation}

In our boundary-layer variables, this becomes

(5.29) \begin{equation} \nu V^G\rho -D^{\textit{eff}}\rho _X=M(\rho -\rho _\infty ) \qquad \text{ at }X=0, \end{equation}

where $M=m\sqrt {\textit{Pe}}$ . If $M\gg 1$ , then we regain the Dirichlet condition $\rho =\rho _\infty$ . The case where $M\ll 1$ corresponds to vapour transport through the air being much slower than vapour diffusion through the porous medium, and thus mass transfer external to the material directly sets the evaporation rate (we do not consider this case further here). If $M=O(1)$ , then we note from the first integral (5.15) that this boundary condition (5.29) may be written as

(5.30) \begin{equation} V^G= M(\rho -\rho _\infty ) \qquad \text{ at }X=0. \end{equation}

Evaluating (5.22) at $X=Y=0$ and applying (5.30), we obtain an algebraic equation satisfied by $V^{G,\textit{out}}=\sqrt {\mathcal{A}D^{\textit{eff}}}\hat {V}^G|_{Y=0}$ , namely

(5.31) \begin{equation} \frac {(V^{G,\textit{out}})^2}{2\mathcal{A}D^{\textit{eff}}}=\frac {1}{\nu }\left (1-\rho _\infty -\frac {V^{G,\textit{out}}}{M}+\left (\frac {1}{\nu }-1\right )\log \left (\frac {1-\nu }{1-\nu \rho _\infty -\nu V^{G,\textit{out}}\kern-2pt/\kern-2pt M}\right )\right ). \end{equation}

In the case where $\nu \ll 1$ , we can solve for the leading-order $V^{G,\textit{out}}$ , finding that

(5.32) \begin{equation} V^{G,\textit{out}}=\sqrt {\mathcal{A}D^{\textit{eff}}}\left (\frac {1-\rho _\infty }{1+\sqrt {\mathcal{A}D^{\textit{eff}}}{\kern-2pt}/{\kern-2pt}M}\right )+O(\nu ) \qquad \text{ if }\nu \ll 1. \end{equation}

As in the case of the Dirichlet boundary condition, the solution $V^{G,\textit{out}}(S)$ of (5.31) then gives the net rate of evaporation, according to (5.14). We plot solutions $V^{G,\textit{out}}$ of (5.31) in figure 6 and, in particular, see that for the physical value $\nu =0.02$ , there is excellent agreement with the $\nu \ll 1$ approximation (5.32).

Figure 6. Net evaporation rate $V^{G,\textit{out}}$ during stage 1 as a function of the mass-transfer rate through the surrounding air, $M$ . Solid curves are the solution of (5.31) and dotted curves are in the large- $M$ limit where we regain (5.24), the case of a Dirichlet boundary condition. We use $\rho =0.1$ , $\mathcal{A}D^{\textit{eff}}=1$ .

5.2. Transition from stage 1 to stage 2

Consistent with experimental and numerical results in the literature (e.g. Pel et al. Reference Pel, Brocken and Kopinga1996), all our numerical simulations of the full model (3.3) shown so far exhibit a sudden end to stage 1 and a rapid transition to stage 2 (analysed in § 5.3). Stage 1 is generally understood to end when the capillary forces are no longer sufficiently strong relative to the viscous drag to draw liquid to the surface of the porous medium. We saw in figure 2 that the transition occurs when first $S=0$ at $z=1$ . In the literature, the transition is observed to occur when the average saturation dips below a critical value (Gupta et al. Reference Gupta, Huinink, Prat, Pel and Kopinga2014; Hall & Hoff Reference Hall and Hoff2021), often seen to be around 0.1–0.3, but strongly dependent on the properties of the porous medium and varying with the evaporation rate (Gupta et al. Reference Gupta, Huinink, Prat, Pel and Kopinga2014).

From our asymptotic analysis, we are in a position to quantify the scaling for this critical saturation, and its dependence on the system parameters and constitutive laws. We note that the reduced model for stage 1 given in § 5.1 does not directly capture the transition, as the solution $S$ is spatially uniform. However, we see that our reduced stage-1 model must break down if $k^L(S)p_c(S)$ approaches zero. Specifically, if at some saturation $S$ we have a small $k^L(S)p_c(S)=O(\textit{Ca}\textit{Pe}^{-1/2})\ll 1$ , then the viscous drag becomes large enough to balance the capillary forces in (3.3a ). Capillary forces then no longer dominate over the order-one domain and so we no longer expect a uniform saturation profile through the porous medium. This slowing of the liquid flux by viscous drag is what eventually leads to a zero liquid saturation at the surface. The balance of capillary and viscous forces in (3.3a ), with the stage-1 scaling $U^L=O(\textit{Pe}^{-1/2})$ , gives an estimate for the critical saturation: viscous forces will become non-negligible, heralding the end of stage 1, when $S=O(S_{\textit{crit}})$ defined by

(5.33) \begin{equation} k^L(S_{\textit{crit}})p_c(S_{\textit{crit}})=\textit{Ca}\textit{Pe}^{-1/2}. \end{equation}

If, for instance, $k^L(S)p_c(S)\sim S^n$ as $S\rightarrow 0$ for some $n\gt 0$ , then the critical saturation is estimated as $S_{\textit{crit}}\approx\textit{Ca}^{1/n}\textit{Pe}^{-1/(2n)}$ . Furthermore, a lower bound on the transition time between stage 1 and stage 2, $t_*$ , is given by the time at which the (spatially uniform) solution of the reduced model (5.24) satisfies $k^L(S)p_c(S)=\textit{Ca}\textit{Pe}^{-1/2}$ , or $S=S^{{crit}}$ . By integrating (5.24), we find this (under) estimate for the transition time $t_*$ to be

(5.34a) \begin{equation} t_*^{\textit{lower}}=\sqrt {\textit{Pe}}\sqrt {\frac {\nu }{2\left (1-\rho _\infty +\left (\dfrac {1}{\nu }-1\right )\log \left (\dfrac {1-\nu }{1-\nu \rho _\infty }\right )\right )}}\int _{S_{\textit{crit}}}^{S_0}\frac {1}{\sqrt {\mathcal{A}(S)D^{\textit{eff}}(S)}}\,\mathrm{d}S. \end{equation}

This lower bound (5.34a ) on $t_*$ is marked with a cross in figure 5(a), for each of the different functional forms for $\mathcal{A}(S)$ and $D^{\textit{eff}}(S)$ .

Conversely, as observed in figure 5(a), allowing our stage-1 reduced model (5.24) to run until $S=0$ may be used to find an upper bound on the time $t_*$ of the transition to stage 2. Specifically, by integrating (5.24), this upper bound is

(5.34b) \begin{equation} t_*^{\textit{upper}}=\sqrt {\textit{Pe}}\sqrt {\frac {\nu }{2\left (1-\rho _\infty +\left (\dfrac {1}{\nu }-1\right )\log \left (\dfrac {1-\nu }{1-\nu \rho _\infty }\right )\right )}}\int _{0}^{S_0}\frac {1}{\sqrt {\mathcal{A}(S)D^{\textit{eff}}(S)}}\,\mathrm{d}S. \end{equation}

The true transition time found in the full numerical solution appears to lie fairly centrally in the interval $[t_*^{\textit{lower}},t_*^{\textit{upper}}]$ for each case shown in figure 5.

Physically, we expect $k^L(S)$ to become small for sufficiently small $S$ , while we expect that $p_c$ should grow as $S$ becomes small. It is very interesting that we require the product $k^L(S)p_c(S)$ to become small for stage 1 to end, as this is not true in general: for instance, if we choose $k^L= S$ and $p_c\sim 1/S$ as $S\rightarrow 0$ , which seem to have the correct behaviour described in § 2.2, then the product $k^{L}\kern-2pt p_{c}\rightarrow 1$ as $S\rightarrow 0$ . As shown in figure 7, there is no stage 2 in this drying scenario: $S$ remains uniform at all times and stage 1 lasts until the drying is complete, so that our stage-1 reduced model (5.24) remains valid for the entire drying process. Similar behaviour to this was seen experimentally by Gupta et al. (Reference Gupta, Huinink, Prat, Pel and Kopinga2014). They observed that the crystallisation of salt within the pore space of building materials during drying impeded the onset of stage 2, so that constant-rate-period drying, with uniform saturation profiles, was seen to last throughout the entire drying process. Using a simplified continuum model, Gupta et al. (Reference Gupta, Huinink, Prat, Pel and Kopinga2014) gave a qualitatively similar (although less general) method to compute the critical saturation. In light of our analysis, we might interpret the results of Gupta et al. (Reference Gupta, Huinink, Prat, Pel and Kopinga2014) to mean that the presence of salt crystals increased the permeability and/or capillary pressure of the material to the extent that $k^{L}\kern-2pt p_{c}$ did not become sufficiently small to trigger stage 2.

Since the qualitative drying behaviour (the existence of a stage 2) depends on the behaviour of these constitutive laws $k^L(S)$ and $p_c(S)$ for small $S$ , it is critical that these are chosen appropriately for the porous medium and drying scenario at hand, if we hope to capture the real drying behaviour.

Figure 7. Drying with no stage 2 when $k^L(S)=S$ , $p_c(S)=1/S-S$ , so that $k^{L}\kern-2pt p_{c}\sim 1$ as $S\rightarrow 0$ . Numerical solutions of the full model (3.3) (solid curves) compared with the stage-1 model (5.24) (green dashed). Here, we take $\mathcal{A}=0.5$ and $D^{\textit{eff}}=1$ constant, $\textit{Pe}=100$ , $\textit{Ca}=0.1$ , and all other parameters as in table 2.

5.3. Stage 2: the ‘receding front period’

We now consider the behaviour during stage 2 of the drying, for $t\gt t_*$ , when we expect the evaporation rate to drop and a drying front to recede into the porous medium.

We look for a receding drying front at $z=h(t)$ , with a wet region ( $S\gt 0$ ) for $0\lt z\lt h$ and a dry region ( $S=0$ ) for $h\lt z\lt 1$ . A schematic diagram for stage 2 is shown in figure 8. Since $\textit{Pe}\gg 1$ , following the same argument as in § 5.1 for stage 1, the evaporation term dominates the vapour density equation (3.3e ) in the wet region $z\lt h$ , so that $\rho =1$ with exponentially small error in that region. We therefore expect that the evaporation will take place only in a thin transition layer at the front $z=h$ . Throughout this section, we take the local evaporation rate $E(\rho )=1-\rho$ to be linear in $\rho$ .

Figure 8. Schematic diagram showing the solution structure during stage 2, with a receding drying front at  $z=h(t)$ .

During stage 2, the saturation $S$ is small, since, as we have seen in § 5.2, the viscous drag must be too strong for the capillary forces to pull liquid fully to the surface of the porous medium (so that stage-1 drying stops). We will actually find that it is appropriate to take different scalings for $S$ in the wet region $z\lt h$ and in the transition layer near to $z=h$ . We use the saturation scalings

(5.35) \begin{equation} S=\eta \hat {S}, \quad \text{in }z\lt h, \qquad S=\bar {\eta } \bar {S}, \quad \text{in the transition layer}, \end{equation}

where $\bar {\eta }\ll \eta \ll 1$ which we specify later. We suppose that $\mathcal{A}\sim S^m$ and $k^{L}\kern-2pt p_{c}\sim S^n$ as $S\rightarrow 0$ (with both powers $n,m\gt 0$ ). Given this behaviour, we see that in the transition layer,

(5.36) \begin{equation} \mathcal{A}(S)=\bar {\eta }^m\bar {\mathcal{A}}(\bar {S}), \qquad k^L(S)p_c(S)=\bar {\eta }^n\bar {k}^L(\bar {S})\bar {p}_c(\bar {S}), \end{equation}

while in the wet region, $z\lt h$ ,

(5.37) \begin{equation} k^L(S)p_c(S)=\eta ^n\hat {k}^L(\hat {S})\hat {p}_c(\hat {S}). \end{equation}

We will assume in the subsequent analysis that $n\gt m$ , so that the capillary flow rate $k^{L}\kern-2pt p_{c}$ approaches zero at a faster rate as the saturation decreases than the liquid–gas interfacial area $\mathcal{A}$ .

To investigate the moving front behaviour, we rescale into a thin transition layer at $z=h$ , by making the change of variables

(5.38a) \begin{align}& t=t_*+\eta\textit{Pe} T, \qquad z=h(T)+\beta X, \qquad \rho =1-\beta \bar {\rho }, \qquad E=\beta \bar {\rho }, \end{align}

(5.38b) \begin{align}&\,\,\, U^G=\frac {1}{\textit{Pe}}W^G, \qquad U^L=\frac {1}{\textit{Pe}}W^L, \qquad p^G=\tilde {p}^G(h)+\frac {{\textit{Ca}}\beta }{\textit{Pe}\xi }\bar {p}^G. \end{align}

Here, the thickness of the transition layer is

(5.39) \begin{equation} \beta =\bar {\eta }^{-m/2}\textit{Pe}^{-1/2}\ll 1, \end{equation}

so that there is a balance of evaporation in the transition layer and (largely diffusive) transport of vapour out of the porous medium through the dry region. If $m\gt 0$ , this transition layer is wider than the $O(\textit{Pe}^{-1/2})$ boundary layer during stage 1, since $\bar {\eta }\ll 1$ . The variation of the vapour density from its saturation value 1, in (5.38), must be of order $\beta$ in the transition layer to match the vapour flux leaving the transition layer with that through dry outer region $z\in (h,1)$ . A balance in the liquid Darcy’s law (3.3a ) in the transition layer requires that

(5.40) \begin{equation} \bar {\eta }^n=\frac {\beta {\textit{Ca}}}{\textit{Pe}}. \end{equation}

The two constraints (5.39) and (5.40) fix the transition-layer width and saturation scaling to be

(5.41) \begin{equation} \beta =\textit{Ca}^{-\frac {m}{m+2n}}\textit{Pe}^{-\frac {(n-m)}{m+2n}}, \qquad \bar {\eta }=\textit{Ca}^{\frac {2}{m+2n}}\textit{Pe}^{-\frac {3}{m+2n}} \end{equation}

although we will retain the $\beta$ and $\bar {\eta }$ notation for simplicity. In particular, since $\textit{Pe}\gg 1$ , we have $\bar {\eta }\ll 1$ and the transition-layer thickness $\beta$ is small so long as $m\lt n$ , which we require to be the case for the transition-layer structure to be valid.

The time scale stated in (5.38) of the stage-2 drying is $\eta\textit{Pe}$ , which is determined by a balance in the liquid mass conservation equation (3.3c ) in the wet region. The velocity scales in (5.38) (which we note are smaller than stage 1) are determined by the liquid and gas mass equations, and the gas pressure scale comes from the gas Darcy law (3.3b ) (with $\tilde {p}^G(h)$ independent of $X$ and set by the gas transport in the outer region, as we will see later). We note that with these scalings, the net evaporation rate (4.4) is

(5.42) \begin{equation} \text{net evaporation rate}=\textit{Pe}^{-1}\int _{X\rightarrow -\infty }^\infty \hat {\mathcal{A}}\hat {\rho }\,\mathrm{d}X, \end{equation}

of order $\textit{Pe}^{-1}$ , and so slower than during stage 1 (where the net evaporation rate was $O(\textit{Pe}^{-1/2})$ ), as we expect.

Since the liquid velocity scale in the wet region $z\lt h$ must be the same as that in the transition layer (i.e. of order $\textit{Pe}^{-1}$ ), a distinguished limit of the model is when the viscous and capillary forces in Darcy’s law balance in the wet region, so that the saturation scaling here is given by

(5.43) \begin{equation} \eta =\left (\frac {\textit{Ca}}{\textit{Pe}}\right )^{1/n}. \end{equation}

At the end of stage 1, we saw in § 5.2 that the saturation scaled with $\textit{Ca}^{1/n}\textit{Pe}^{-1/2n}$ , which is larger than $\eta$ defined in (5.43). Thus, we might actually expect saturations larger than $O(\eta )$ (defined in (5.43)) at the start of stage 2. However, we will see that the large- $\hat {S}$ case is contained within the distinguished limit and thus may use the distinguished-limit saturation scaling (5.43) for stage 2. We also note that $\bar {\eta }=\beta ^{1/n}\eta$ , so that the transition-layer saturation is smaller than the wet-region saturation, as expected.

With these changes of variable (5.38), we find that to leading order in the transition layer, the model (3.3) becomes

(5.44a) \begin{align}W^L&= \bar {k}^L\bar {p}_{cX}, \end{align}

(5.44b) \begin{align}W^G&=-k^G\bar {p}^G_X, \end{align}

(5.44c) \begin{align}W^L_X&=-\bar {\mathcal{A}}\bar {\rho }, \end{align}

(5.44d) \begin{align}W^G_X&=\bar {\mathcal{A}}\bar {\rho }, \end{align}

(5.44e) \begin{align}0&=-\big (D^{\textit{eff}}_0\bar {\rho }_X\big)_X+ (1-\nu )\bar {\mathcal{A}}\bar {\rho }, \end{align}

where $D^{\textit{eff}}_0=D^{\textit{eff}}(0)$ . These transition-layer equations (5.44) are closed by matching with the surrounding $O(1)$ -length scale regions.

With the same time and velocity scales, and scaling $p^G=\textit{Ca}\textit{Pe}^{-1}\xi ^{-1}\tilde {p}^G$ , we see that to leading order in the dry region $h(T)\lt z\lt 1$ ,

(5.45) \begin{equation} W^G_z=0, \qquad W^G=-\tilde {p}^G_z, \qquad \nu W^G\rho _z=D^{\textit{eff}}_0\rho _{zz}. \end{equation}

At the surface $z=1$ , we impose the Dirichlet condition $\rho =\rho _\infty$ , noting that, at the slower evaporation rates during stage 2 than in stage 1, we expect that the alternative mass flux condition (3.5c ) (left) reduces to this Dirichlet form at leading order.

In the wet region $0\lt z\lt h(T)$ , since $E=0$ , we must have no gas flow $W^G=0$ and constant gas pressure $p^G=\tilde {p}^G(h)$ . The liquid equations are

(5.46) \begin{equation} \hat {S}_T+W^L_z=0, \qquad W^L=\hat {k}^L\hat {p}_{cz}, \end{equation}

which combine to a single diffusion equation for $S$ , namely

(5.47) \begin{equation} \hat {S}_T+\big(\hat {k}^L(\hat {S})\hat {p}_c^{\prime}(\hat {S})\hat {S}_z\big)_z=0. \end{equation}

5.3.1. Gas and vapour flow

We first study the gas and vapour flow through the dry region $h\lt z\lt 1$ and the transition layer. In this dry region, from (5.45), we see that $W^G=W^G_{\textit{dry}}(T)$ is spatially uniform. By combining the transition-layer equations (5.44d )–(5.44e ), we see that the quantity

(5.48) \begin{equation} (1-\nu )W^G-D^{\textit{eff}}_0\bar {\rho }_X=0 \end{equation}

is conserved over the transition layer, and equal to zero by matching with the wet region $z\lt h$ where there is no gas flow or vapour flux. Matching also with the dry region $z\gt h$ , we find that

(5.49) \begin{equation} (1-\nu )W^G_{\textit{dry}}=D^{\textit{eff}}_0\bar {\rho }_X\Big|_{X=0}=-D^{\textit{eff}}_0\rho _z\Big|_{z=h}. \end{equation}

Solving the vapour equation (5.45) (right) subject to the matching conditions (5.49) and $\rho =1$ at $z=h$ , and the boundary condition $\rho =\rho _\infty$ at $z=1$ , we obtain the solution in the dry region $z\in (h,1)$ :

(5.50a) \begin{equation} \rho =\frac {1}{\nu }\left (1-(1-\nu \rho _\infty )\left (\frac {1-\nu }{1-\nu \rho _\infty }\right )^{\frac {1-z}{1-h}}\right ), \qquad W^G_{\textit{dry}}=\frac {D^{\textit{eff}}_0\log \left (\dfrac {1-\nu \rho _\infty }{1-\nu }\right )}{\nu (1-h)}. \end{equation}

Then, solving (5.45) (centre), with $\tilde {p}^G=0$ at the surface $z=1$ , we see that the gas pressure through the dry region is

(5.50b) \begin{equation} \tilde {p}^G=W^G_{\textit{dry}}(1-z). \end{equation}

The constant gas pressure in the wet region is therefore

(5.51) \begin{equation} \tilde {p}^G(h)=W^G_{\textit{dry}}(1-h)=\frac {D^{\textit{eff}}_0\log \left (\dfrac {1-\nu \rho _\infty }{1-\nu }\right )}{\nu }, \end{equation}

independent of both space and time (as observed numerically in figure 2 d). As an aside, for small $\nu \ll 1$ , we note that (5.50) reduce to

(5.52) \begin{equation} \rho =\rho _\infty +(1-\rho _\infty )\frac {1-z}{1-h}+O(\nu ), \qquad W^G_{\textit{dry}}=\frac {D^{\textit{eff}}_0}{1-h}(1-\rho _\infty )+O(\nu ) \qquad \text{ as }\nu \rightarrow 0. \end{equation}

We have thus found the vapour and gas solutions in the order-1 wet and dry regions, in terms of the (as yet unknown) position, $h$ , of the drying front. Additionally, by integrating the transition-layer equation (5.44d ) over the entire transition layer and matching with the behaviour in the wet and dry regions on either side, we see that the net evaporation rate (5.42) over the transition layer – and so over the entire material since the transition layer is the only region with non-negligible evaporation – is

(5.53) \begin{equation} \text{(net evaporation rate)}\,\textit{Pe}=\int _{X\rightarrow -\infty }^\infty \hat {\mathcal{A}}\hat {\rho }\,\mathrm{d}X= W^G_{\textit{dry}}= \frac {D^{\textit{eff}}_0\log \left (\dfrac {1-\nu \rho _\infty }{1-\nu }\right )}{\nu (1-h)}. \end{equation}

Thus, given the position $h(T)$ of the receding front, the vapour transport through the dry region determines the net evaporation rate. However, $h(T)$ is still to be determined: to understand how the receding front moves, we must consider the liquid problem in the wet region and transition layer.

5.3.2. Liquid flow and the receding front

In the wet region $z\in (0,h(T))$ , the saturation $\hat {S}$ satisfies (5.47), with $W^L=\hat {k}^L\hat {p}_{cz}=0$ at $z=0$ , while in the transition layer, combining (5.44a ) and (5.44c ), we find

(5.54) \begin{equation} \big(\bar {k}^L(\bar {S})\bar {p}_c^{\prime}(\bar {S})\bar {S}_X\big)_X=-\bar {\mathcal{A}}\bar {\rho }, \end{equation}

with $\bar {S}\rightarrow 0$ and $W^L=\bar {k}^L\bar {p}_{cX}\rightarrow 0$ as $X\rightarrow \infty$ . The matching conditions between the wet and transition-layer regions are

(5.55) \begin{equation} \lim _{z\rightarrow h}(\hat {S})=0 \quad \text{and} \quad \lim _{z\rightarrow h}\big(\hat {k}^L\hat {p}_c^{\prime}\hat {S}_z\big)=\lim _{X\rightarrow -\infty }\big(\bar {k}^L\bar {p}_c^{\prime}\bar {S}_X\big), \end{equation}

equivalent to imposing continuity of saturation and liquid flux $W^L$ .

Integrating (5.54) in $X$ over the entire transition layer, using (5.53) and the fact that $\bar {S}\rightarrow 0$ and $W^L\rightarrow 0$ as $X\rightarrow \infty$ , we find that

(5.56) \begin{equation} W^G_{\textit{dry}}=\int _{X\rightarrow -\infty }^\infty \bar {\mathcal{A}}\bar {\rho }\,\mathrm{d}X=\big(\bar {k}^L(\bar {S})\bar {p}_c^{\prime}(\bar {S})\bar {S}_X\big)\big|_{X\rightarrow -\infty }. \end{equation}

Matching with the wet region using (5.55) (right), this gives

(5.57) \begin{equation} \big(\hat {k}^L(\hat {S})\hat {p}_c^{\prime}(\hat {S})\hat {S}_z\big)\big|_{z=h}=W^G_{\textit{dry}}, \end{equation}

which we can view as a boundary condition on the problem (5.47) in the wet region.

Using the expression (5.50a ) for $W^G$ , we therefore obtain a free-boundary problem on the wet region, namely

(5.58a) \begin{align} \hat {S}_T+\big(\hat {k}^L\hat {p}_c^{\prime}\hat {S}_z\big)_z&=0 &\text{ for }0\lt z\lt h(T),\, \end{align}

(5.58b) \begin{align} \hat {S}_z&=0 &\text{on } z=0,\qquad\quad\;\,\, \end{align}

(5.58c) \begin{align} \hat {k}^L\hat {p}_c^{\prime}\hat {S}_z&=\frac {D^{\textit{eff}}_0\log \left (\dfrac {1-\nu \rho _\infty }{1-\nu }\right )}{\nu (1-h)} & \text{on }z=h(T),\qquad \end{align}

(5.58d) \begin{align} \hat {S}&=0 &\text{on }z=h(T).\qquad \end{align}

This system (5.58) is a closed problem for $\hat {S}(z,T)$ and $h(T)$ . We note that $\hat {S}\rightarrow 0$ and so $\hat {k}^L\hat {p}_c^{\prime}\rightarrow 0$ as $z\rightarrow h$ . The flux remains bounded, but the gradient $\hat {S}_z$ must blow up here. The initialisation of the system (5.58) is not trivial, since the flux in (5.58c ) blows up as $h\rightarrow 1$ . We note that if $h$ is close to 1 and $\hat {S}$ is large, then the capillary flow term $\hat {k}^L\hat {p}_c^{\prime}\hat {S}_z$ dominates and liquid may easily be drawn to the evaporation front at the rate given by (5.58c ). In practice, we may therefore initialise (5.58) at some $h$ close to 1 and with a large $\hat {S}$ : over a fast initial transient, $\hat {S}$ reduces to be of order 1 and $h(T)$ will not begin to vary appreciably until this point. Thus, the system (5.58) is independent of the initialisation, up to the translation in time to account for the initial transient.

Given the solution of (5.58), we have the saturation, fluid velocity and pressure, and vapour density profiles throughout the order-1 length scale wet and dry regions, as well as the evolution of the receding drying front $h(T)$ and the net drying rate via (5.53).

Similarly to our analysis in stage 1, by integrating over the thin evaporating layer (here the transition layer moving with the receding front, rather than the stationary boundary layer at the surface as in stage 1), we have derived a reduced problem which gives the leading-order evaporation rate. We note that the dependence of the flux boundary condition (5.58c ) on $1-h$ precludes the possibility of a similarity solution and, therefore, expect that the reduced model (5.58) must be solved numerically. We will find numerical solutions for an example form of $\hat {k}^L\hat {p}_c$ in § 5.3.3.

Given the solution of (5.58) and thus the leading-order drying behaviour, we might then solve the transition-layer problem (5.44) to find the behaviour within the transition layer. We note that several of the equations (5.44) decouple, so that we need only solve the coupled diffusion equations for $\bar {S}$ and $\bar {\rho }$ , namely

(5.59a) \begin{align} \big(\bar {k}^L\bar {p}_c^{\prime}\bar {S}_X\big)_X&=-\bar {\mathcal{A}}(\bar {S})\bar {\rho }, \end{align}

(5.59b) \begin{align} D^{\textit{eff}}_0\bar {\rho }_{XX}&=(1-\nu )\bar {\mathcal{A}}(\bar {S})\bar {\rho }, \end{align}

closed by matching to the leading-order solution on either side as $X\rightarrow \pm \infty$ . Given the solution of this transition-layer problem (5.59), we may then straightforwardly compute $W^L$ , $W^G$ and $\bar {p}^G$ in the transition layer from (5.44). The solution of this transition-layer problem (5.59) will also motivate an $O(\beta )$ correction to the position $h$ of the receding front. We explore this for the example drying scenario in § 5.3.3.

In this section, we have derived a reduced model describing the motion of a drying front into the porous medium during stage-2 drying, in the case of $\textit{Pe}\gg 1$ . The model is valid so long as $\hat {S}$ remains of order 1. It breaks if the wet-region saturation $\hat {S}=O(\beta ^{1/n})$ shrinks to the same size as the saturation ( $\bar {S}$ , remaining order 1) in the transition layer. The model also ceases to be valid when the receding front $h$ approaches within an $O(\beta )$ length scale of the end of the domain at $z=0$ .

5.3.3. Example and comparison with numerical solution

To illustrate the analysis of the previous section, we consider the example constitutive laws

(5.60) \begin{equation} k^L(S)=S^3, \qquad p_c(S)=\frac {1}{S}-S, \qquad \mathcal{A}=\begin{cases}\text{constant} & \text{for }S\gt 0,\\ 0 & \text{for }S=0.\end{cases} \end{equation}

Thus, for this example, we have $m=0$ and $n=2$ , so that

(5.61) \begin{equation} \beta =\textit{Pe}^{-1/2}, \qquad \qquad \bar {\eta }=\textit{Ca}^{1/2}\textit{Pe}^{-3/4}, \qquad \eta =\textit{Ca}^{1/2}\textit{Pe}^{-1/2}. \end{equation}

From (5.60), we have $k^{L}\kern-2pt p_{c}^{\prime}=-S-S^3\approx -S$ since $S$ is small. The wet-region reduced model (5.58) is therefore

(5.62a) \begin{align}\hat {S}_T-\big(\hat {S}\hat {S}_z\big)_z=0 \qquad\text{ for }0&\lt z\lt h(T), \end{align}

(5.62b) \begin{align}\hat {S}_z=0 \qquad\text{ on } z&=0, \end{align}

(5.62c) \begin{align}-\hat {S}\hat {S}_z=\frac {D^{\textit{eff}}_0\log \left (\dfrac {1-\nu \rho _\infty }{1-\nu }\right )}{\nu (1-h)} \qquad\text{ on }z&=h(T), \end{align}

(5.62d) \begin{align}\hat {S}=0, \qquad\text{ on }z&=h(T). \end{align}

To solve (5.62) numerically for $\hat {S}$ and $h$ , we make a change of variables $z=h(T)y$ so that the model (5.62) holds on the fixed domain $y\in (0,1)$ . We discretise using the method of lines, upwinding the spatial derivatives in the artificial advection terms introduced by the change of variables and with central differences otherwise. The discretised system is solved in MATLAB using the same solver ode15s as for the full system. As discussed for (5.58), the reduced model (5.62) is difficult to initialise numerically since (5.62c ) becomes unbounded as $h\rightarrow 1$ . We initialised close to $h=1$ , but found that up to a shift in $T$ , the solution obtained was independent of the choice of initialisation.

Figure 9. Stage-2 drying. Numerical solutions of (3.3) (solid curves, colour indicates stage 1 versus stage 2, with parameters as in figure 2), are compared with the solutions of the leading-order Stage-2 reduced model (5.58) (magenta dashed curves). (a–b) Direct solutions $S$ and $h$ of (5.58); (c) liquid velocity, $U^L$ , for the wet region; (d–f) the evaporation rate, (5.53), the dry-region vapour density, and gas velocity (5.50), all computed using the solution $h$ from panel (b). We additionally show the composite solutions (5.69) (magenta dotted lines) in each of panels (a–c,e–f).

A solution computed in this way is shown in figure 9 (dashed magenta lines) and compared with solutions of the full model (3.3) (solid lines). We find excellent agreement up to the expected error, the width of the transition layer, which is of order $\beta =\textit{Pe}^{-1/2}\approx 0.03$ in this example. From figure 9(a), we see that the shapes of the $S(z)$ profiles are well captured through the wet region, until we reach the transition region where $S$ approaches zero. Since the reduced model (5.62) prescribes $S=0$ at the edge of the wet region, the position, $h$ , of the drying front is therefore under-predicted, as seen in figure 9(b), although still within the expected $O(\beta )$ error. Although the position, $h$ , is under-predicted by the reduced model, the net evaporation rate in figure 9(d) is captured very precisely (except for the very start and end of stage 2, when we do not expect our reduced model to be valid).

From figure 9(b), we observe (for both the full numerical solution and the reduced model solution) that the drying front velocity increases over time. This is because the saturation $S$ decreases over time, meaning both that the capillary flow rate decreases (so less liquid can be transported to the drying front) and also that there is simply less liquid at the front to vaporise. To meet the required vaporisation flux, $W^G_{\textit{dry}}$ , (which depends only on $h$ , not on $S$ ), the front must move faster. We emphasise that this increasing-speed drying front is contrary to the decreasing $t^{-1/2}$ front velocity we would find if there was no liquid capillary flow and the front motion was purely determined by the diffusion of vapour out of the material (as in e.g. Luckins et al. Reference Luckins, Breward, Griffiths and Please2023).

The liquid flux $W^L=-\hat {S}\hat {S}_z$ in the wet region captures the magnitude and gradient of the numerical velocity field well in figure 9(c), until the transition layer is reached. Similarly in figures 9(e)– 9(f), we see that the dry-region vapour density $\rho$ and gas-flux $W^G$ are well captured by (5.50a ), up to the expected $O(\beta )$ error.

Due to the simple, piecewise-constant form (5.60) of $\mathcal{A}$ in this example, we may easily solve the transition-layer equation (5.59b ) for $\bar {\rho }$ to find

(5.63) \begin{equation} \bar {\rho }=W^G_{\textit{dry}}\sqrt {\frac {1-\nu }{\bar {\mathcal{A}}D^{\textit{eff}}_0}}\exp \left (\sqrt {\frac {\bar {\mathcal{A}}(1-\nu )}{D^{\textit{eff}}_0}}X\right ) \quad \text{ for }\quad X\leqslant0, \end{equation}

having used the matching conditions $\bar {\rho }\rightarrow 0$ as $X\rightarrow -\infty$ and (5.49) at $X=0$ , and where $W^G_{\textit{dry}}(h(T))$ is given by (5.50a ). From (5.44c ) and (5.44d ), we then find that the Darcy velocities in the transition layer are, for $X\leqslant0$ ,

(5.64) \begin{align} W^G=W^G_{\textit{dry}}\exp \left (\sqrt {\frac {\bar {\mathcal{A}}(1-\nu )}{D^{\textit{eff}}_0}}X\right ), \qquad W^L=W^G_{\textit{dry}}\left (1-\exp \left (\sqrt {\frac {\bar {\mathcal{A}}(1-\nu )}{D^{\textit{eff}}_0}}X\right )\right ), \end{align}

and by integrating (5.59a ) and imposing that $\bar {S}=0$ at $X=0$ , we find the saturation is

(5.65) \begin{equation} \bar {S}=\sqrt {2W^G_{\textit{dry}}\left (-X-\sqrt {\frac {D^{\textit{eff}}_0}{\mathcal{A}(1-\nu )}}\left (1-\exp \left (\sqrt {\frac {\mathcal{A}(1-\nu )}{D^{\textit{eff}}_0}}X\right )\right )\right )} \quad \text{ for }X\leqslant0. \end{equation}

In particular, we see that the gradient $\bar {S}_X\gt 0$ at $X=0$ is finite at the drying front, unlike the outer (wet) problem solution $\hat {S}$ .

Figure 10. Schematic showing the matching between the outer (wet) and transition-layer regions. While $\sqrt {\beta }\bar {S}(X)$ and $\hat {S}(z)$ from (5.66) and (5.67) match at leading order, the constant translation of $h_1$ in $X$ improves the agreement. (This is not a formal $O(\beta )$ correction, which would require higher order terms in the expansions of $\bar {S}$ and $\hat {S}$ in each region).

We note that

(5.66) \begin{equation} \sqrt {\beta }\bar {S}\sim \sqrt {\beta }\sqrt {2W^G_{\textit{dry}}\left (-X-\sqrt {\frac {D^{\textit{eff}}_0}{\mathcal{A}(1-\nu )}}\right )} \quad \text{ as }X\rightarrow -\infty , \end{equation}

while from the boundary condition (5.62c ) on the wet problem,

(5.67) \begin{equation} \hat {S}\sim \sqrt {2W^G_{\textit{dry}}(h-z)} \quad \text{ as }z\rightarrow h. \end{equation}

While (by construction, since $z=h+\beta X$ ) these clearly match to leading order, there is a discrepancy arising from the constant term in (5.66). This is illustrated schematically in figure 10. This discrepancy therefore motivates a shift of the transition-layer variable, replacing $X$ by $X-h_1$ , where

(5.68) \begin{equation} h_1=\sqrt {\frac {D^{\textit{eff}}_0}{\mathcal{A}(1-\nu )}} \end{equation}

may be viewed as a (constant) $O(\beta )$ correction to the drying front position $h$ , and is chosen so that the limiting behaviours (5.66) and (5.67) fully agree (as shown in figure 10). We note that this is not a formal $O(\beta )$ correction to the interface position $h$ , which would require going to higher order in both the wet and transition-layer regions, and a more careful matching. Nevertheless, using this shift in $X$ , we may define composite solutions (by summing the inner (transition layer) and outer (wet region) solutions, and subtracting the common limit), valid throughout the domain $z\in (0,1)$ ,

(5.69a) \begin{align} h^{\textit{comp}}(T)&=h(T)+\beta h_1, \end{align}

(5.69b) \begin{align} S^{\textit{comp}}(z,T)&=\begin{cases}\eta \bigg (\hat {S}(z,T)+\sqrt {\beta }\sqrt {2W^G_{\textit{dry}}\left (-\dfrac {z-h}{\beta }+h_1\exp \left (\dfrac {z-h}{\beta h_1}-1\right )\right )}& \\ \def\lumina\hspace{5.5cm} -\sqrt {2W^G_{\textit{dry}}(h-z)}\bigg )&\kern-20pt\text{for }z\lt h^{\textit{comp}},\\ 0 &\kern-19.7pt\text{ for }z\gt h^{\textit{comp}},\end{cases} \end{align}

(5.69c) \begin{align} \rho ^{\textit{comp}}(z,T)&=\begin{cases} 1-\beta \dfrac {W^G_{\textit{dry}}}{\mathcal{A}h_1}\exp \left (\dfrac {z-h}{\beta h_1}\right ) &\text{ for }z\lt h^{\textit{comp}},\\ \frac {1}{\nu }\left (1-(1-\nu \rho _\infty )\left (\dfrac {1-\nu }{1-\nu \rho _\infty }\right )^{\frac {1-z}{1-h}}\right ) &\text{for }z\gt h^{\textit{comp}}, \end{cases} \end{align}

(5.69d) \begin{align} U^{L,{comp}}&=\begin{cases}\textit{Pe}^{-1}\left (-\hat {S}\hat {S}_z-W^G_{\textit{dry}}\exp \left (\dfrac {z-h}{\beta h_1}-1\right )\right ) &\text{ for }z\lt h^{\textit{comp}},\\ 0 &\text{ for }z\gt h^{\textit{comp}}, \end{cases} \end{align}

(5.69e) \begin{align} U^{G,{comp}}&=\begin{cases}\textit{Pe}^{-1}W^G_{\textit{dry}}\exp \left (\dfrac {z-h}{\beta h_1}-1\right ) &\text{ for }z\lt h^{\textit{comp}},\\\textit{Pe}^{-1}W^G_{\textit{dry}} &\text{ for }z\gt h^{\textit{comp}}. \end{cases} \end{align}

Here, $\hat {S}(z,T)$ and $h(T)$ are the solution of the wet-region problem (5.62), $W^G_{\textit{dry}}$ is as given by (5.50a ) (evaluated at $h(T)$ ), and we have used the definition (5.68) of $h_1$ to simplify the notation. (We have not applied the shift of $h_1$ for $\rho ^{\textit{comp}}$ in (5.69c ) for the composite solution to remain continuous.)

These composite solutions are plotted as magenta dotted curves in figure 9. We see in figure 9(b) that our ad hoc correction $h_1$ is good early in stage 2, but does not capture the correction fully at later times. Nevertheless, the composite solutions $S^{\textit{comp}}$ , $U^{L,{comp}}$ , $\rho ^{\textit{comp}}$ and $U^{G,{comp}}$ capture the profiles very well in figure 9(a,c,e,f), respectively, with the best agreement at earlier times (carrying through from the fact that $h^{\textit{comp}}$ is most accurate at earlier times).

Table 3. Summary of our $\textit{Pe}\gg 1$ analysis for the different stages of the drying process. Here, $k^L(S)p_c(S)\sim S^n$ as $S\rightarrow 0$ .

5.4. Summary of the drying stages

In table 3, we summarise our analysis of the limit $\textit{Pe}\gg 1$ . During stage 1, the saturation is uniform since capillary forces dominate, drawing liquid to the surface of the porous medium, with evaporation only occurring in a thin $O(\sqrt {\textit{Pe}})$ layer at the surface. The transition to stage 2 is reached when viscous drag limits the capillary flow of liquid; our (under) estimate is defined as the time where $S$ becomes small enough for capillary and viscous forces to balance. This occurs at the critical saturation scaling of $S=O(\textit{Pe}^{-1/2n}\textit{Ca}^{1/n})$ . Subsequently, in stage 2, we found a receding drying front at $z=h$ , with evaporation only taking place within a thin layer at this front. The net evaporation rate is determined by the vapour transport from the front to the surface, and capillary forces pull liquid to the front, slowing its retreat. Since capillary forces are reduced over time as $S$ becomes smaller, the receding front velocity increases over time.