2.1. Pore-scale model

In the pore space occupied by the vapour–gas mixture (behind the evaporating front, i.e. $y< h(x,t)$), we expect the Reynolds number to be small (Luckins et al. Reference Luckins, Breward, Griffiths and Please2023) and so we assume that the mixture satisfies the Stokes equations

(2.1a,b)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \quad -\boldsymbol{\nabla} p+\mu\nabla^2\boldsymbol{u}=\boldsymbol{0}, \end{equation}

where $\boldsymbol {u}$ is the mass-averaged velocity of the mixture, $p$ is the pressure and $\mu$ is the viscosity (assumed constant). The vapour contained within the mixture is advected with the flow, and also diffuses through the mixture with diffusivity $D_v$, and thus, the density of the vapour, $\rho _v$ [${\rm kg}\ {\rm m}^{-3}$], satisfies

(2.2)\begin{equation} (\rho_v)_t+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\rho_v=D_v\nabla^2\rho_v, \end{equation}

where the subscript $t$ denotes partial derivative. The overall density of the inert gas–vapour mixture, $\rho _G$, is given by

(2.3)\begin{equation} \rho_G=\rho_g+\rho_v. \end{equation}

Wherever the gas–vapour mixture meets the solid walls of the pore space we suppose there is no flux or slip of the gas–vapour mixture, and no flux vapour into the solid material, so that on the solid–liquid or dirt–liquid boundary, with normal $\boldsymbol {n}_s$,

(2.4a,b)\begin{equation} \boldsymbol{u}=\boldsymbol{0}, \quad \left(\rho_v\boldsymbol{u}-D_v\boldsymbol{\nabla}\rho_v\right)\boldsymbol{\cdot}\boldsymbol{n}_s=0. \end{equation}

In the liquid–dirt mixture in $y>h(x,t)$, we assume that there is no net flow of the mixture, and the suspended dirt and liquid diffuse against one another in an ideal mixture, due to Brownian motion. The suspended dirt volume fraction, $\theta$, therefore satisfies

(2.5)\begin{equation} \theta_t=D_d\nabla^2\theta, \end{equation}

where $D_d$ is the diffusivity of suspended dirt in liquid. As discussed in Luckins et al. (Reference Luckins, Breward, Griffiths and Please2023), the assumption that the liquid does not flow means that capillary effects are neglected from the model.

At the solid walls of the pore space, we suppose that the suspended dirt can deposit onto the solid microstructure, forming a layer that may also then be eroded away. We suppose that the deposited layer has a dirt volume fraction $\theta _*$ which is the packing volume fraction of the dirt particles. We expect this to be close to one, as only a small amount of liquid (volume fraction $1-\theta _*$) is trapped within the deposited dirt layer. Conservation of dirt across the interface is given by

(2.6)\begin{equation} \theta V_n +D_d\boldsymbol{\nabla}\theta\boldsymbol{\cdot}\boldsymbol{n}_s = \theta_*V_n, \end{equation}

where $V_n$ is the normal velocity of the depositing/eroding interface. We note that, in order that there is no flow generated at the depositing interface, we assume that the dirt and liquid have the same mass density, so that the total mixture density is the same on either side of the depositing/eroding interface, while the dirt and liquid fractions can jump (see, for instance, Geng, Kamilova & Luckins Reference Geng, Kamilova and Luckins2023).

We suppose that the dirt is deposited at a rate dependent on the local suspended dirt volume fraction, while the erosion rate depends on the (constant) packing volume fraction $\theta _*$. (If there was a flow of the fluid, we might extend this model and allow the erosion rate to depend on the local shear stress.) Thus, we prescribe

(2.7)\begin{equation} V_n=k_+\theta -k_-\theta_*, \end{equation}

where the constants $k_\pm$ have units ${\rm m}\ {\rm s}^{-1}$. This type of law-of-mass-action deposition rate, in which the deposition rate is linear in the quantity of suspended dirt, is common in the phenomenological bed-filtration literature (Zamani & Maini Reference Zamani and Maini2009; Dressaire & Sauret Reference Dressaire and Sauret2017), and is also used by Ji & Sanaei (Reference Ji and Sanaei2023) as a model for adsorption of particles onto the deposit layer.

At the evaporating interface $y=h(x,t)$, we suppose that the inert gas and the dirt do not pass through the interface, while liquid turns into vapour. We thus impose conservation of mass of each of the liquid/vapour, gas, and suspended dirt, namely

(2.8a)$$\begin{gather} -\rho_l(1-\theta)V_m+\rho_dD_d\boldsymbol{\nabla}\theta\boldsymbol{\cdot}\boldsymbol{m}= \rho_v\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{m}-V_m\right)-D_v\boldsymbol{\nabla}\rho_v\boldsymbol{\cdot}\boldsymbol{m}, \end{gather}$$

(2.8b)$$\begin{gather}0=\rho_g\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{m}-V_m\right)+D_v\boldsymbol{\nabla}\rho_v\boldsymbol{\cdot}\boldsymbol{m}, \end{gather}$$

(2.8c)$$\begin{gather}-\rho_d\theta V_m-\rho_dD_d\boldsymbol{\nabla}\theta\boldsymbol{\cdot}\boldsymbol{m}=0, \end{gather}$$

where $\rho _l$ and $\rho _d$ are the densities of pure liquid and dirt, respectively. Combining these, we derive the more helpful form

(2.9a)$$\begin{gather} -\rho_LV_m=\rho_G\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{m}-V_m\right), \end{gather}$$

(2.9b)$$\begin{gather}-\rho_LV_m=\rho_v\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{m}-V_m\right)-D_v \boldsymbol{\nabla}\rho_v\boldsymbol{\cdot}\boldsymbol{m}, \end{gather}$$

(2.9c)$$\begin{gather}\theta V_m+D_d\boldsymbol{\nabla}\theta\boldsymbol{\cdot}\boldsymbol{m}=0, \end{gather}$$

interpretable as a condition on each of $\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {m}$, $\rho _v$ and $\theta$, respectively, where $\rho _L=\rho _l(1-\theta )+\rho _d\theta$ is the (assumed constant) liquid–dirt mixture density. The normal velocity of the interface and unit normal to the interface are given by

(2.10a,b)\begin{equation} V_m=\frac{h_t}{\sqrt{1+h_x^2}}, \quad \boldsymbol{m}=\frac{(h_x,-1)}{\sqrt{1+h_x^2}}, \end{equation}

where subscripts $t$ and $x$ denote partial derivatives. In addition to (2.9), we also impose a no-slip condition for the gas-mixture velocity

(2.11)\begin{equation} u+vh_x=0 \quad \text{on }y=h(x,t). \end{equation}

Finally, we must also incorporate a condition that describes the chemistry governing the evaporation. In Luckins et al. (Reference Luckins, Breward, Griffiths and Please2023) the effect of different chemistry conditions were considered, and these were shown to affect the form of macroscale boundary conditions derived through a homogenisation analysis. For simplicity, we assume that the liquid and vapour are in chemical equilibrium at the evaporating interface. The chemical potential on the liquid side of the interface is dependent on the amount of liquid ($1-\theta$) at the interface, while the chemical potential on the gas-mixture side depends on the density of vapour, $\rho _v$, at the interface. In general, we may express this chemical equilibrium as

(2.12)\begin{equation} \rho_v=\rho^*f(\theta) \quad \text{on }y=h(x,t), \end{equation}

where $\rho ^*$ is the (constant) saturation vapour density when $\theta =0$ and there is no suspended dirt, and the function $f(\theta )$ captures the dirt dependence of the saturation vapour density. (We note that therefore $f(0)=1$.) The presence of particles at the interface are expected to hinder the vaporisation; in both Ji & Sanaei (Reference Ji and Sanaei2023) and Karapetsas et al. (Reference Karapetsas, Sahu and Matar2016) the evaporative flux is modelled as decreasing with increased particles on the fluid surface. We keep $f(\theta )$ general as far as possible, and in § 5 we investigate the effect of different functional dependencies $f(\theta )$ on the drying rate. However, in our numerical simulations we use the simple linear form

(2.13)\begin{equation} f(\theta)=1-\theta \end{equation}

to capture the effect of the dirt inhibiting vaporisation. We choose this form so that the saturation vapour density scales with the amount of liquid at the interface.

At the surface of the porous material, we impose a constant atmospheric vapour density and atmospheric pressure

(2.14)\begin{equation} \rho_v=\rho_a, \quad p=p_a, \quad \text{on }y=0. \end{equation}

Dirt cannot diffuse through the impermeable boundary and thus so impose that

(2.15)\begin{equation} \theta_y=0\quad \text{at }y=l. \end{equation}

This depth $l$ is assumed to be much greater than the typical pore-length scale, so that there is separation between the pore- and macro-length scales.

2.2. Non-dimensionalisation

We non-dimensionalise the vapour/gas problem in a similar way to Luckins et al. (Reference Luckins, Breward, Griffiths and Please2023), making the rescalings

(2.16a–f)\begin{align} \boldsymbol{x}=d\hat{\boldsymbol{x}}, \quad h=d\hat{h}, \quad t=\frac{d^2}{\epsilon\delta D_v}\hat{t}, \quad \rho_v=\rho^*\hat{\rho}, \quad \boldsymbol{u}=\frac{D_v\nu\epsilon}{d}\hat{\boldsymbol{u}}, \quad p=p_a+\frac{\mu D_v\nu}{d^2}\hat{p}, \end{align}

where

(2.17a–c)\begin{equation} \delta=\frac{\rho^*}{\rho_L}, \quad \epsilon=\frac{d}{l}, \quad \nu=\frac{\rho^*}{\rho_G}\end{equation}

are dimensionless parameters representing the ratio of vapour density to liquid density, the ratio of pore- to macro-length scales and the ratio of vapour to gas densities, respectively. In particular, we note that we have chosen the time scale associated with the speed of the motion of the evaporating interfaces on the microscale, i.e. the time scale over which sufficient vapour is removed by diffusion to empty the microscale pore space of liquid. Making these rescalings and dropping the hat notation, the dimensionless microscale model is, in $y< h(x,t)$,

(2.18a)\begin{equation} \epsilon\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \quad \epsilon\nabla^2\boldsymbol{u}=\boldsymbol{\nabla} p,\quad \epsilon\left(\delta\rho_t+\nu\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\rho\right)=\nabla^2\rho,\end{equation}

while in $y>h(x,t)$,

(2.18b)\begin{equation} \epsilon\sigma\theta_t=\nabla^2\theta.\end{equation}

At gas–solid interfaces in $y< h(x,t)$ (which are stationary),

(2.18c)\begin{equation} \boldsymbol{u}=\boldsymbol{0}, \quad \boldsymbol{\nabla}\rho\boldsymbol{\cdot}\boldsymbol{n}_s=0,\end{equation}

and at liquid–solid interfaces in $y>h(x,t)$ (which move with velocity $V_n$),

(2.18d)\begin{equation} \boldsymbol{\nabla}\theta\boldsymbol{\cdot}\boldsymbol{n}_s=\epsilon\sigma\left(\theta_*-\theta\right) V_n, \quad V_n=\epsilon\kappa\left(\theta-\beta\right),\end{equation}

while at the evaporating interfaces $y=h(x,t)$,

(2.18e)$$\begin{gather} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{m}-\delta\nu^{{-}1} \frac{h_t}{\sqrt{1+h_x^2}}={-}\frac{h_t}{\sqrt{1+h_x^2}}, \end{gather}$$

(2.18f) $$\begin{gather}\epsilon \rho\Big(\nu\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{m}-\delta \frac{h_t}{\sqrt{1+h_x^2}}\Big)-\boldsymbol{\nabla}\rho\boldsymbol{\cdot} \boldsymbol{m}={-}\epsilon\frac{h_t}{\sqrt{1+h_x^2}}, \end{gather}$$

(2.18g)$$\begin{gather}u+vh_x=0, \end{gather}$$

(2.18h)$$\begin{gather}\epsilon\sigma\theta V_m+ \boldsymbol{\nabla}\theta\boldsymbol{\cdot}\boldsymbol{m}=0, \end{gather}$$

(2.18i)$$\begin{gather}\rho=f(\theta). \end{gather}$$

At the surface of the porous material,

(2.18j)\begin{equation} \rho=\alpha \quad \text{at }y=0,\end{equation}

while at the impermeable surface (or centre of a symmetric porous material),

(2.18k)\begin{equation} \theta_y=0 \quad \text{at }y=\epsilon^{{-}1}.\end{equation}

Here we have introduced the additional dimensionless parameters

(2.19a–d)\begin{equation} \sigma=\delta\frac{D_v}{D_d}, \quad \kappa=\frac{k_+d}{\epsilon^2\delta D_v}, \quad \beta=\frac{k_-\theta_*}{k_+}, \quad \alpha=\frac{\rho_a}{\rho^*} \end{equation}

that appear in the dirt problem, representing the ratio of the suspended dirt diffusion time scale to the time scale of the evaporation-front motion, the ratio of the dirt-deposition rate to the evaporation rate, the ratio of the dirt erosion rate to deposition rate and the ratio of the atmospheric vapour density to the maximum saturation vapour density, respectively.

The micro- to macro-length scale ratio $\epsilon$ (defined in (2.17a–f)) is the small parameter we take advantage of in order to homogenise (2.18). As in Luckins et al. (Reference Luckins, Breward, Griffiths and Please2023), we take $\delta <1$ and $\nu <1$ to be order one parameters relative to $\epsilon$ for the homogenisation analysis. We note that $\delta \approx 10^{-3}$ for water, so will later consider the additional limit of $\delta \rightarrow 0$ (which is equivalent to taking this limit before performing the homogenisation analysis). We require $\alpha <1$ but expect $\alpha \ll 1$ to be reasonable. Although the diffusion of vapour in air is generally much faster than the diffusion of any kind of molecule through a liquid, so that $D_v\gg D_d$, we note that since $\delta \approx 10^{-3}\unicode{x2013}10^{-4}$ is small, $\sigma$ is likely to be order one. For instance, if we take $D_v\approx 2.5\times 10^{-5}\ {\rm m}^2\ {\rm s}^{-1}$ and $D_d\approx 10^{-9}\ {\rm m}^2\ {\rm s}^{-1}$ (Cussler Reference Cussler2009), then with $\delta =10^{-4}$ we find that $\sigma \approx 2.5$. We consider the distinguished limit of $\sigma =O(1)$ in this paper, but we note there is an alternative, slow-dirt-diffusion, distinguished limit with $\sigma =O(\epsilon ^{-1})$. We discuss this alternative case further in Appendix A.2, and briefly in § 2.3 below. In summary, all of $\sigma, \kappa, \beta$ and $\alpha$ are taken to be order one relative to $\epsilon$ to homogenise the model.

2.3. Summary of the homogenised drying model

The homogenisation analysis is described in Appendix A. The result of this analysis is a macroscale model for the vapour density $\rho$, suspended dirt volume fraction $\theta$, deposited dirt-layer thickness, $R$, and position, $Y=H(T)$ of the evaporation front, namely

(2.20a)\begin{equation} \delta\phi\rho_T-(\nu-\delta)\phi|_H H_T\rho_Y= \left(\mathcal{D}\rho_Y\right)_Y\end{equation}

for $Y< H(T)$, and

(2.20b)$$\begin{gather} \sigma\phi\theta_T=(\mathcal{D}\theta_Y)_Y-\sigma\kappa\mathcal{C}(\theta_*-\theta)(\theta-\beta\chi_R), \end{gather}$$

(2.20c)$$\begin{gather}R_T=\kappa(\theta-\beta\chi_R) \end{gather}$$

for $Y>H(T)$. At $Y=H(T)$,

(2.20d)$$\begin{gather} \mathcal{D}\rho_Y=(1-\nu\rho)\phi H_T, \end{gather}$$

(2.20e)$$\begin{gather}\rho=f(\theta), \end{gather}$$

(2.20f)$$\begin{gather}\sigma\phi H_T\theta+\mathcal{D}\theta_Y=0, \end{gather}$$

while

(2.20g)$$\begin{gather} \rho=\alpha \quad\text{on }Y=0, \end{gather}$$

(2.20h)$$\begin{gather}\theta_Y=0 \quad\text{on }Y=1. \end{gather}$$

The porosity, $\phi (R)=1-{\rm \pi} (r_0+R)^2$, surface area, $\mathcal {C}(R)=2{\rm \pi} (r_0+R)$, and effective diffusivity, $\mathcal {D}(R)$ (given by (A4)), all vary with the thickness of the deposited dirt layer.

We assume that, initially, the porous material is entirely saturated with a uniform liquid–dirt mixture, none of which has yet deposited (i.e. the time scale of deposition is assumed longer than the time scale over which the liquid–dirt mixture flooded the material). Thus, at $T=0$,

(2.21a–c)\begin{equation} R=0, \quad H=0, \quad \theta=\theta_{IC}.\end{equation}

Our homogenised model (2.20) is similar in structure to those proposed in Breward et al. (Reference Breward2020), Sanaei et al. (Reference Sanaei2022) and Ji & Sanaei (Reference Ji and Sanaei2023), with the suspended dirt satisfying a reaction–diffusion equation ahead of a moving evaporation front. However, through the systematic homogenisation analysis, we have found the correct form for the diffusion term in (2.20b), which was erroneously given as $(D(\phi \theta )_Y)_Y$ (in our notation) by Ji & Sanaei (Reference Ji and Sanaei2023). Additionally, we have quantified the effective parameters $\mathcal {D}$, $\mathcal {C}$ and $\phi$, which all vary with the deposited dirt-layer thickness $R$. We also impose different boundary conditions to Ji & Sanaei (Reference Ji and Sanaei2023), which will result in different drying behaviours, as discussed in § 3 below.

A key assumption of our homogenisation analysis was that $\sigma =O(1)$, which ensured that $\theta$ (and therefore $R$) is uniform to leading order on the microscale. As discussed further in Appendix A.2, the extremely slow suspended dirt diffusion limit of $\sigma =O(\epsilon ^{-1})$ is not captured by this model: in this case we would expect non-periodic behaviour on the microscale at the evaporating interface, and the homogenisation analysis would break down. We do not consider this situation here.

4.1. Early time analysis

To study the early time behaviour of the system (2.20), we suppose $T=b\tau$ where $b\ll 1$ is the smallest parameter in the system, and $\tau =O(1)$. From (2.20c), on this time scale we see that $R_\tau =O(b\kappa )\ll 1$, and so $R$ is small, hence, all of $\mathcal {D},\phi$ and $\mathcal {C}$ are constant to leading order in $b$.

We first consider the vapour problem in $Y< H$. On the short time scale the interface only moves a short distance, and so we rescale

(4.1a,b)\begin{equation} H=\sqrt{\frac{b\mathcal{D}}{\phi}}\bar{H}, \quad Y=\sqrt{\frac{b\mathcal{D}}{\phi}}\bar{Y} \end{equation}

in order to balance the mass-flux boundary condition (2.20d). The vapour problem is therefore self-similar in that we regain the same system at early time with these rescalings as the full system (2.20a), (2.20d)–(2.20e) and (2.20g), namely

(4.2a)$$\begin{gather} \delta\rho_\tau-(\nu-\delta) \bar{H}_\tau\rho_{\bar{Y}}=\rho_{\bar{Y}\bar{Y}}, \quad\text{for } \bar{Y}\in(0,\bar{H}(\tau)), \end{gather}$$

(4.2b)$$\begin{gather}\rho=\alpha \quad\text{on }\bar{Y}=0, \end{gather}$$

(4.2c)$$\begin{gather}\rho_{\bar{Y}}=(1-\nu\rho) \bar{H}_\tau \quad\text{on }\bar{Y}=\bar{H}(\tau), \end{gather}$$

(4.2d)$$\begin{gather}\rho=f(\theta) \quad\text{on }\bar{Y}=\bar{H}(\tau). \end{gather}$$

We have already noted that $\delta \ll 1$ in general, and we take this limit now to make analytical progress. As in § 3, we find that

(4.3)\begin{equation} \rho=\frac{1}{\nu}\left(1-(1-\nu \alpha)\exp\left(-\nu \bar{H}_\tau \bar{Y}\right)\right),\end{equation}

where $\bar {H}(\tau )$ is the solution of

(4.4)\begin{equation} \bar{H}\bar{H}_\tau=\frac{1}{\nu}\log\Big(\frac{1-\nu \alpha}{1-\nu f(\theta|_{\bar{Y}=\bar{H}(\tau)})}\Big).\end{equation}

The value of $\theta$ at $\bar {Y}=\bar {H}(\tau )$ depends on the solution of the suspended dirt problem in the domain $Y\in (H,1)=(\sqrt {b\mathcal {D}/\phi }\, \bar {H}(\tau ),1)$. On this short time scale, the full dirt problem (2.20b), (2.20f) and (2.20h), with the initial condition (2.21a–c), becomes

(4.5a)$$\begin{gather} \sigma\phi\theta_\tau=b\left(\mathcal{D}\theta_{YY}-\sigma\kappa\mathcal{C}(\theta_*-\theta)(\theta-\beta\chi_R)\right) \quad\text{for }Y\in(\sqrt{b\mathcal{D}/\phi}\, \bar{H}(\tau),1), \end{gather}$$

(4.5b)$$\begin{gather}\sigma \phi\theta \bar{H}_\tau+\sqrt{b}\mathcal{D}\theta_Y=0 \quad\text{on }Y=\sqrt{b\mathcal{D}/\phi}\, \bar{H}(\tau), \end{gather}$$

(4.5c)$$\begin{gather}\theta_Y=0 \quad\text{on }Y=1, \end{gather}$$

(4.5d)$$\begin{gather}\theta=\theta_{IC} \quad\text{at }\tau=0. \end{gather}$$

To leading order in $b$, we see that $\theta _\tau =0$, so that $\theta =\theta _{IC}$ is independent of time over the domain. However, in a boundary layer at $Y=\sqrt {b\mathcal {D}/\phi } \bar {H}$, suspended dirt accumulates due to the motion of the evaporation front. To examine this region, we make the change of variables $Y=\sqrt {b\mathcal {D}/\phi }(\bar {H}(\tau )+z)$, so that, at leading order in $b$, the equations are

(4.6a)$$\begin{gather} \sigma\left(\theta_\tau-\bar{H}_\tau\theta_z\right)=\theta_{zz} \quad\text{for }z>0, \end{gather}$$

(4.6b)$$\begin{gather}\sigma\theta \bar{H}_\tau+\theta_z=0 \quad\text{on }z=0, \end{gather}$$

(4.6c)$$\begin{gather}\theta\rightarrow\theta_{IC} \quad\text{as }z\rightarrow\infty, \end{gather}$$

(4.6d)$$\begin{gather}\theta=\theta_{IC} \quad\text{at }\tau=0. \end{gather}$$

This system (4.6) must be solved with (4.4) to determine $\theta$ and $\bar {H}$.

We look for a similarity solution of the form

(4.7a,b)\begin{equation} \bar{H}=\frac{2\lambda}{\sqrt{\sigma}}{\sqrt{\tau}}, \quad \theta=\varTheta\left(\frac{\sqrt{\sigma}z}{\sqrt{ \tau}}\right),\end{equation}

for some constant $\lambda$ to be determined. In particular, from (4.4) we see that the suspended dirt volume fraction at the evaporating interface, $\theta |_{\bar {H}}=\varTheta (0)$, must be constant in time for such a similarity solution to exist. Substituting into (4.6), we find the solution

(4.8)\begin{equation} \varTheta=\theta_{IC}+\left(\varTheta(0)-\theta_{IC}\right)\frac{\text{erfc}\left(\lambda+ \dfrac{\sqrt{\sigma}z}{2\sqrt{ \tau}}\right)}{\text{erfc}(\lambda)} ,\end{equation}

where $\lambda$ and the constant $\varTheta (0)$ satisfy

(4.9)$$\begin{gather} \lambda^2=\frac{\sigma}{2\nu}\log\left(\frac{1-\nu \alpha}{1-\nu f(\varTheta(0))}\right), \end{gather}$$

(4.10)$$\begin{gather}\varTheta(0)=\theta_{IC}+\sqrt{\rm \pi}\lambda\varTheta(0){\rm e}^{\lambda^2}\,\text{erfc}(\lambda). \end{gather}$$

Solutions of (4.9)–(4.10) may be computed numerically, and are shown for various $\sigma$ and $\theta _{IC}$ in figure 2.

Figure 2.

Early time solution behaviour. Coloured lines show the variation of the solution of (4.9) and (4.10) with the suspended dirt diffusion time scale $\sigma$. Green dashed lines are the small-$\sigma$ approximations (4.11), while black dashed lines are the large-$\sigma$ approximations (4.13). Here we use the form $f(\theta )=1-\theta$ and set $\alpha =0$, so that $\hat {\theta }=1$. We additionally set $\nu =0.5$.

To establish some intuitive understanding of this early time behaviour, we now consider the sublimits $\sigma \ll 1$, $\theta _{IC}\ll 1$ and $\sigma \gg 1$ in turn. We show our early time analytic solutions for each case in figure 3 (alongside numerical solutions for comparison, computed using the method described in § 4.2 below), all with excellent agreement. In each, we see that the vapour density $\rho$ varies from $f(\theta |_H)=1-\theta |_H$ at $Y=H$ to $\alpha =0$ at the surface of the material, according to (4.3). The evaporation front moves with the expected $\sqrt {T}$ behaviour, faster if there is a steeper vapour-density gradient. Suspended dirt accumulates in the liquid ahead of the evaporation front, with a spatial maximum in $\theta$ at $Y=H$. The size of the boundary layer at $H$ over which $\theta$ varies is dependent on $\sigma$, which quantifies suspended dirt diffusion. At early times, we expect little dirt deposition, so that the dirt-layer thickness $R\approx 0$ throughout the porous material.

Figure 3.

Early time solutions (4.3), (4.7a,b) and (4.8) (dotted lines) compared with numerical solutions (solid lines) of the full model (2.20), for $\sigma =0.01$ (a,b) $\sigma =1$, (c,d) and $\sigma =10$ (e,f). The profiles of $\rho$, $\theta$ and $R$ in figures a, c, and e are at times for which $H=0.3$ (towards the end of what we would consider ‘early’ time, especially in the small-$\sigma$ case). Throughout the figure we take $f(\theta )=1-\theta$, $\kappa =1$, $\theta _{IC}=0.1$, $\nu =0.5$, $\delta =10^{-3}$, $r_0=0.2$, $\alpha =0$ and $\beta =0$.

If $\sigma \ll 1$ so that suspended dirt diffusion is fast relative to the motion of the evaporation front, then we see from (4.9) that $\lambda =O(\sqrt {\sigma })$, and so from (4.10) that $\theta |_h\sim \theta _{IC}$. Specifically, we find that

(4.11)\begin{equation} \left.\begin{array}{c} \varTheta(0)=\theta_{IC}+O(\sqrt{\sigma}),\\[3pt] \lambda=\sqrt{\sigma}\left(\sqrt{\dfrac{1}{2\nu}\log\left(\dfrac{1-\nu \alpha}{1-\nu f(\theta_{IC})}\right)}+O (\sqrt{\sigma})\right)\end{array}\!\!\right\} \quad \text{as }\sigma\rightarrow 0.\end{equation}

Thus, reverting to our original variables, the early time evaporating interface is given by

(4.12)\begin{equation} H=\sqrt{T}\sqrt{\frac{\mathcal{D}}{2\nu\phi}\log\left(\frac{1-\nu \alpha}{1-\nu f(\theta_{IC})}\right)}+O(\sqrt{\sigma}) \quad \text{as }\sigma\rightarrow 0. \end{equation}

We also note from the form (4.8) of the solution that the spatial region over which $\theta$ varies is wide, $O(1/\sqrt {\sigma })$. In this small-$\sigma$ limit, the diffusion of dirt is fast relative to the motion of the evaporation front, and so the suspended dirt volume fraction $\theta$ remains close to its initial value $\theta _{IC}$, only deviating by a small, $O(\sqrt {\sigma })$, amount. For conservation of overall dirt, the region over which the accumulating suspended dirt is spread is wide, of $O(1/\sqrt {\sigma })$ relative to the early time boundary layer. This may be seen in figure 3(a,b): since $\sigma \ll 1$, the $\theta$ profile is approximately uniform in $Y$, and so close to its initial value of $\theta _{IC}=0.1$ at early times. The suspended dirt is accumulating due to the evaporation, but spread almost evenly through the domain.

Next, we suppose that $\sigma =O(1)$ but the initial suspended dirt volume fraction $\theta _{IC}\ll 1$ is small. In this case, for a balance in both of (4.9) and (4.10), we must have $\varTheta (0)=O(\theta _{IC})$ and $\lambda =O(1)$, as we might expect. The solution shown in figure 3(c,d) is for this case, with $\theta _{IC}=0.1$: we indeed observe that $\theta |_H=O(0.1)$ (and this effect becomes increasingly clear for smaller $\theta _{IC}$).

Finally, if $\sigma \gg 1$, so that the diffusion of suspended dirt is slow relative to the motion of the evaporation front, then from (4.9) we see that we must have $f(\varTheta (0))=\alpha$ to leading order in $\sigma ^{-1}\ll 1$. At this value,

(4.13a)\begin{equation} \varTheta(0)=\hat{\theta}, \end{equation}

there is no evaporation at leading order, as the vapour density at the atmospheric value is in equilibrium with the liquid–dirt interface, and there is no transport of vapour out of the porous material. Indeed, we see from (4.10) that when $\theta =\hat {\theta }$, $\lambda$ is the solution of

(4.13b)\begin{equation} \lambda {\rm e}^{\lambda^2}\,\text{erfc}(\lambda)=\frac{\hat{\theta}-\theta_{IC}}{\hat{\theta}\sqrt{\rm \pi}}, \end{equation}

which is independent of $\sigma$ and of $O(1)$, so that the position of the evaporating interface, given by

(4.14)\begin{equation} H=\frac{2\lambda}{\sqrt{\sigma}}\sqrt{\frac{\mathcal{D}}{\phi}}\sqrt{T}, \end{equation}

is of order $\sigma ^{-1/2}\ll 1$ away from its initial position. Thus, when the suspended dirt diffusion is slow, the diffusion of dirt away from $H$ limits the speed of the evaporation front, so that there is a slower, $O(\sigma )$, drying time scale. We also note from (4.8) that the region over which $\theta$ varies is narrow, with width $O(1/\sqrt {\sigma })$ relative to the early time boundary layer. The boundary layer is narrow for large $\sigma$, so that the early time solution actually remains valid for the majority of the drying process. Indeed, in figure 3(e,f) we see excellent agreement between the early time analytic solution and the numerical solution for $\rho, \theta$ and $H$ up to the time when the evaporating front is halfway through the domain. This is because the boundary layer at $H$ over which $\theta$ varies is narrow, and so the effect of the boundary at $Y=1$ is not felt until $H$ is close to 1. However, we notice in figure 3(e) that the early time approximation $R=0$ ceases to be accurate at these late times. The early time solution would remain valid so long as $R$ remains relatively small (e.g. if $\kappa$ and $\theta _{IC}$ are fairly small). We note that, since the boundary layer width scales with $\sqrt {T}$, it quickly becomes numerically impractical to resolve the solution at small times for large $\sigma$. The early time asymptotic solution is therefore very valuable in initialising the simulations accurately for large $\sigma$.

Finally, we note that our early time analysis in this section is equivalent to studying the original problem on a semi-infinite domain $Y\in (0,\infty )$, in the combined limit $\kappa,\delta \rightarrow 0$.

4.2. Numerical method

We solve the model (2.20) numerically using the method of lines. Specifically, we first transform the model onto two separate fixed domains, setting $\eta =Y/H(T)$ for the gas–vapour problem, which then holds in $\eta \in (0,1)$, and setting $\xi =(Y-H(T))/(1-H(T))$ for the liquid–dirt problem, which then also holds in $\xi \in (0,1)$. We discretise spatially on these transformed domains, with a uniform mesh, using central differences for diffusive terms and first-order upwinding for advective terms, so that the scheme is overall first order. (The advection for the vapour problem (2.20a), including the artificial advection terms due to the change of variables, is negative; the purely artificial advection in the liquid–dirt problem (2.20b)–(2.20c) is also negative. Upwinding these terms therefore requires forward differences in both cases.) We then use the inbuilt ODE solver ode15s in Matlab for the time stepping. We note that the model is stiff in certain parameter regimes of interest ($\delta \ll 1$ and/or $\sigma \ll 1$), and that ode15s is specifically designed for stiff systems. Ode15s is a multistep solver, using numerical differentiation formulas of order 1–5 (Shampine & Reichelt Reference Shampine and Reichelt1997). We make use of our early time asymptotic solution of § 4.1 to initialise our numerical simulations. In particular, the spatial mesh must be sufficiently fine to resolve the boundary layer in the suspended dirt problem at $Y=H$ at early times. Our analysis in § 4.1 suggests we require the number of spatial mesh points $N$ to scale like

(4.15)\begin{equation} N=O\left(\sqrt{\frac{\mathcal{D}(0)\sigma}{\phi(0)T}}\right) \quad \text{for }T\ll1. \end{equation}

More efficient solvers might take further advantage of the asymptotic structure of the system and distribute mesh points unevenly through the domain in order to ensure good resolution of the boundary layer while maintaining computational efficiency. However, by making use of our early time asymptotic solution we do not require simulations at particularly small $T$, and our uniform-mesh formulation suffices.

5.1. Infinitely slow evaporation when $\kappa =0$

Taking $\kappa =0$, we see from (2.20c) that we have $R=0$ everywhere, and thus, $\mathcal {D}=\mathcal {D}_0$, $\mathcal {C}=\mathcal {C}_0$ and $\phi =\phi _0$ are all constant, equal to their values at $R=0$. Taking the $\delta \ll 1$ limit as in § 3, the motion of the evaporation front is therefore governed by (3.4), i.e.

(5.1)\begin{equation} HH_T={-}\frac{\mathcal{D}_0}{\nu\phi_0}\log\left(\frac{1-\nu f(\theta|_H)}{1-\nu \alpha}\right),\end{equation}

while $\theta$ satisfies

(5.2a)$$\begin{gather} \frac{\sigma\phi_0}{\mathcal{D}_0}\theta_T=\theta_{YY} \quad\text{for }Y\in(H(T),1), \end{gather}$$

(5.2b)$$\begin{gather}\frac{\sigma\phi_0}{\mathcal{D}_0}\theta H_T+\theta_Y=0 \quad\text{on }Y=H(T), \end{gather}$$

(5.2c)$$\begin{gather}\theta_Y=0 \quad\text{on }Y=1, \end{gather}$$

(5.2d)$$\begin{gather}\theta=\theta_{IC} \quad\text{at }T=0. \end{gather}$$

To investigate how the accumulation of suspended dirt affects the evaporation rate, we consider the additional limit of $\sigma \ll 1$ so that the diffusion of suspended dirt is fast. In this case, we see from (5.2) that $\theta (T)$ is uniform, and so, for overall conservation of dirt, we must have

(5.3)\begin{equation} \theta(T)=\frac{\theta_{IC}}{1-H(T)}.\end{equation}

Substituting (5.3) into (5.1) we obtain the single equation for $H(T)$,

(5.4)\begin{equation} HH_T={-}\frac{\mathcal{D}_0}{\nu\phi_0}\left[\log\left(1-\nu f\left(\frac{\theta_{IC}}{1-H}\right)\right)-\log(1-\nu \alpha)\right].\end{equation}

As discussed previously, the evaporation shuts down when $\theta =\hat {\theta }$ so that $f(\hat {\theta })=\alpha$, since then $H_T=0$. At this point we see from (5.3) that $H=1-\theta _{IC}/\hat {\theta }$.

Numerical solutions of the model (2.20) (with $\kappa =0$, $\delta =10^{-3}$) are shown in figures 4(a) and 4(b), and compared with the solution of (5.4) for the limit of $\sigma \rightarrow 0$, with good agreement for $\sigma =0.1$ and smaller. We take the functional form $f(\theta )=1-\theta$ for these simulations, and fix $\alpha =0$ (so that $\hat {\theta }=1$). In figure 4(c) we show solutions of (5.4) for various $\theta _{IC}$. We see that, for larger $\theta _{IC}$, the evaporation is slower, with the evaporating interface $H$ moving more slowly into the domain. In particular, we note that when $\theta _{IC}=0$, the evaporation is completed (with $H=1$) in finite time

(5.5)\begin{equation} T={-}\frac{\nu\phi_0}{2\mathcal{D}_0\log(1-\nu)}\approx 0.36, \end{equation}

(from (5.4) with $\theta _{IC}=0$), whereas for $\theta _{IC}>0$, we see that $H$ appears to take infinite time to reach $1-\theta _{IC}/\hat {\theta }$.

Figure 4.

The effect of suspended dirt accumulation on the evaporation, for $\kappa =0$ and $\sigma \ll 1$. Numerical solutions of (2.20) shown with $\kappa =0$ (no dirt deposition), alongside the solution of (5.4) in the limit $\sigma \rightarrow 0$, (taking $f(\theta )=1-\theta$ and $\alpha =0$). Throughout the figure we take $\nu =0.5$ and $r_0=0.2$. (a) Motion of the evaporation front $H(T)$ for various $\sigma$, with $\theta _{IC}=0.3$. (b) Vapour density and suspended dirt volume fraction profiles at the same time $T=0.2$ for various $\sigma$, with $\theta _{IC}=0.3$. (c) Effect of $\theta _{IC}$ on the evaporating interface motion. Dashed lines are $1-\theta _{IC}$.

We investigate this late time behaviour (within the $\sigma \ll 1$ limit) by considering the expansion

(5.6)\begin{equation} H=1-\frac{\theta_{IC}}{\hat{\theta}}(1+c\mathcal{H}) \quad \text{so that}\ \theta=\frac{\hat{\theta}}{1+c\mathcal{H}},\end{equation}

where $c\ll 1$ is small and $\mathcal {H}=O(1)$ is positive. Assuming that $f$ is continuous at $\theta =\hat {\theta }$, on substitution of (5.6) into (5.4) we find that (retaining only leading-order terms on either side)

(5.7) \begin{align} c\left(1-\frac{\theta_{IC}}{\hat{\theta}}\right)\frac{\theta_{IC}}{\hat{\theta}}\mathcal{H}_T&={-} \frac{\mathcal{D}_0}{\phi_0(1-\nu \alpha)}\biggl(f\biggl(\frac{\hat{\theta}}{1+c\mathcal{H}}\biggr)-\alpha\biggr)\nonumber\\ &={-}\frac{\mathcal{D}_0}{\phi_0(1-\nu \alpha)}( f(\hat{\theta}(1-c\mathcal{H}+O(c^2)))-\alpha). \end{align}

So long as the gradient of the function $f$ is bounded at $\hat {\theta }$, we may Taylor expand the right-hand side of (5.7) and, thus, find that, to leading order in $c$,

(5.8)\begin{equation} \mathcal{H}_T=\frac{\mathcal{D}_0\hat{\theta}^3}{\phi_0(1-\nu \alpha)\theta_{IC}(\hat{\theta}-\theta_{IC})}f'(\hat{\theta})\mathcal{H}, \end{equation}

since $f(\hat {\theta })=\alpha$ by definition. Thus, $\mathcal {H}$ decays exponentially to zero as $T\rightarrow \infty$ (since $f'(\hat {\theta })<0$) and the evaporation takes infinite time. (Similarly, if $f'(\hat {\theta })=0$ but higher derivatives are non-zero then $\mathcal {H}$ has polynomial, and still infinite-time, decay.)

We might ask if there are sensible choices for $f(\theta )$ that result in finite-time completion of the evaporation. In order for the evaporation to complete in finite time, we see that the gradient of $f$ must be unbounded at $\hat {\theta }$. For instance, if $f(\hat {\theta }(1-c\mathcal {H}))=\alpha +A(c\mathcal {H})^{1/n}$, for $n> 1$ and some constant $A$, then by substituting this expansion for $f$ into (5.7) and integrating the resulting ODE for $\mathcal {H}$ in time $T$, we find that

(5.9) \begin{equation} \mathcal{H}\sim\Big(\text{const.}-\frac{D_0\hat{\theta}^2A(n-1)}{\phi_0\theta_{IC}(\hat{\theta}-\theta_{IC})(1-\nu \alpha) n}c^{-(n-1)/n}T\Big)^{n/(n-1)}, \end{equation}

and so $\mathcal {H}$ reaches zero in finite time. However, this finite-time drying requires that the gradient of $f$ is unbounded at $\hat {\theta }$, which is physically unreasonable, not least because $\hat {\theta }$ (defined as the value of $\theta$ for which $f=\alpha$) depends on the atmospheric vapour density $\rho _a$ via the value of $\alpha$. For physically reasonable functional forms $f(\theta )$, we therefore expect a bounded derivative at $\hat {\theta }$, and thus, that the evaporation becomes unboundedly slow as $\theta \rightarrow \hat {\theta }$.

This analysis is for the case $\sigma \ll 1$. For larger $\sigma$, we see from figure 4(a) that the evaporation rate is slower. As we see in figure 4(b), this is because suspended dirt accumulates near the evaporating interface rather than quickly diffusing through the domain, and with higher values of $\theta |_H$, we see from (5.1) that the evaporation rate is reduced. Numerically, we see that, for non-negligible $\sigma$, we still have $H\rightarrow 1-\theta _{IC}/\hat {\theta }$ in infinite time. Indeed, although for larger $\sigma$ the evaporation is additionally slowed by the diffusion of the suspended dirt away from the evaporating interface, at late times when the evaporation becomes infinitely slow, the diffusion of dirt does not limit the drying process.

In summary, for $\kappa =0$, the drying takes infinite time to complete and, as drying occurs, the suspended dirt concentrates in a layer at $y=1$. The drying never fully stops (although it becomes infinitely slow). By contrast, we see in § 5.2 that if $\kappa \ll 1$ is non-zero then the dirt begins to deposit at late time, and this causes the drying to completely stop in finite time (which we refer to as clogging).

5.2. Dry-clogging behaviour for small but non-zero $\kappa$

The analysis in § 5.1 assumed that $\kappa =0$ so that there was no deposition of dirt at all, and we saw that, in this case, there is an infinitely long drying time as $\theta \rightarrow \hat {\theta }$. However, in reality we might instead have $\kappa \ll 1$ small but non-zero. In this case, we would expect the same behaviour as in § 5 initially (over an $O(1)$ time), but then during the slow evolution towards $\theta \rightarrow \hat {\theta }$, the deposition will begin to become important. Dirt is slowly deposited, and the deposited dirt layer, thickness $R$, grows. From the microscale geometry, we see that, when $R=R_{clog}:=0.5-r_0$, the dirt layers on neighbouring solid circles meet, and the pore-scale liquid region ceases to be connected. This means that the effective diffusivity $\mathcal {D}(R_{clog})=0$. In particular, if $R=R_{clog}$ then vapour cannot be transported through the porous material, and thus, evaporation ceases. We define ‘clogging’ to be this situation when $R=R_{clog}$ at $Y=H(T)$ at finite time $T$, and thus, the evaporation is stopped.

To investigate this, we look at the behaviour of the system on the long time $T=\tilde {T}/\kappa$ over which $R$ varies. With this change of variables in (5.1), we see that $H$ satisfies

(5.10)\begin{equation} \kappa H_{\tilde{T}}\int_{0}^H\frac{1}{\mathcal{D}}\,\mathrm{d}Y={-}\frac{1}{\nu\phi}\log\left(\frac{1-\nu f(\theta|_H)}{1-\nu \alpha}\right),\end{equation}

while, for $Y>H$,

(5.11)$$\begin{gather} R_{\tilde{T}}=\theta-\beta\chi_R, \end{gather}$$

(5.12)$$\begin{gather}\kappa\sigma\phi\theta_{\tilde{T}}=\left(\mathcal{D}\theta_Y\right)_Y-\sigma\kappa\mathcal{C}(\theta_*-\theta)(\theta-\beta\chi_R), \end{gather}$$

along with the boundary conditions (2.20f) and (2.20h),

(5.13)$$\begin{gather} \kappa\sigma\phi H_{\tilde{T}}\theta+\mathcal{D}\theta_Y=0 \quad\text{on }Y=H(\tilde{T}), \end{gather}$$

(5.14)$$\begin{gather}\theta_Y=0\quad\text{on }Y=1. \end{gather}$$

At leading order in $\kappa$, we see from (5.12)–(5.14) that $\theta =\hat {\theta }$ is uniform, and thus, in $Y>H$, from (5.11) we find that

(5.15)\begin{equation} R=(\hat{\theta}-\beta){\tilde{T}} \end{equation}

is spatially uniform. Clearly $R=R_{clog}$ after time ${\tilde {T}}=R_{clog}/(\hat {\theta }-\beta )$, or in the original time variable, at

(5.16)\begin{equation} T_{end}=\frac{R_{clog}}{\kappa(\hat{\theta}-\beta)}+O(1).\end{equation}

We note that, during this late, $O(1/\kappa )$ time, the position of the evaporating interface, $H$, and the $O(\kappa )$ deviation of $\theta$ from $\hat {\theta }$ may be found by going to next order in $\kappa$ in (5.10) and (5.12), and matching to the early time behaviour in ${\tilde {T}}$ (or $O(1)$ time behaviour in $T$). Thus, $R$ reaches the clogging point $R_{clog}$ in finite time, and thus, the system clogs. We term this type of clogging ‘dry’ clogging because, to leading order, $\theta =\hat {\theta }$ everywhere ahead of the clogging front. Since we expect $\hat {\theta }\approx 1$ (indeed we take $\hat {\theta }=1$ in our numerical simulations), there is therefore a negligible amount of liquid left trapped in the porous material by the clogging. This dry-clogging behaviour is illustrated in figure 5(a).

Figure 5.

The two clogging behaviours: dry clogging (a) for which $\theta =\hat {\theta }$ for $Y>H$ when the system clogs, so that (almost) pure dirt remains with negligible liquid trapped; and wet clogging (b) for which a mixture of liquid and dirt is trapped in $Y>H$ by the clogging.

Results of a numerical simulation of (2.20) for $\kappa =0.04$ are shown in figure 6. As expected, we see in figure 6(b) that $H$ varies from zero to around $1-\theta _{IC}/\hat {\theta }$ over an $O(1)$ time, while $R$ at the interface $Y=H$ (where $R$ is maximised in space at that time $T$) remains closer to zero during the time for which $H$ varies. Subsequently, there is a longer $O(1/\kappa )$ time over which $H$ remains nearly stationary, since $\theta \approx \hat {\theta }$ everywhere in $Y>H$, while $R$ (at $Y=H$) increases linearly to $R_{clog}=0.3$. The prediction (5.16) gives $T_{end}=7.5$ for the parameter values used in figure 6(b), which is seen to be fairly accurate, although a slight underestimate, as this does not take into account the early time (in ${\tilde {T}}$, or $O(1)$ in $T$) stage.

Figure 6.

Numerical solution of (2.20) for small $\kappa$ showing dry-clogging behaviour, using parameter values $\kappa =0.04$, $\sigma =1$, $\theta _{IC}=0.3$. (a) Profiles of $\rho$, $\theta$ and $R$ very close to the clogging time $T_{end}$. (b) The motion $H(T)$ of the evaporation front.

The dry clogging that we have described in this section always occurs for $\kappa \ll 1$, but it may also occur for $\kappa =O(1)$, when the deposition rate is of the same order as the evaporation rate. Indeed, the model (2.20) must always dry clog for $\kappa >0$, even if $\beta =0$. This is because, if $\beta =0$, $\theta$ can never reach zero even for large $\kappa$, and can only decay exponentially towards it. However, if $\kappa \gg 1$ is sufficiently large that $\theta$ is very close to zero by the end stages of the drying, the dry clogging occurs at a negligible distance from the end of the domain, $H=1$, and is not physically meaningful. We discuss this more in § 7 below. Furthermore, if $\beta >0$ then we expect $\theta \geq \beta$ for all time (so long as this is true initially). In this case we would certainly anticipate much more prominent dry-clogging behaviour, for a wider range of parameters, although to investigate this thoroughly is beyond the scope of the present study.

6.1. Large-$\kappa$ behaviour

To understand the deposition (and potential clogging) behaviour when $\kappa \gg 1$, we change to the fast time scale over which $R$ varies, by setting

(6.1)\begin{equation} T=\frac{1}{\kappa}\bar{t}, \end{equation}

where $\bar {t}=O(1)$. At such early times the evaporating interface is close to the surface of the porous material, and we rescale

(6.2a,b)\begin{equation} H=\frac{1}{\sqrt{\sigma\kappa}}\bar{h}, \quad Y=\frac{1}{\sqrt{\sigma\kappa}}\bar{y}, \end{equation}

for $Y< H(T)$ (assuming that $\sigma \kappa \gg 1$), so that (3.3) becomes

(6.3)\begin{equation} \frac{1}{\sigma}\bar{h}_{\bar{t}}\int_0^{\bar{h}}\frac{1}{\mathcal{D}(R(\bar{y}))}\,\mathrm{d}\bar{y}={-}\frac{1}{\nu\phi|_{\bar{h}}} \log\left(\frac{1-\nu f(\theta|_h)}{1-\nu \alpha}\right).\end{equation}

For $Y>H=(1/\sqrt {\sigma \kappa })\bar {h}$, we see that (2.20b)–(2.20c) become

(6.4)$$\begin{gather} \phi\theta_{\bar{t}}=\frac{1}{\sigma\kappa}\left(\mathcal{D}\theta_Y\right)_Y-\mathcal{C}(\theta_*-\theta) (\theta-\beta\chi_R), \end{gather}$$

(6.5)$$\begin{gather}R_{\bar{t}}=\theta-\beta\chi_R. \end{gather}$$

To leading order in $(\kappa \sigma )^{-1}\ll 1$, we have a plane-autonomous deposition system:

(6.6a)$$\begin{gather} \phi(R)\theta_{\bar{t}}={-}\mathcal{C}(R)(\theta_*-\theta)(\theta-\beta\chi_R), \end{gather}$$

(6.6b)$$\begin{gather}R_{\bar{t}}=\theta-\beta\chi_R. \end{gather}$$

We refer to the system (6.6) as the outer problem, which holds away from the evaporation front in the majority of the domain. (At the evaporation front there must be a boundary layer, which we discuss later.) Specifically, the initial conditions $\theta =\theta _{IC}>\beta$ and $R=0$ at $\bar {t}=0$, imply that both $\theta$ and $R$ are independent of $Y$ for all $\bar {t}$, and, recalling that $\phi (R)=1-{\rm \pi} (r_0+R)^2$ and $\mathcal {C}(R)=-\phi '(R)=2{\rm \pi} (r_0+R)$, a first integral from (6.6) is

(6.7)\begin{equation} \phi(R)(\theta_*-\theta)=\phi(0)(\theta_*-\theta_{IC}),\end{equation}

which is independent of time $\bar {t}$. Equation (6.7) may be interpreted as an expression of overall conservation of dirt: since there is no transport of dirt on this time scale, the total suspended dirt in the liquid and deposited dirt in the layer, must remain constant in time. We could use (6.7) in (6.6) to find implicit expressions for the spatially uniform $R$ and $\theta$, although this is not particularly illustrative. Instead, we use (6.7) to plot the phase plane in the outer region in figure 7. The system begins at $R=0$, $\theta =\theta _{IC}$. If $\theta _{IC}>\beta$, we see from (6.6b) that $R$ increases in time and, from (6.7), that $\theta$ decreases towards $\theta =\beta$. If instead $\theta _{IC}\leq \beta$, then $R$ remains zero and $\theta$ remains at its initial value (as any deposited dirt immediately re-suspends, from (6.6b)). We note that the system clogs if $R$ reaches $R_{clog}$ before $\theta$ reaches $\beta$. We see that this occurs for sufficiently large initial suspended dirt concentrations, namely (from (6.7)) if

(6.8)\begin{equation} \theta_{IC}>\theta_{IC}^{crit}:=\theta_*-(\theta_*-\beta)\Big(\frac{1-{\rm \pi}/4}{1-{\rm \pi} r_0^2}\Big) .\end{equation}

Although this analysis shows that the system certainly clogs for $\theta _{IC}>\theta _{IC}^{crit}$ (in this limit $\kappa \sigma \gg 1$), we expect that the system will in fact clog for lower values of $\theta _{IC}$ too, since we expect there will be higher $\theta$ and, therefore, faster deposition near to the evaporation front at $\bar {h}$.

Figure 7.

Phase plane dynamics in the outer region when $\kappa \gg 1$. Black curves show trajectories, with the system moving down these curves from $(0,\theta _{IC})$ to the attractor at $\theta =\beta$ (direction shown by arrows). For $\theta _{IC}>\theta _{IC}^{crit}$, the trajectory hits $R=R_{clog}$ before reaching $\theta =\beta$. Here we take $\theta _*=1$, $\beta =0.15$, $r_0=0.2$.

Indeed, in a boundary layer of width $O(1/\sqrt {\sigma \kappa })$, dirt accumulates due to the motion of the evaporation front. By making the change of variables

(6.9)\begin{equation} Y=\frac{1}{\sqrt{\sigma\kappa}}\left(\bar{h}+z\right), \end{equation}

we find that, in the boundary layer $z\in (0,\infty )$,

(6.10a)$$\begin{gather} \phi\left(\theta_{\bar{t}}-\bar{h}_{\bar{t}}\theta_z\right)= \left(\mathcal{D}\theta_z\right)_z-\mathcal{C}(\theta_*-\theta)(\theta-\beta\chi_R), \end{gather}$$

(6.10b)$$\begin{gather}R_{\bar{t}}-\bar{h}_{\bar{t}}R_z=\theta-\beta\chi_R, \end{gather}$$

while at the evaporation interface $z=0$,

(6.10c)\begin{equation} \phi\theta\bar{h}_{\bar{t}}+\mathcal{D}\theta_{z}=0, \end{equation}

and, as $z\rightarrow \infty$,

(6.10d)\begin{equation} R\rightarrow R^{out}(\bar{t}), \quad \theta\rightarrow\theta^{out}(\bar{t}), \end{equation}

where $R^{out}$ and $\theta ^{out}$ are the values in the outer region, as described above. This boundary layer system (6.10), coupled with (6.3), describes the motion of the evaporation front, accumulation of suspended dirt ahead of it and the deposition of dirt. We note that dirt accumulation, transport and deposition all balance in the boundary layer.

We show a solution with $\kappa =100$ in figure 8. Figures 8(a)–8(d) show the spatial profiles for $\rho$, $\theta$ and $R$ at four successive times, while the motion of the evaporation front $H(T)$ is shown in figure 8(e), with the time points of the plots (a–d) marked as red circles. At early times, in figures 8(a) and 8(b), we clearly see the boundary layer structure in the $\theta$ and $R$ plots in $Y>H$, with both $R$ and $\theta$ uniform in the rest of the domain. As time progresses, we see that $R$ increases while $\theta$ decreases, as suspended dirt becomes deposited onto the solid structure. By the time shown in figure 8(c), we see that $\theta$ is close to zero everywhere: the deposition has nearly finished. Since $\kappa \gg 1$, this occurs when the evaporation front is still only a short $O(\sqrt {\kappa })$ distance into the domain. After this time, $\theta ({\approx }0)$ and $R$ are both constant in $Y>H$, and the evaporation front travels to the bottom of the domain, resulting in a fully dry porous material with non-uniformly deposited dirt. Indeed, we note that the combination of dirt accumulation, diffusion and deposition in the boundary layer results in an internal spatial peak in the final thickness, $R$, of the deposited dirt layer (for $Y< H$) near the top of the domain. Since $R$ is higher here, the effective diffusivity of the vapour is correspondingly reduced. This can be seen in the non-monotone gradient of $\rho$ in figure 8(d), where the $\rho$ profile is steepest at the peak value of $R$. The reduced diffusivity here limits the drying rate for the remainder of the process.

Figure 8.

Numerical solution of (2.20) for large deposition rate $\kappa =100$. We also take $f(\theta )=1-\theta$, $\alpha =0$, $\beta =0$, $\nu =0.5$, $\sigma =1$, $\theta _{IC}=0.4$, $\delta =10^{-3}$ and $r_0=0.2$. The green dashed lines in figures (a–d) show the clogging point, $R_{clog}=0.5-r_0$, which is not attained in this simulation. Results are shown for (a) $T=0.0025$, (b) $T=0.005$, (c) $T=0.01$, (d) $T=0.1$, (e) Motion of the evaporation front.

The fact that we obtain this internal peak in $R$ at early time due to the boundary layer accumulation, diffusion and deposition of dirt means that our estimate for the clogging criterion (6.8), which assumes that dirt deposits in a spatially uniform way, must be an upper bound on the true critical $\theta _{IC}$: we expect to have clogging at $\theta _{IC}$ lower than the critical value given by (6.8). Indeed, in figure 9 we show a simulation for the same parameter values as in figure 8, except that we take a larger initial suspended dirt volume fraction $\theta _{IC}=0.6$, which is still lower than the estimate of $\theta _{IC}^{crit}=0.755$ given by (6.8), for the chosen parameter values. Nevertheless, we see that the system indeed clogs early in the domain. Both liquid and suspended dirt are trapped ahead of the clogged point, since $\theta <1$ in $Y>H$. The clogging happens at an $O(1/\sqrt {\kappa })$ distance into the domain, at $H\approx 0.07$, and after an $O(1/\kappa )$ time, as predicted by our analysis above.

Figure 9.

Wet clogging: numerical solution of (2.20) for large deposition rate $\kappa =100$ and $\theta _{IC}=0.6$. We also take $f(\theta )=1-\theta$, $\alpha =0$, $\beta =0$, $\sigma =1$, $\nu =0.5$, $\delta =10^{-3}$ and $r_0=0.2$. The green dashed line shows the clogging point, $R_{clog}=0.5-r_0$, which is attained at $T\approx 5.1\times 10^{-3}$ in this simulation. The profiles in (a) are at the time when the system clogs with $R=R_{clog}$ at $Y=H$.

We note that the motion of the evaporation front shown in figure 9(b) no longer follows a $\sqrt {T}$ behaviour. Instead, we see that the speed $H_T$ of the evaporation front is not smoothly decreasing. Evaporation is fast to start with as $R$ is small and the vapour has only a short way to travel to the surface of the porous material. Then the evaporation front begins to slow down, as both $\theta |_H$ increases due to accumulation and the effective diffusivity $\mathcal {D}$ decreases, as $R$ begins to increase. As the clogging point is approached, $H_T$ increases again, because the porosity $\phi$ is decreasing, so that there is less liquid in the pore space to be evaporated since so much of the pore space is occupied by the deposited and suspended dirt. Similar time-varying behaviour of $H_T$ is visible at early times in the case shown in figure 8(e) (inset), when the system did not clog.

Evidently, wet clogging can occur at initial suspended dirt volume fractions $\theta _{IC}$ significantly below the estimate for the critical value given by (6.8). To better understand the clogging criterion, we must better understand the behaviour in the boundary layer at $Y=H$. However, the nonlinearity of the boundary layer equations (6.10) makes them intractable to further analytical progress. Instead, we study a paradigm problem in the next section, which is a simplified version of the large-$\kappa$ problem, but which still captures the essence of the wet-clogging behaviour.

6.2. The clogging criterion for a paradigm problem

In § 6.1 we saw numerically that, for large $\kappa$, the deposited dirt-layer thickness $R$ is non-uniform near $Y=0$, increasing from zero to $R_{max}:=\max _T(R|_H)$ that is larger than the final value of $R$ attained in the outer region. This means that wet clogging occurs for lower values of $\theta _{IC}$ than predicted by the outer region estimate (6.8). Since the spatially uniform deposition dynamics in the outer region do not capture this non-uniformity in $R$ near to $Y=0$, we consider the behaviour in the boundary layer at $Y=0$, where there is a balance between all of the accumulation, diffusion and deposition of dirt processes, in order to derive an improved criterion for wet clogging.

However, the boundary layer equations (6.10) are intractable analytically, due to the coupling between variables and the nonlinear terms. In order to build intuition for what determines the size of the peak $R_{max}$, in this section we investigate a paradigm problem, with a different functional form for $f(\theta )$, and unphysical simplifications of $\mathcal {D}$ and the bulk deposition term. With these (non-systematic) simplifications, we solve the boundary layer equations (6.10) explicitly, and hence, determine $R_{max}$ analytically. In this way, we determine a criterion for wet clogging (given by $R_{max}\geq R_{clog}$), and build intuition for the more general case.

We assume, for simplicity, that $\beta =0$ and $\alpha =0$. Additionally, we suppose that $r_0$ is close to $1/2$, so that $R\ll r_0$ (since the circular inclusions are close together, and only thin deposited dirt layers are possible before clogging occurs). Then $\phi \approx \phi _0$ and $\mathcal {C}\approx \mathcal {C}_0$ are approximately constant, even for $R$ approaching $R_{clog}$. We make the additional approximations

(6.11a,b)\begin{equation} \mathcal{D}=\begin{cases}\mathcal{D}_0 & \text{if }R< R_{clog},\\ 0 & \text{if }R=R_{clog},\end{cases} \quad f(\theta)=\begin{cases} 1 & \text{if }\theta<1,\\ 0 & \text{if }\theta=1.\end{cases}\end{equation}

Finally, we alter the deposition term in (6.10) by replacing the factor $\theta _*-\theta$ with the constant value $\theta _*$. With these choices of functional forms (which we emphasise are not a true limit of the full problem), the paradigm boundary layer problem is, while $\theta <1$,

(6.12a)\begin{equation} \frac{1}{\sigma \mathcal{D}_0}\bar{h}\bar{h}_{\bar{t}}={-}\frac{1}{\nu\phi_0} \log\left(1-\nu \right),\end{equation}

along with, in $z\in (0,\infty )$ and while $R< R_{clog}$,

(6.12b)$$\begin{gather} \phi_0\left(\theta_{\bar{t}}-\bar{h}_{\bar{t}}\theta_z\right)= \mathcal{D}_0\theta_{zz}-\mathcal{C}_0\theta_*\theta, \end{gather}$$

(6.12c)$$\begin{gather}R_{\bar{t}}-\bar{h}_{\bar{t}}R_z=\theta, \end{gather}$$

while at the evaporation interface $z=0$,

(6.12d)\begin{equation} \phi_0\theta\bar{h}_{\bar{t}}+\mathcal{D}_0\theta_{z}=0, \end{equation}

and, as $z\rightarrow \infty$,

(6.12e)\begin{equation} R\rightarrow R^{out}(\bar{t}), \quad \theta\rightarrow\theta^{out}(\bar{t}),\end{equation}

where the outer solution $R^{out}$, $\theta ^{out}$ satisfy

(6.12f)$$\begin{gather} \phi_0\theta^{out}_{\bar{t}}={-}\mathcal{C}_0\theta_*\theta^{out}, \end{gather}$$

(6.12g)$$\begin{gather}R^{out}_{\bar{t}}=\theta^{out}. \end{gather}$$

Initially, at $\bar {t}=0$,

(6.12h)\begin{equation} \theta=\theta_{IC}, \quad R=0, \quad \bar{h}=0. \end{equation}

Solving (6.12), we see from (6.12a) that the evaporation front is simply

(6.13)\begin{equation} \bar{h}=B\sqrt{D\bar{t}}, \end{equation}

where, for simplicity of notation, we define

(6.14a,b)\begin{equation} B:=\sqrt{-\frac{2\sigma}{\nu}\log(1-\nu)}, \quad D:=\frac{\mathcal{D}_0}{\phi_0}.\end{equation}

Furthermore, the outer region phase plane equations (6.12f)–(6.12g) decouple, with the explicit solution

(6.15a,b)\begin{equation} \theta^{out}=\theta_{IC}\exp\left({-}C\bar{t}\right), \quad R^{out}=\frac{\theta_{IC}}{C}\left(1-\exp\left({-}C\bar{t}\right)\right), \end{equation}

where, again for simplicity of notation, we define

(6.16)\begin{equation} C:=\frac{\mathcal{C}_0\theta_*}{\phi_0}.\end{equation}

The boundary layer equations (6.12b)–(6.12e) are therefore

(6.17a)\begin{equation} \theta_{\bar{t}}-\frac{B\sqrt{D}}{\sqrt{\bar{t}}}\theta_z=D\theta_{zz}-C\theta, \quad R_{\bar{t}}-\frac{B\sqrt{D}}{\sqrt{{\bar{t}}}} R_z=\theta, \end{equation}

while, at the evaporation interface $z=0$,

(6.17b)\begin{equation} \frac{B\sqrt{D}}{\sqrt{{\bar{t}}}}\theta +D\theta_{z}=0,\end{equation}

and, as $z\rightarrow \infty$,

(6.17c)\begin{equation} \theta\rightarrow\theta_{IC}\exp\left({-}C\bar{t}\right),\quad R\rightarrow \frac{\theta_{IC}}{C}\left(1-\exp\left({-}C\bar{t}\right)\right). \end{equation}

Thus, under all of our simplifying assumptions, the equations for $\theta$ and $R$ decouple: we may solve the linear system (6.17a) (left), (6.17b) and (6.17c) (left) for $\theta$, before then computing $R$.

Specifically, we look for a solution for $\theta$ of the form

(6.18)\begin{equation} \theta=\theta_{IC}{\rm e}^{{-}C\bar{t}}F(\eta), \quad \text{where } \eta=\frac{z}{\sqrt{D\bar{t}}} \end{equation}

is a similarity variable. Thus, in the boundary layer the far-field deposition solution $\theta ^{out}=\theta _{IC}\textrm {e}^{-C\bar {t}}$ is modified by an ‘accumulation factor’ $F$, which describes how the evaporation causes suspended dirt to accumulate near the evaporating interface. Substituting this form of $\theta$ into (6.17a) (left equation), (6.17b) and (6.17c) (left equation), and solving the resulting ODE system for $F(\eta )$, we find that

(6.19) \begin{equation} \theta=\theta_{IC}{\rm e}^{{-}C\bar{t}}\Bigg(1+\frac{\sqrt{\rm \pi}B{\rm e}^{B^2}}{1-\sqrt{\rm \pi} B{\rm e}^{B^2}\,\text{erfc}(B)}\,\text{erfc}\Bigg(B+\frac{z}{2\sqrt{D\bar{t}}}\Bigg)\Bigg).\end{equation}

The factor $1-\sqrt {{\rm \pi} }B\textrm {e}^{B^2}\,\textrm {erfc}(B)$ is always positive (although it approaches zero as $B$, or equivalently $\sigma$, approaches $\infty$), and so the accumulation factor $F>1$ (and $F\rightarrow 1$ as $\eta$ (or $z$)$\rightarrow \infty$). Thus, as we would expect, the value of $\theta$ is higher in the boundary layer than in the far field.

The large $\theta$ in the boundary layer results in greater deposition occurring there, and so $R$ increases compared with the solution as $z\rightarrow \infty$. The problem (6.17a) (right equation) for $R$ is first-order hyberbolic, and – given the form (6.19) – may be solved using the method of characteristics. Specifically, we find the characteristic curves take the form $z=z_0-2B\sqrt {D\bar {t}},$ where $z_0$ parameterises the initial data $R=0$ at $\bar {t}=0$. (The shape of the characteristic curves mean we do not require data for $R$ at $\bar {z}=0$.) The solution $R$, in terms of $z$ and $\bar {t}$, is given by

(6.20)\begin{align} R(z,\bar{t})&=\frac{\theta_{IC}}{C}\left(1-{\rm e}^{{-}C\bar{t}}\right) +\frac{\theta_{IC}\sqrt{\rm \pi}B{\rm e}^{B^2}}{2C(1-\sqrt{\rm \pi}B{\rm e}^{B^2}\,\text{erfc}(B))}\nonumber\\ &\quad \times\left[\exp({-(\sqrt{C/D}z+2B\sqrt{C\bar{t}})})\,\text{erfc}\left(\frac{z}{2\sqrt{D\bar{t}}}+B- \sqrt{C\bar{t}}\right)\right.\nonumber\\ &\quad +\exp({(\sqrt{C/D}z+2B\sqrt{C\bar{t}})})\,\text{erfc}\left(\frac{z}{2\sqrt{D\bar{t}}}+B+\sqrt{C\bar{t}}\right)\nonumber\\ &\quad \left.-\,2{\rm e}^{{-}C\bar{t}}\,\text{erfc}\left(\frac{z}{2\sqrt{D\bar{t}}}+B\right)\right]. \end{align}

The solutions (6.19) and (6.20) of the paradigm model (6.17) are shown in figure 10. Clearly the system bears qualitative resemblance to the full model. The deposited dirt-layer thickness is maximised (spatially in the liquid–dirt domain, at any given time) at the evaporation front. Evaluating (6.20) at $z=0$, we have

(6.21)\begin{align} R|_{\bar{h}}&=\frac{\theta_{IC}}{C}\left(1-{\rm e}^{{-}C\bar{t}}\right) +\frac{\theta_{IC}\sqrt{\rm \pi}B{\rm e}^{B^2}}{2C(1-\sqrt{\rm \pi}B{\rm e}^{B^2}\,\text{erfc}(B))}\nonumber\\ &\quad \times\left[{\rm e}^{{-}2B\sqrt{C\bar{t}}}\,\text{erfc}\left(B-\sqrt{C\bar{t}}\right) -2{\rm e}^{{-}C\bar{t}}\,\text{erfc}\left(B\right) +{\rm e}^{2B\sqrt{C\bar{t}}}\,\text{erfc}\left(B+\sqrt{C\bar{t}}\right)\right]. \end{align}

In order to find the maximum value $R$ attains, we simply maximise $R|_{\bar {h}}$ over $\bar {t}$. We find that there is a single maximum for positive $\bar {t}$, at the critical time $\bar {t}_*=\tau ^2_*/C$, where $\tau _*=\tau _*(B)$ satisfies

(6.22)\begin{equation} \frac{2\tau_*}{\sqrt{\rm \pi}B^2}={\rm e}^{(B-\tau_*)^2}\,\text{erfc}(B-\tau_*)-{\rm e}^{(B+\tau_*)^2}\,\text{erfc}(B+\tau_*).\end{equation}

The solution $\tau _*$ of (6.22) is shown as a function of $B$ in figure 11(a). We see $O(1)$ variation of $\tau _*$ when $B$ varies over six orders of magnitude. Indeed, $\tau _*\rightarrow \sqrt {3/2}$ as $B\rightarrow \infty$ (or $\sigma \rightarrow \infty$, say, the limit of slow dirt diffusion), whilst $\tau _*$ grows very slowly as $B\rightarrow 0$, behaving like the growing root of $\tau _* \textrm {e}^{-\tau _*^2}=\sqrt {{\rm \pi} }B^2$ for $B\ll 1$, namely

(6.23)\begin{equation} \tau_*\sim\sqrt{\frac{-W_{{-}1}({-}2{\rm \pi} B^4)}{2}} \sim\sqrt{\log\left(\frac{1}{\sqrt{\rm \pi}B^2}\right)}+\frac{\log\left(\log\left(\dfrac{1} {\sqrt{\rm \pi}B^2}\right)\right)}{\sqrt{\log\left(\dfrac{1}{\sqrt{\rm \pi}B^2}\right)}}+\cdots, \end{equation}

where $W_{-1}$ is the lower branch of the real Lambert $W$ function.

Figure 10.

The black curves are the analytic solution (a) $\theta$ given by (6.19) and (b) $R$ for $Y>H$ given by (6.20), of the paradigm problem (6.17), at various times through the drying, with arrows showing increasing time. The red curve in (b) shows the final deposited dirt layer after the drying is completed. We have taken parameter values $\kappa =100$, $\nu =0.5$, $\sigma =1$, $\theta _{IC}=0.1$ and $r_0=0.4$.

Figure 11.

The maximum of $R$ for the paradigm problem. (a) The solution $\tau _*$ of (6.22) as a function of $B$, with the large $B$ limit $\sqrt {3/2}$ (red dashed) and the small-$B$ limit, namely the (non-negligible) root of $\tau _* \textrm {e}^{-\tau _*^2}=\sqrt {{\rm \pi} }B^2$ (green dashed). (b) The (scaled) maximum value of $R_{max}$, given by (6.25), against $B$. For small $B$, we have $G\sim 1$ while, for large $B$, we have $G\sim \textrm {e}^{-3/2}(\sqrt {6}-1) B^2$ (red dashed).

From the critical time $\bar {t}_*$, we can compute the value that $R$ attains at its maximum, namely

(6.24)\begin{equation} R_{max}=\frac{\theta_{IC}}{C}G(B),\end{equation}

where

(6.25)\begin{align} G(B)=1+\frac{1}{1-\sqrt{\rm \pi}B{\rm e}^{B^2}\,\text{erfc}(B)} \Big({\rm e}^{-\tau_*^2}(\tau_*-1)+\frac{\sqrt{\rm \pi}B{\rm e}^{B^2}}{2}{\rm e}^{2B\tau_*}\, \text{erfc}(B+\tau_*)\Big),\end{align}

with $\tau _*$ the solution of (6.22). We show $G=R_{max}C/\theta _{IC}$ as a function of $B$ in figure 11(b). For $B\ll 1$, we see from (6.25) that $G\approx 1$, while for $B\gg 1$, we find, using the large-$B$ value $\tau _*\approx \sqrt {3/2}$, along with the facts that $1-\sqrt {{\rm \pi} }B\textrm {e}^{B^2}\,\textrm {erfc}(B)\sim 1/(2B^2)$ and $\sqrt {{\rm \pi} }B\textrm {e}^{B^2}\textrm {e}^{\sqrt {6}B}\,\textrm {erfc}(B+\sqrt {3/2})\rightarrow \textrm {e}^{-3/2}$ as $B\rightarrow \infty$, that

(6.26)\begin{equation} G\sim {\rm e}^{{-}3/2}(\sqrt{6}-1) B^2 \quad \text{as }B\rightarrow\infty. \end{equation}

We note from (6.14a,b) that $B=O(\sqrt {\sigma })$, and so $R_{max}$ and $G$ increase linearly with $\sigma$ as $\sigma \rightarrow \infty$. Thus, we see that, no matter how small the initial suspended dirt level $\theta _{IC}$, for sufficiently slow dirt diffusion (sufficiently large $\sigma$), we will find that $R_{max}$ exceeds $R_{clog}$ and the system will clog.

Indeed, this expression (6.24) for $R_{max}$ gives the clogging condition for this paradigm setting: the system clogs when $R_{max}\geq R_{clog}$. Using (6.16), this clogging criterion is

(6.27)\begin{equation} \theta_{IC}\geq \theta_{IC}^{crit}:=\frac{\mathcal{C}_0\theta_*R_{clog}}{\phi_0G(B(\sigma, \nu))}.\end{equation}

We show this critical initial suspended dirt volume fraction (6.27) as a function of the suspended dirt diffusion rate $\sigma$ in figure 12, along with the small- and large-$\sigma$ limits (which follow directly from the small- and large-$B$ behaviour of $G$ discussed above), namely

(6.28a)$$\begin{gather} \theta_{IC}^{crit}\rightarrow \frac{\mathcal{C}_0\theta_*R_{clog}}{\phi_0} \quad\text{as }\sigma\rightarrow 0, \end{gather}$$

(6.28b)$$\begin{gather}\theta_{IC}^{crit}\rightarrow \frac{{\rm e}^{3/2}\nu\mathcal{C}_0\theta_*R_{clog}}{2\phi_0(\sqrt{6}-1) \log(1/(1-\nu))}\sigma^{{-}1} \quad\text{as }\sigma\rightarrow \infty. \end{gather}$$

We note that, since $G\geq 1$ and $R_{clog}>0$, this $\theta _{IC}^{crit}$, given by (6.27), for the paradigm problem is strictly smaller than the upper bound (6.8) (that essentially assumes a uniform deposited dirt layer), since, in the case $\beta =0$ that we have assumed for the paradigm case, the upper bound (6.8) may be rearranged to be written in terms of $\phi _0=1-{\rm \pi} r_0^2$, $\mathcal {C}_0=2{\rm \pi} r_0$ and $R_{clog}=1/2-r_0$, becoming

(6.29)\begin{equation} \theta_{IC}^{crit}=\frac{\theta_*(\mathcal{C}_0R_{clog}+{\rm \pi} R_{clog}^2)}{\phi_0}. \end{equation}

Thus, as expected, the non-uniformity of the dirt deposit in the boundary layer at the evaporation front results in clogging for lower initial suspended dirt volume fractions than if the dirt were to deposit uniformly.

Figure 12.

The critical initial suspended dirt volume fraction (6.27) predicted in the paradigm case, as a function of $\sigma$. Dashed curves are the low- and high-$\sigma$ limits (6.28). We take parameter values $r_0=0.45$, $\nu =0.1$ and $\theta _*=1$.

We note that our expression (6.27) for $\theta _{IC}^{crit}$ in this paradigm case is independent of $\kappa$, since we have taken the leading-order behaviour as $\kappa \rightarrow \infty$. Since we require the boundary layer structure, our paradigm estimate for the clogging criterion is only valid when the width $1/\sqrt {\sigma \kappa }$ of the boundary layer is small, and so only for $\sigma$ sufficiently large that $\sigma \gg 1/\kappa$ (although since we assume $\kappa \gg 1$ this is not particularly restrictive).

Additionally, we see that the expression (6.27) is independent of $D$, and hence, of $\mathcal {D}_0$. From this, we learn that it is the relative diffusivity of the suspended dirt and the liquid vapour, captured through $\sigma$ that is important, and not how both are equally affected by the presence of the porous material (recall that $\mathcal {D}_0$ is the effective diffusivity). Of course, the depth $y$ into the porous material at which the clogging occurs does depend on $D$. Specifically, the clogging depth when $\theta _{IC}=\theta _{IC}^{crit}$ takes its critical value is

(6.30)\begin{equation} H_{clog}^{crit}=\sqrt{\frac{2\mathcal{D}_0\log(1/(1-\nu))}{\nu\mathcal{C}_0\theta_*}}\frac{\tau_*}{\sqrt{\kappa}}. \end{equation}

Curiously, although this position given by (6.30) depends on $\kappa$ and $\mathcal {D}_0$, it is only weakly dependent on $\sigma$, since $\tau _*$ (the solution of (6.22)) was seen to have very weak dependence on $B\propto \sqrt {\sigma }$.

We considered a simplified, paradigm version of the problem in this section, in order to make analytic progress. Although this analysis does not directly relate to the full drying model, we may extrapolate some general conclusions. Firstly, the qualitative wet-clogging behaviour seen in the full model, characterised by an internal peak in $R$, does not rely on the variation of $\mathcal {D}$, $\phi$ and the evaporation-front speed $h_T$ with $R$ and $\theta$ (in the paradigm case we supposed all these were constant). Instead, the important mechanism, captured by the paradigm model, is the variation of the rate of dirt deposition with $\theta$, along with the fact that $\theta$ is spatially non-uniform in the boundary layer at the evaporation front, determined by a balance of all three mechanisms of diffusion, accumulation due to liquid evaporation and the deposition.

In the following section we generate bifurcation diagrams showing parameter regimes for which the system clogs. We compare these numerical results with the predictions of the paradigm model as appropriate.