A homogenised model for the motion of evaporating fronts in porous media

Abstract

Evaporation within porous media is both a multiscale and interface-driven process, since the phase change at the evaporating interfaces within the pores generates a vapour flow and depends on the transport of vapour through the porous medium. While homogenised models of flow and chemical transport in porous media allow multiscale processes to be modelled efficiently, it is not clear how the multiscale effects impact the interface conditions required for these homogenised models. In this paper, we derive a homogenised model, including effective interface conditions, for the motion of an evaporation front through a porous medium, using a combined homogenisation and boundary layer analysis. This analysis extends previous work for a purely diffusive problem to include both gas flow and the advective–diffusive transport of material. We investigate the effect that different microscale models describing the chemistry of the evaporation have on the homogenised interface conditions. In particular, we identify a new effective parameter, $\mathcal{L}$, the average microscale interface length, which modifies the effective evaporation rate in the homogenised model. Like the effective diffusivity and permeability of a porous medium, $\mathcal{L}$ may be found by solving a periodic cell problem on the microscale. We also show that the different microscale models of the interface chemistry result in fundamentally different fine-scale behaviour at, and near, the interface.

1. Introduction

The drying of a porous or granular material is of great importance in the construction and ceramics industries [22, 23], in soil mechanics [13, 20], in the food manufacturing industry [1, 2], and in the production of pharmaceuticals [28]. In the case of evaporation from soils and hygroscopic porous materials (in which moisture is bound within the solid matrix), the evaporation process may be complicated by the release of moisture adsorbed by the grains even when the large pores of the material, between grains, are no longer saturated with liquid [3, 22]. The evaporation process is simpler when the porous material consists of an impermeable solid structure, such as an array of glass beads [22, 26].

Liquid and gas may be distributed in multiple ways within the pores of a porous medium. For example, both liquid and gas phases may co-exist within the pore space, or the two phases may occupy distinct macroscale regions of the porous medium, separated by an evaporating interface. Detailed modelling of mass, momentum, and heat transfer effects, and a derivation of equations averaged over the microscale in the former case is undertaken in [29], while the latter is studied in, for instance, [27].

In fact, different configurations of liquid and gas within the pore space can occur at different stages of the evaporating process [10, 20]. Liquid that saturates a porous medium initially evaporates from the surface of the porous medium. As liquid is removed a region containing both liquid and gas appears between the surface of the porous medium and the region saturated with liquid, and liquid is drawn to the surface by capillary action. This stage of drying is known as the ‘constant rate period’, as the evaporation rate from the surface is approximately constant. Eventually, the connections between the liquid reservoir and the surface of the porous medium break, and evaporation then occurs within the pores of the porous medium, at the interfaces between the gas and liquid phases. The evaporation process is then affected by the rate at which vapour is transported out of the porous material, which decreases as the gas–liquid interface advances further into the porous medium, and so this stage of drying is known as the ‘falling rate period’. Capillary effects may still be important during the falling rate period [30, 31]. However, in the case of low surface tension, capillary effects are very small, and we expect to find a sharp interface between the liquid- and gas-occupied regions. Our focus in this paper will be the case in which a macroscale evaporation front divides a region of pore space saturated with liquid from a region saturated with gas.

The motion of a one-dimensional evaporation front through porous media is studied in [27], assuming a priori a volume-averaged or homogenised model consisting of Darcy flow and the advection–diffusion of chemicals. A similarity solution is found describing the motion of the evaporating interface, limited by the transport of vapour through the porous medium. The model does not explicitly take the microstructure of the porous medium into account. Indeed, the authors note ‘It is certainly acknowledged that pore structure will affect diffusion, the formulation of the problem, and the motion and dimensionality of the interface. However, these effects are not in the scope of this investigation’.

The interaction of microscale and macroscale processes makes evaporation within porous media complicated and computationally expensive to simulate. Simulations on relatively small areas of porous medium have been performed, for instance in [10]. Pore-network approaches have been taken as a simplification of continuum models, including in [18, 26]; the approach is reviewed in [22]. In this method, the porous structure is modelled as a network with the nodes as pores and the edges as the narrow throats between pores. Another approach is to view evaporation as an ‘immiscible displacement’ problem [24]; here, the gas phase is injected at high pressure and forces evaporation of the liquid (a technique used in oil recovery [27]). The stability of such evaporation fronts, due to the transport of vapour away from the front, may be demonstrated experimentally [24].

In numerous applications, homogenisation and other upscaling methods have been used to great effect in deriving effective partial differential equations (PDEs) to describe flow and transport processes in porous media. For instance, Darcy’s law for fluid flow through a porous medium may be derived by homogenising the Stokes flow equations within the pores [14]. The equations describing the advection and diffusion of chemicals, particles or heat through a porous medium may be homogenised to give effective advection–diffusion equations, in terms of the Darcy velocity [5, 6]. A homogenised model for a multi-phase flow is derived in [8], while flows in deformable porous media are homogenised in [25]. In each case, the result of the homogenisation analysis is a set of averaged, or effective, PDEs that hold over the domain of the entire porous medium, which are significantly simpler to solve than the original model (which is valid within the complicated domain of the pore space) while rigorously taking the microstructure of the porous medium into account.

However, evaporation within a porous medium is an example of a multiscale process that is driven by processes at an interface, in this case, the evaporating interface separating the gas- and liquid-filled regions of the porous medium. Another example is the reactive decontamination of porous media, as studied in [7, 16]. Here, a chemical contaminant within the pore space is neutralised by a cleanser fluid, with the neutralisation reaction occurring at the interface between the two immiscible fluids. Stefan problems within porous media also fall into this category of interface-driven behaviour, for example, those involving the freezing or melting of a liquid within a porous structure [11, 15, 17]. While the homogenised PDEs for the macroscale flow and transport of material through porous media incorporate the microscale processes and structure through effective parameters, it is not clear how and if the microstructure affects the interface conditions. The key question is: what interface conditions should be imposed on the homogenised PDEs that accurately capture the effect the microstructure has on the interface processes? This question was addressed for the simple reactive decontamination process described in [16], where effective interface conditions at the reacting interface are derived using a combined homogenisation and boundary layer approach. A key result is that, in order to smooth out temporal oscillations in the diffusive flux at the interface (due to the solid obstructions on the pore-scale), a nested boundary layer structure is required.

To describe an evaporating process, we must understand the mechanism by which liquid molecules become vapour at the evaporating interface. A commonly used concept is that the evaporative flux is determined by the difference between the local rate of evaporation and condensation of liquid molecules; this is modelled using the Hertz–Knudsen law [12]. There are some concerns with this model, in particular since the evaporation and condensation coefficients must be determined experimentally [21], but it is nevertheless widely used. Some authors (see [9], for example) instead assume that the surface is in thermodynamic equilibrium, with the chemical potential continuous across the interface. Thus, equivalently, the vapour density is at the saturation density. We note that this is a sublimit of the Hertz–Knudsen law in the case where the rate at which the vapour pressure at the interface reaches the saturation pressure is much faster than the rate of interface motion. For slowly moving evaporation fronts, the Dirichlet condition where the vapour density is at its prescribed saturation point is therefore likely to be accurate.

In this paper, we aim to use homogenisation methods combined with boundary layer analysis to derive an effective model for the motion of an evaporation front through a porous medium consisting of an impermeable, non-reactive, solid matrix, including effective interface conditions at the evaporating boundary. We assume for simplicity that the liquid is stationary, and that there is no variation in temperature. Our evaporation model will describe the flow of gas and the transport of the liquid vapour through an inert gas, away from the gas–liquid interface. Our analysis will build on the method of [16], which was for a purely diffusive system, incorporating the flow of the gas mixture and the advective–diffusive transport of liquid vapour. In particular, we will investigate the effect of different assumptions about the surface chemistry at the evaporating interface. We will use a Hertz–Knudsen law but will also consider the sublimit where the vapour at the interface is at its saturation point. Indeed, we will find that these two scenarios give very different fine-scale behaviour near the evaporating interface.

The remainder of the paper is structured as follows. In Section 2, we introduce the model for the flow and transport of vapour in the gas-filled pores of the porous medium, and nondimensionalise, taking care to ensure the correct balance of terms on the microscale. In Section 3, we then homogenise this model to derive the effective PDEs and interface conditions at the macroscale evaporating interface. The analysis in Section 3 does not require us to specify the chemistry condition at the interface. However, an effective chemistry condition is required to close the homogenised model; we derive this in Section 4, in two sublimits of the Hertz–Knudsen relation. We then explore how the different chemistry boundary conditions result in very different asymptotic structures in the boundary layer(s) around the evaporating interface. Discussion of the homogenisation process and the results of our analysis are given in Section 5.

2. The pore-scale model

In this section, we present a model for the evaporation of a liquid from within the pores of a porous medium. We consider a two-dimensional, locally periodic porous medium, initially saturated with stationary liquid, occupying the region $y\leq 0$ , with the surface of the porous medium at $y=0$ . As the liquid begins to evaporate, we assume that there is an evaporation front, located at $y=h(x,t)$ for $h\leq 0$ , moving down into the porous medium in the direction $-y$ . This front $y=h(x,t)$ is the interface or interfaces between the gas and liquid phase: we note that due to the solid microstructure, the front may be formed of several disjoint interfaces, and so h(x, t) may be defined piecewise. Above the evaporation front, in $y\in[h,0]$ , there is a mixture of inert gas (eg: air) and liquid vapour. We make the simplifying assumption that the liquid is stationary so that there is no capillary action or any other effects due to viscous or surface tension forces. This no-liquid flow case is valid in the limit of very small forces due to surface tension compared with liquid inertia. Although the liquid is stationary, we study the flow of the gas mixture and the advection–diffusion of the vapour through the pores, as well as the motion of the evaporating front down through the pores. A schematic is shown in Figure 1.

Figure 1. Schematic showing the regions of pore space saturated with liquid (blue) and gas mixture (yellow). The solid microstructure is illustrated by the grey circles. The pore lengthscale, d, and macro-lengthscale, l, are shown; in reality, $\epsilon=d/l$ would be much smaller than in the schematic.

We assume that the total density $\rho_G$ of the gaseous mixture is constant so that the gas mixture is incompressible. In the gas-occupied pore space, we estimate a Reynolds number of around $10^{-4}$ (taking the approximate values for water evaporating into air $\rho_G\approx 1\,$ kg $\,\textit{m}^{-3}$ and the kinematic viscosity $\mu\approx 10^{-5}\,$ kg $\,\textit{m}^{-1}\textit{s}^{-1}$ , using the pore lengthscale $d\approx 10^{-6}\,$ m typical for brick or concrete, and using the approximation $[u]=D/l$ for the velocity scaling introduced in Section 2.1 below, where the diffusivity of water vapour in air is around $D\approx 10^{-5}\,\textit{m}^2\,\textit{s}^{-1}$ , and we assume the evaporation lengthscale is around $l\approx 10^{-2}\,$ m). Since $\text{Re}\ll1$ , it is reasonable to assume that the gas mixture satisfies the Stokes flow equations:

(2.1a) \begin{align} \nabla\cdot{\textbf{u}}=0, \end{align}

(2.1b) \begin{align} \mu\nabla^2{\textbf{u}}-\nabla p={\textbf{0}}, \end{align}

within the pores, where ${\textbf{u}}$ , p and $\mu$ are the velocity, pressure and viscosity of the gas mixture, respectively. We assume Fickian diffusion of liquid vapour through the inert gas so that the liquid vapour, with local density $\rho_v$ , is transported within the mixture according to the advection–diffusion equation:

(2.1c) \begin{equation} (\rho_v)_t+{\textbf{u}}\cdot\nabla\rho_v=D\nabla^2\rho_v,\end{equation}

where D is the diffusivity of the vapour in the inert gas, and the subscript t denotes the partial derivative with respect to time.

At interfaces between the gas mixture and solid microstructure, we assume that there is no flux of the inert gas or vapour into the solid material, and no slip on the solid surface, so that

(2.2a) \begin{align} {\textbf{u}}={\textbf{0}}, \end{align}

(2.2b) \begin{align} D\nabla\rho_v\cdot{\textbf{n}}_s=0, \end{align}

where ${\textbf{n}}_s$ is the normal to these interfaces.

At interfaces between the gas mixture and the stationary liquid, we impose the total conservation of mass of both liquid and gas over the interface, and also the conservation of mass of the liquid molecules across the interface, which may be expressed as:

(2.3a) \begin{align} -\rho_LV_m&=\rho_G({\textbf{u}}\cdot{\textbf{m}}-V_m), \end{align}

(2.3b) \begin{align} -\rho_LV_m&=\rho_v({\textbf{u}}\cdot{\textbf{m}}-V_m)-D\nabla\rho_v\cdot{\textbf{m}}, \end{align}

on $y=h(x,t)$ , where ${\textbf{m}}$ is normal to the gas–liquid interface (pointing into the gas region) and $V_m$ is the velocity of the interface in the normal direction ${\textbf{m}}$ . The (constant) density of the liquid phase is $\rho_L$ . The mixture velocity ${\textbf{u}}=u{\textbf{e}}_x+v{\textbf{e}}_y$ also satisfies a kinematic condition on the interface, which reduces to

(2.3c) \begin{equation} u+vh_x=0,\end{equation}

on $y=h(x,t)$ , in the case of interest where the liquid is stationary. We note that the normal velocity of the interface and unit normal may be expressed as:

(2.4) \begin{equation} V_m=\frac{h_t}{\sqrt{1+h_x^2}}, \qquad {\textbf{m}}=\frac{(h_x,-1)}{\sqrt{1+h_x^2}}.\end{equation}

The interface h may be defined piecewise as necessary due to the solid microstructure of the porous medium.

The Stokes flow and advection–diffusion model presented thus far are very general – a statement of the physical principles of conservation of mass and momentum – and could apply to a variety of different processes. To close the system and determine the motion of the gas–liquid interface, we need to impose an additional condition at the interface that determines the evaporation rate. This extra condition encapsulates the chemistry driving the evolution of the system. As discussed in Section 1, we use a Hertz–Knudsen law, which, assuming a constant temperature, states that the mass flux J across the interface, is proportional to the difference between the vapour density at the interface and the (constant) saturation density of the liquid vapour, denoted $\rho_v^*$ , and thus we write

(2.5) \begin{equation} J=a(\rho_v^*-\rho_v).\end{equation}

Here, a [m $\,\textit{s}^{-1}$ ] is a mass transfer coefficient, the value of which will depend in general on the temperature, as well as the chemistry of the fluids. We will assume that a is constant. We note that the mass flux across the interface may be given in terms of the interface velocity, according to

(2.6) \begin{equation} J=-\rho_LV_m.\end{equation}

2.1. Nondimensionalisation

We nondimensionalise our pore-scale model (2.1)–(2.6) using the scalings:

(2.7) \begin{align} \textbf{x}=d\hat{\textbf{x}}, \quad h=d\hat{h}, \quad \rho_v=\rho_v^*\hat{\rho}, \quad t=\frac{d^2}{\epsilon\delta D}\hat{t}, \quad \textbf{u}=\frac{D\nu\epsilon}{d}\hat{\textbf{u}}, \quad p=\frac{\mu D \nu\epsilon^2}{d^2}\hat{p},\end{align}

where we have introduced the dimensionless parameters:

(2.8) \begin{equation} \delta=\frac{\rho_v^*}{\rho_L}, \qquad \epsilon=\frac{d}{l}\ll1, \qquad \nu=\frac{\rho_v^*}{\rho_G}.\end{equation}

These are, respectively, the ratio of the saturated vapour density to the liquid density, the ratio of the pore lengthscale, d, to the lengthscale of the drying porous medium, l, (these lengths are shown schematically in Figure 1), and the ratio of the saturated vapour density to the gas mixture density (the dilution of the vapour within the gas mixture).

We have chosen the lengthscale to be d, the lengthscale of the pore space within the porous medium, and the vapour density scaling to be the saturation vapour density, $\rho_v^*$ . The timescale is chosen to be that over which the liquid–gas interfaces move over the pore lengthscale d. To see this, we note that the motion of these microscale interfaces is due to transport of the vapour away from the interfaces. Since the vapour is much less dense than the liquid, we will see that evaporation drives a gas flow moving away from the interface. The lengthscale of the vapour transport, that is, that over which $\hat{\rho}$ has O(1) variation, must therefore be l, the lengthscale of the entire porous medium. The timescale has therefore been chosen to balance the diffusive flux of vapour out of the interface with the left-hand side of (2.3b). The velocity scaling for the gas is determined by balancing the terms in (2.3a), and the pressure scaling is chosen to give the correct balance in order to derive Darcy’s law, as discussed in [6].

Substituting from (2.7) into (2.1)–(2.6), and dropping the hat notation for dimensionless variables, we find that for $y>h(x,t)$ in the pore space:

(2.9a) \begin{align} \epsilon\nabla\cdot{\textbf{u}}&=0, \end{align}

(2.9b) \begin{align} \epsilon\nabla^2{\textbf{u}}&=\nabla p, \end{align}

(2.9c) \begin{align} \epsilon\left(\delta\rho_t+\nu{\textbf{u}}\cdot\nabla\rho\right)&=\nabla^2\rho, \end{align}

and at gas–solid interfaces,

(2.9d) \begin{equation} {\textbf{u}}={\textbf{0}}, \qquad \nabla\rho\cdot{\textbf{n}}_s=0,\end{equation}

while at $y=h(x,t)$ ,

(2.9e) \begin{align} {\textbf{u}}\cdot{\textbf{m}}-\delta \frac{h_t}{\sqrt{1+h_x^2}}&=-\frac{h_t}{\sqrt{1+h_x^2}}, \end{align}

(2.9f) \begin{align} \epsilon \nu\rho\left({\textbf{u}}\cdot{\textbf{m}}-\delta \frac{h_t}{\sqrt{1+h_x^2}}\right)-\nabla\rho\cdot{\textbf{m}}&=-\epsilon\frac{h_t}{\sqrt{1+h_x^2}}, \end{align}

(2.9g) \begin{align} u+vh_x&=0, \end{align}

along with the chemistry condition:

(2.9h) \begin{align} 1-\rho&=-\alpha\frac{h_t}{\sqrt{1+h_x^2}}.\end{align}

Here, the dimensionless parameter

(2.10) \begin{equation} \alpha\;:\!=\;\frac{\rho_L[x]}{a[t]\rho_v^*}=\frac{D}{a l}\end{equation}

describes the ratio of the speed of the interface motion to the rate at which $\rho$ approaches its equilibrium value (the saturation density). We note that $\alpha$ decreases as the lengthscale l of the porous material increases.

There are four dimensionless parameters in the model (2.9), namely $\epsilon$ , $\delta$ , $\nu$ and $\alpha$ . We note that all of $\delta,\epsilon$ and $\nu$ are less than one by definition. We require $\epsilon\ll 1$ in order to perform the homogenisation analysis. This is a good assumption since pores in concrete or brick, say, are on the order of microns, whereas we might be interested in the drying of centimetres of the porous medium. We expect that $\delta\approx 10^{-3}$ , which is also small (this value is for a water and air combination). The parameter $\nu$ gives the dilution of the liquid vapour in the gas mixture, and may be small, but is not expected to be as small as $\delta$ . For the homogenisation analysis in the following section, we assume that $\delta$ and $\nu$ are both O(1), whilst $\epsilon\ll1$ . We may subsequently take either of the sublimits $\delta\rightarrow 0$ or $\nu\rightarrow 0$ after the homogenisation, if necessary. We note that the reduced version of the homogenised model in the case(s) $\delta\ll1$ and/or $\nu\ll 1$ remains valid in the case $\delta\ll\epsilon$ or $\nu\ll\epsilon$ .

The final parameter, $\alpha$ , quantifies how far the vapour density at the interface is from the saturation density or, alternatively, how far the interface is from being in thermodynamic equilibrium. We note that we have assumed for our nondimensionalisation that $\alpha$ is at most O(1), or equivalently, that the lengthscale l of the porous material is sufficiently large. If $\alpha\gg 1$ then we would require a different, longer, timescale of the interface motion, $t\sim\alpha$ , by balancing (2.9h), since the evaporation process is in that case limited by the rate at which liquid molecules escape into the gas phase. Conversely, in the case that $\alpha$ is small, the evaporating interface moves slowly enough that equilibrium is reached at the surface, and so the chemistry interface condition (2.9h) reduces to a Dirichlet condition, $\rho=1$ . This statement of ‘continuous chemical potential’ over the interface is contained within the Hertz–Knudsen law in the $\alpha\ll1$ limit. In Section 4, we will investigate in detail how the size of $\alpha$ (relative to powers of $\epsilon$ ) affects both the effective interface conditions for the homogenised model and the fine-scale processes near the evaporating interface.

3. Homogenisation

In this section, we homogenise the pore-scale model (2.9), deriving effective PDEs that hold over the entire porous medium (rather than only within the pore space) and also the effective boundary conditions that hold at the gas–liquid interface. Although the microscale behaviour is necessarily two-dimensional due to the pore structure, for the homogenisation analysis we restrict our attention to the simplest case in which the macroscale behaviour is one-dimensional, with a flat macroscale evaporating interface, parallel to the surface of the porous material. We assume that the microstructure of the porous medium is spatially periodic.

After giving an overview of the homogenisation process in Section 3.1, we will then derive an effective model for the gas flow in Section 3.2, before analysing the advection–diffusion system in Section 3.3. We note that for neither of these analyses do we need to use the form of the chemistry interface condition, which will be incorporated in Section 4, for various ranges of $\alpha$ .

3.1. Overview and problem set-up

We perform the homogenisation in terms of the small parameter $\epsilon$ and assume the distinguished limit in which $\delta,\nu=O(1)$ . Since our dimensionless variables ${\textbf{x}}=(x,y)$ and t are O(1) on the microscale, we also introduce macroscale variables:

(3.1) \begin{equation} {\textbf{X}}\;:\!=\;\epsilon {\textbf{x}}, \qquad T\;:\!=\;\epsilon t,\end{equation}

in which the lengthscale is O(1) on the entire porous medium, and the macroscale timescale is that of the interface motion over the entire porous medium. We change variables to these macroscale variables X, Y, T and denote the macroscale interface $Y=H(T)$ , where H(T) is the average interface position and is to be determined as part of the solution. We note that H(T) is independent of the microscale (since it is an average value). We assume that the macroscale problem is one-dimensional, so the interface is flat on the macroscale with H independent of X. We then employ a second change of variables to the frame moving with the macroscale interface, by setting $Z=Y-H(T)$ .

To average our model, we perform a homogenisation analysis in both space and time. We will find that, far from the gas–liquid interface, the fluid flow problem homogenises to Darcy flow, and we will find an effective advection–diffusion equation that holds on the macroscale. This homogenisation procedure is almost exactly as presented in [6], although we must take care with the time dependence, since we have changed to a travelling coordinate system and the solid microstructure is continuously moving up through the microscale domain. Near the macroscale interface, we combine a homogenisation analysis with a boundary layer analysis, in a similar way to that presented in [16], in order to incorporate effects due to the microstructure into the effective interface conditions. We refer to the homogenisation away from the interface as the outer problem, and the processes within an $O(\epsilon)$ distance of the interface as the inner problem. As in [16], we will find that, in general, an additional intermediate boundary layer is required, with lengthscale $O(\epsilon^{1/2})$ , over which temporal oscillations in the vapour density (due to the oscillations in the evaporation rate at the interface) are smoothed out, allowing the inner problem to be matched with the outer problem. These regions are shown in the schematic of Figure 2. We state the problem in the outer and inner regions in Sections 3.1.1 and 3.1.2 below but introduce the intermediate problem in Section 3.3, since it is only needed for the advection–diffusion problem.

Figure 2. Schematic showing the macroscale evaporating interface, at $Z=0$ in the moving coordinate system, with the outer, intermediate and inner regions for the boundary layer analysis.

The asymptotic analysis presented in the following sections forms a general framework for an evaporation process. For certain ranges of $\alpha$ , we will show (see Section 4) that this framework may be significantly simplified: when $\alpha \gtrsim O(\epsilon^{1/2})$ the intermediate layer at the evaporating interface is necessary to smooth oscillations in the vapour flux, while for $\alpha\ll\epsilon^{1/2}$ , the intermediate region is shown to be redundant and only a single boundary layer (the inner region) about the evaporating front is required.

3.1.1. Statement of outer problem

Away from the evaporating interface (an O(1) macro-lengthscale distance away), we (re)-introduce microscale variables $\xi=X/\epsilon$ and $\eta=Z/\epsilon$ (these move with the macroscale interface, hence differ from the static coordinates x and y that we had previously) and work on the timescale $t=T/\epsilon$ . We introduce a microscale cell containing a periodic section of the microstructure, denoted by $\omega_c$ , and illustrated in Figure 3 (left), with fluid-occupied region $\omega_f(t)$ , and solid inclusion $\omega_s(t)$ . We treat the dependent variables $\xi,\eta,t$ , and Z and T as independent variables and so (2.9a)–(2.9c) in the outer region become

(3.2a) \begin{align} \epsilon\left(\nabla_\xi+\epsilon{\textbf{e}}_2\frac{\partial}{\partial Z}\right)\cdot{\textbf{u}}&=0, \end{align}

(3.2b) \begin{align} \epsilon\left(\nabla^2_\xi+\epsilon\left(\nabla_\xi\cdot{\textbf{e}}_2\frac{\partial}{\partial Z}+{\textbf{e}}_2\frac{\partial}{\partial Z}\cdot\nabla_\xi\right)+\epsilon^2\frac{\partial^2}{\partial Z^2}\right){\textbf{u}}&=\left(\nabla_\xi+\epsilon{\textbf{e}}_2\frac{\partial}{\partial Z}\right)p,\end{align}

(3.2c) \begin{align} \epsilon\delta\left(\rho_t+\epsilon\rho_T-H_T\left(\rho_\eta+\epsilon\rho_Z\right)\right)+\left(\nabla_\xi+\epsilon{\textbf{e}}_2\frac{\partial}{\partial Z}\right)&\cdot\left(\epsilon\nu\rho{\textbf{u}}-\left(\nabla_\xi+\epsilon{\textbf{e}}_2\frac{\partial}{\partial Z}\right)\rho\right)=0.\end{align}

The vapour density $\rho$ , gas velocity ${\textbf{u}}$ and pressure p in this outer region are all functions of $\xi,$ $\eta,$ t, Z, and T. On the boundary of the cell $\omega_c$ , we assume that the dependent variables are periodic, with period 1, in $\xi,\eta$ . We note that the geometry of the microscale cell changes over time, since the solid inclusions move up through the cell at speed $H_T$ . Since the microstructure is spatially periodic, this means that the cell geometry is also periodic in t, with period $1/H_T$ . We also therefore assume that the dependent variables are periodic in t, with period $1/H_T$ . On the surface $\partial\omega_s(t)$ of the solid material (2.9d) becomes

(3.2d) \begin{align} {\textbf{u}}={\textbf{0}}, \end{align}

(3.2e) \begin{align} \left(\nabla_\xi+\epsilon{\textbf{e}}_2\frac{\partial}{\partial Z}\right)\rho\cdot{\textbf{n}}_s=0, \end{align}

where ${\textbf{n}}_s(t)$ is the unit normal to the solid surface. In (3.2), we have introduced the notation:

(3.3) \begin{equation} \nabla_\xi={\textbf{e}}_1\frac{\partial}{\partial \xi}+{\textbf{e}}_2\frac{\partial}{\partial \eta},\end{equation}

for the microscale gradient operator, where ${\textbf{e}}_1$ and ${\textbf{e}}_2$ are the unit vectors in the X or $\xi$ , and Z or $\eta$ directions, respectively (we use the same notation for unit vectors in both coordinate systems).

Figure 3. Schematics of the microscale cells regions. Left: the cell $\omega_f(t)=[-1/2,1/2]^2\setminus \omega_s(t)$ for the outer and intermediate regions. Right: the inner cell $\omega_{{\mathcal{H}}}$ formed of the region between $z={{\mathcal{H}}}(x,t)$ and $\Gamma$ , $x\in[-1/2,1/2]$ , excluding the solid structure $\omega_s(t)$ . In the analysis we consider the limit as the height of $\Gamma$ increases, so that $\omega_{{\mathcal{H}}}$ becomes semi-infinite in z. The solid inclusions $\omega_s(t)$ (shown as grey circles) move up through both cells at speed $H_T$ , so that the cells are periodic in time t with period $1/H_T$ .

3.1.2. Statement of the inner problem

As in [16], we will find that there must be a boundary layer near the gas–liquid interface in which oscillations due to the motion of the interface are felt at leading order. In our moving coordinate system, the macroscale interface is at $Z=0$ . We look for a boundary layer of the macroscale problem by changing to the microscale variables $Z=z/\epsilon$ , $X=x/\epsilon$ centred on the interface. We define the interface to be located at $z={\mathcal{H}}$ on the microscale, where $\mathcal{H}(x,t)$ is the microscale variation of the evaporating interface about its mean position $z=0$ . We assume the problem is periodic in x and t, and so study the inner region on the microscale ‘cell’ region, $\omega_{\mathcal{H}}(t)$ , illustrated on the right of Figure 3. This cell is the region $x\in[-1/2,1/2]$ between $z={{\mathcal{H}}}(x,t)$ and some curve $\Gamma$ , excluding the solid structure $\omega_s(t)$ . However, we cannot assume periodicity in z. Instead, we will match with the outer region in the ${\textbf{e}}_2$ direction by sending the height z of the curve $\Gamma$ to $\infty$ so that the inner cell approaches a semi-infinite region. We call this region a ‘cell’ by analogy with the cell region $\omega_c(t)$ for the outer problem. We will homogenise along the macroscale interface in the x direction and in time t.

In this inner region, we use an overbar for the dependent variables $\bar{\rho}$ , $\bar{p}$ and $\bar{{\textbf{u}}}$ , to differentiate from the variables in the outer region. These inner region variables are functions of x, z, t and T (via the value of $H_T$ ) but are independent of Z. We also use the gradient operator:

(3.4) \begin{equation} \nabla_x={\textbf{e}}_1\frac{\partial}{\partial x}+{\textbf{e}}_2\frac{\partial}{\partial z}.\end{equation}

The equations (2.9a)–(2.9c) in $\omega_{{\mathcal{H}}}(t)$ are

(3.5a) \begin{align} \epsilon\nabla_x\cdot\bar{{\textbf{u}}}&=0, \end{align}

(3.5b) \begin{align} \epsilon\nabla_x^2\bar{{\textbf{u}}}&=\nabla_x \bar{p}, \end{align}

(3.5c) \begin{align} \epsilon\delta(\bar{\rho}_t+\epsilon\bar{\rho}_T-H_T\bar{\rho}_z)+\nabla_x\cdot\left(\epsilon\nu\bar{\rho}\bar{{\textbf{u}}}-\nabla_x\bar{\rho}\right)&=0.\end{align}

On the solid surface $\partial\omega_s(t)$ , (2.9d) becomes

(3.5d) \begin{align} \bar{{\textbf{u}}}={\textbf{0}}, \end{align}

(3.5e) \begin{align} \nabla_x\bar{\rho}\cdot{\textbf{n}}_s=0,\end{align}

and (2.9e)–(2.9g) become, on $z={{\mathcal{H}}}(x,t)$ ,

(3.5f) \begin{align} \bar{{\textbf{u}}}\cdot{\textbf{m}}&=-(1-\delta)\frac{H_T+{{\mathcal{H}}}_t}{\sqrt{1+{{\mathcal{H}}}_x^2}},\end{align}

(3.5g) \begin{align} \nabla_x\bar{\rho}\cdot{\textbf{m}}&=\epsilon(1-\nu\bar{\rho})\frac{H_T+{{\mathcal{H}}}_t}{\sqrt{1+{{\mathcal{H}}}_x^2}}, \end{align}

(3.5h) \begin{align} \bar{u}+\bar{v}{{\mathcal{H}}}_x&=0.\end{align}

The chemistry interface condition, (2.9h), applies at the microscale interface $z={{\mathcal{H}}}$ . Although we will not use this condition until Section 4, it reads

(3.6) \begin{equation} 1-\bar{\rho}=-\alpha \left(\frac{H_T+{{\mathcal{H}}}_t}{\sqrt{1+{{\mathcal{H}}}_x^2}}\right) \qquad \text{ on }z={{\mathcal{H}}}(x,t).\end{equation}

3.2. Gas flow problem

We note that the equations for the gas velocity and pressure, namely (3.2a)–(3.2b) and (3.2d) in the outer region, and (3.5a)–(3.5b) with (3.5d), (3.5f) and (3.5h) in the inner region, decouple from the advection–diffusion problem for the vapour density. We therefore begin by analysing this gas flow problem (assuming that we know $H_T$ and ${{\mathcal{H}}}$ ). We start with the outer region, in which the Stokes flow equations are homogenised to Darcy’s law, via an analysis that is largely standard, although we take care to understand how temporal oscillations in the dependent variables are accounted for. We then consider the inner region of the domain and derive an effective flux boundary condition for the homogenised Darcy’s law, using boundary layer and homogenisation analysis.

3.2.1. Outer region

The homogenisation of Stokes flow to Darcy’s law has been performed in, for instance, [6, 14]. We follow the analysis presented in [6] closely, but note that, in our evaporation setting, the velocity and the pressure gradients in general all vary in time t as well as in space. We will therefore need to average the dependent variables in time as well, although this does not significantly change the analysis or results of the homogenisation.

We expand the dependent variables ${\textbf{u}}$ and p in powers of $\epsilon$ , using superscript notation to denote the order, for example, $ p=p^{(0)}+\epsilon p^{(1)}+\dots$ . At O(1) we find that (3.2b) becomes

(3.7) \begin{equation} \nabla_\xi p^{(0)}={\textbf{0}},\end{equation}

so that $p^{(0)}$ is independent of the microscale spatial variables $\xi$ and $\eta$ . At $O(\epsilon)$ we find that (3.2a)–(3.2b) become

(3.8a) \begin{align} \nabla_\xi\cdot{\textbf{u}}^{(0)}&=0, \end{align}

(3.8b) \begin{align} \nabla_\xi^2{\textbf{u}}^{(0)}&=p^{(0)}_Z{\textbf{e}}_2+\nabla_\xi p^{(1)}, \end{align}

with the boundary condition ${\textbf{u}}^{(0)}={\textbf{0}}$ on $\partial\omega_s(t)$ . As in [6], the problem is linear and so has solution of the form:

(3.9) \begin{align} {\textbf{u}}^{(0)}=- {\textbf{K}}(\xi,\eta,t)p^{(0)}_Z{\textbf{e}}_2, \qquad p^{(1)}=-{\boldsymbol\Pi}(\xi,\eta,t)\cdot {\textbf{e}}_2 p^{(0)}_Z+\pi_1(t,T),\end{align}

where ${\textbf{K}}$ is a tensor, ${\boldsymbol\Pi}$ a vector, and $\pi_1$ is a scalar. Substituting these forms into (3.8) we obtain a ‘cell problem’ for ${\textbf{K}}$ and ${\boldsymbol\Pi}$ , namely

(3.10a) \begin{align} \nabla_\xi\cdot{\textbf{K}}&={\textbf{0}}, \end{align}

(3.10b) \begin{align} \nabla_\xi^2{\textbf{K}}&=\nabla_\xi{\boldsymbol\Pi}-{\textbf{I}}, \end{align}

where ${\textbf{I}}$ is the identity matrix, along with the boundary condition ${\textbf{K}}={\textbf{0}}$ on $\partial\omega_s(t)$ , and the assumption of periodicity of ${\textbf{K}}$ and ${\boldsymbol\Pi}$ in $\xi$ , $\eta$ and t. We note that ${\textbf{K}}$ and ${\boldsymbol\Pi}$ depend on the short timescale t, since the position of the solid inclusion $\omega_s(t)$ depends on t.

We define the Darcy velocity to be the leading-order velocity, averaged over the cell, and over the micro-timescale:

(3.11) \begin{equation} {\textbf{q}}\;:\!=\; H_T\int_{t=0}^{1/H_T}\iint_{\omega_f(t)}{\textbf{u}}^{(0)}\,\mathrm{d}\xi\mathrm{d}\eta\,\mathrm{d}t,\end{equation}

where the average over the cell has been reduced to an integral over the fluid region since there is no flow in the solid region and the cell has unit area. Since we are assuming a one-dimensional macroscale flow with no X dependence, we note that ${\textbf{q}}$ is unidirectional, ${\textbf{q}}=q{\textbf{e}}_2$ .

While $p^{(0)}$ is independent of $\xi$ and $\eta$ , it is not necessarily independent of t, and so we define the time-averaged pressure to be

(3.12) \begin{equation} P(Z,T)\;:\!=\;H_T\int_{t=0}^{1/H_T}p^{(0)}\,\mathrm{d}t.\end{equation}

Substituting the definition (3.9) of ${\textbf{u}}^{(0)}$ into (3.11), we obtain Darcy’s law:

(3.13) \begin{equation} q=-H_T\int_{t=0}^{1/H_T}\iint_{\omega_f(t)}p^{(0)}_Z{\textbf{K}}_{22}\,\mathrm{d}\xi\mathrm{d}\eta\,\mathrm{d}t=-kP_Z,\end{equation}

using the definition of the average pressure P, and where we have defined the average permeability (for the 1D flow) to be

(3.14) \begin{equation} k=\iint_{\omega_f(t)}{\textbf{K}}_{22}\,\mathrm{d}\xi\mathrm{d}\eta.\end{equation}

The equations (3.10) are quasi-steady; the dependence of ${\textbf{K}}$ and ${\boldsymbol\Pi}$ on t is through the moving solid inclusions only. Due to the periodicity of ${\textbf{K}}$ in both time t and over the cell $\omega_c(t)$ , we see that the permeability, k, given by (3.14), must be independent of t.

We have found the form of Darcy’s law but are yet to derive an expression of overall conservation of mass. To do this we now look at the $O(\epsilon^2)$ problem, where we find that (3.2a) becomes

(3.15) \begin{equation} \nabla_\xi\cdot{\textbf{u}}^{(1)}+{\textbf{e}}_2\frac{\partial}{\partial Z}\cdot{\textbf{u}}^{(0)}=0.\end{equation}

We also require periodicity of ${\textbf{u}}^{(1)}$ and impose ${\textbf{u}}^{(1)}={\textbf{0}}$ on the solid surface $\partial \omega_s(t)$ . Integrating (3.15) over the fluid-occupied region of the cell, $\omega_f$ , using the divergence theorem, periodicity and boundary condition for ${\textbf{u}}^{(1)}$ , we find that

(3.16) \begin{align} 0&=\iint_{\omega_f(t)}{\textbf{e}}_2\frac{\partial}{\partial Z}\cdot{\textbf{u}}^{(0)}\,\mathrm{d}\xi\mathrm{d}\eta={\textbf{e}}_2\frac{\partial}{\partial Z}\cdot\left(\iint_{\omega_f(t)}{\textbf{u}}^{(0)}\,\mathrm{d}\xi\mathrm{d}\eta\right).\end{align}

Averaging over t (and using the fact that $H_T$ is independent of Z) gives the macroscale expression of conservation of mass:

(3.17) \begin{equation} \frac{\partial q}{\partial Z}=0.\end{equation}

This completes the analysis for the flow problem in the outer region: we have derived Darcy’s law (3.13) and the average conservation of mass equation (3.17).

3.2.2. Inner region

We now study the boundary layer near to the macroscale evaporating interface in order to derive the boundary conditions for the averaged flow problem in the outer region. We analyse the model (3.5a)–(3.5b) with the boundary conditions (3.5d), (3.5f) and (3.5h).

At leading order, we find, in $\omega_{{{\mathcal{H}}}^{(0)}}$ (which is bounded by the curve $z={{\mathcal{H}}}^{(0)}(x,t)$ , rather than $z={{\mathcal{H}}}(x,t)$ ),

(3.18a) \begin{align} \nabla_x\cdot\bar{{\textbf{u}}}^{(0)}&=0, \end{align}

(3.18b) \begin{align} \nabla_x\bar{p}^{(0)}&={\textbf{0}}, \end{align}

along with the boundary conditions, on $z={{\mathcal{H}}}^{(0)}(x,t)$ ,

(3.18c) \begin{align} \bar{{\textbf{u}}}^{(0)}\cdot{\textbf{m}}^{(0)}&=-(1-\delta)\frac{H_T+{{\mathcal{H}}}^{(0)}_t}{\sqrt{1+\left({{\mathcal{H}}}^{(0)}_x\right)^2}}, \end{align}

(3.18d) \begin{align} \bar{u}^{(0)}+\bar{v}^{(0)}{{\mathcal{H}}}^{(0)}_x&=0, \end{align}

where ${\textbf{m}}^{(0)}=\big(1/\sqrt{1+({{\mathcal{H}}}^{(0)})_x^2}\big)({{\mathcal{H}}}^{(0)}_x,-1)$ is the leading-order unit normal to the interface. On the solid surface $\partial\omega_s(t)$ ,

(3.18e) \begin{align} \bar{{\textbf{u}}}^{(0)}={\textbf{0}},\end{align}

along with periodicity of $\bar{p}^{(0)}$ and $\bar{{\textbf{u}}}^{(0)}$ in x and t. We recall that $\bar{{\textbf{u}}}^{(0)}$ and $p^{(0)}$ are allowed to be functions of x, z, t and T, but not of Z. However, we see from (3.18b) that $\bar{p}^{(0)}$ is independent of x and z.

Integrating equation (3.18a) for ${\textbf{u}}^{(0)}$ over the cell $\omega_{{{\mathcal{H}}}^{(0)}}(t)$ and using the divergence theorem, we find that

(3.19) \begin{align} 0=\int_{\partial\omega_{{{\mathcal{H}}}^{(0)}}}{\textbf{u}}^{(0)}\cdot{\textbf{n}}\,\mathrm{d}s.\end{align}

Using the periodicity in x and boundary conditions (3.18c) on $z={{\mathcal{H}}}^{(0)}(x,t)$ and (3.18e) on the solid, we see that

(3.20) \begin{align} \lim_{z_\Gamma\rightarrow\infty}\int_{\Gamma}\bar{{\textbf{u}}}^{(0)}\cdot{\textbf{n}}\,\mathrm{d}s=-\int_{z={{\mathcal{H}}}^{(0)}}\frac{(1-\delta)(H_T+{{\mathcal{H}}}_t^{(0)})}{\sqrt{1+\left({{\mathcal{H}}}_x^{(0)}\right)^2}}\,\mathrm{d}s ,\end{align}

which relates the flux of material through the evaporating interface to the far-field flux leaving the inner layer through the curve $\Gamma$ , as the height, $z_\Gamma$ of this curve is sent to infinity. We note (again by the divergence theorem, by constructing a closed region bounded by two choices $\Gamma_1$ and $\Gamma_2$ for $\Gamma$ , and by the curves $x=\pm1/2$ connecting them) that the left-hand side of (3.20) is independent of the choice of curve $\Gamma$ , with normal ${\textbf{n}}$ , so long as this connects $x=-{1/2}$ to $1/2$ , passing through the gas-filled pore space $\omega_{{{\mathcal{H}}}^{(0)}}$ . We average (3.20) over the period $[0,1/H_T]$ , say, of t to find that

(3.21) \begin{equation} \lim_{z_\Gamma\rightarrow\infty}H_T\int_{t=0}^{1/H_T}\int_{\Gamma}\bar{{\textbf{u}}}^{(0)}\cdot{\textbf{n}}\,\mathrm{d}s\,\mathrm{d}t =-(1-\delta)\mathcal{F}H_T,\end{equation}

where we have defined

(3.22) \begin{align} \mathcal{F} &=\int_{t=0}^{1/H_T} \int_{z={{\mathcal{H}}}^{(0)}}\frac{H_T+{{\mathcal{H}}}_t^{(0)}}{\sqrt{1+\left({{\mathcal{H}}}_x^{(0)}\right)^2}}\,\mathrm{d}s\,\mathrm{d}t =\int_{t=0}^{1/H_T}\int_{z={{\mathcal{H}}}^{(0)}}H_T+{{\mathcal{H}}}_t^{(0)}\,\mathrm{d}x\,\mathrm{d}t\nonumber\\ &=H_T\int_{t=0}^{1/H_T}\int_{z={{\mathcal{H}}}^{(0)}}\,\mathrm{d}x\,\mathrm{d}t,\end{align}

where we have used the fact that $H_T$ is independent of x, and that we may choose H and ${\mathcal{H}}$ such that the average interface velocity perturbation is zero. We interpret this parameter $\mathcal{F}$ as the total flux of material through the interface in a time period.

The value of $\mathcal{F}$ is not clear from considering only the dynamics of the gas. However, we may determine $\mathcal{F}$ by considering the conservation of mass of the liquid. Since we have assumed that the liquid does not flow, the total flux of liquid that passes through the evaporating interface over the time period $1/H_T$ (the time taken for the interface to move down through one repetition of the periodic pore-scale geometry) must equal the volume of liquid contained in that periodic cell. This is precisely the porosity, $\phi=|\omega_f|/|\omega_c|$ , of the porous medium, which is equivalently the cell volume occupied by fluid (and which we assume is constant).

More formally, we consider the liquid occupied region $\omega_l$ , illustrated in Figure 4, between $z=-1$ and the leading-order evaporating interface at $z={{\mathcal{H}}}^{(0)}$ . The rate of change of liquid mass contained in $\omega_l$ is given by:

(3.23) \begin{equation} R=\int_{\{z=-1\}\setminus\omega_s}\,\mathrm{d}x-\int_{z={{\mathcal{H}}}^{(0)}}\frac{H_T+{{\mathcal{H}}}_t^{(0)}}{\sqrt{1+({{\mathcal{H}}}_x^{(0)})^2}} \,\mathrm{d}s,\end{equation}

which is the rate at which liquid is gained as the liquid and solid move up through $z=-1$ (due to our moving coordinate system) minus the rate at which it is lost through the evaporating interface $z={{\mathcal{H}}}^{(0)}$ . The total liquid lost or gained over the time period $1/H_T$ must be zero as the process is periodic, and so

(3.24) \begin{align} 0&=\int_{t=0}^{1/H_T}R\,\mathrm{d}t=\int_{t=0}^{1/H_T}\left[\int_{\{z=-1\}\setminus\omega_s}\,\mathrm{d}x-\int_{z={{\mathcal{H}}}^{(0)}}\frac{H_T+{{\mathcal{H}}}_t^{(0)}}{\sqrt{1+({{\mathcal{H}}}_x^{(0)})^2}} \,\mathrm{d}s\right]\,\mathrm{d}t.\end{align}

The second integral term here is $\mathcal{F}$ by the definition (3.22), while the first is precisely the porosity, $\phi$ , of the porous medium. We therefore conclude that

(3.25) \begin{equation} \mathcal{F}=\phi=\iint_{\omega_f(t)}\,\mathrm{d}\xi\mathrm{d}\eta,\end{equation}

is the porosity of the porous medium.

Figure 4. Schematic showing the liquid region $\omega_l$ , consisting of the pore space between $z=-1$ and $z={{\mathcal{H}}}^{(0)}$ .

We now return to our expression (3.21) of the average velocity flux leaving the inner region. Since the left-hand side of (3.21) is independent of the choice of $\Gamma$ , we may match directly with the outer problem by letting the height, $z_{\Gamma}$ of $\Gamma$ become large, so that

(3.26) \begin{equation} \lim_{z_\Gamma\rightarrow\infty}\left( \int_{\Gamma}\bar{{\textbf{u}}}^{(0)}\cdot{\textbf{n}}\,\mathrm{d}s \right)=\left( \int_{\Gamma}{\textbf{u}}^{(0)}\cdot{\textbf{n}}\,\mathrm{d}s \right)\bigg|_{Z=0},\end{equation}

where on the right-hand side we consider $\Gamma$ to be a curve within the outer cell $\omega_f$ , joining $\xi=-1/2$ to $1/2$ . We choose $\Gamma$ to be the flat curve at $\eta=\eta_*$ , say, for some $\eta_*\in[-1/2,1/2]$ , wherever this lies within $\omega_f$ , and to follow the solid boundary $\partial\omega_s$ wherever the solid obstructs the curve $\eta=\eta_*$ . Since the integral along $\Gamma$ on the right-hand side of (3.26) is independent of the value of $\eta_*$ , this integral must equal its average over all possible $\eta_*\in[-1/2,1/2]$ . Thus, the left-hand side of (3.21) is

(3.27) \begin{align} \lim_{z_\Gamma\rightarrow\infty}\left(H_T\int_{t=0}^{1/H_T}\int_{\Gamma}\bar{{\textbf{u}}}^{(0)}\cdot{\textbf{n}}\,\mathrm{d}s\,\mathrm{d}t\right)&=\left(H_T\int_{t=0}^{1/H_T}\int_{\Gamma}{\textbf{u}}^{(0)}\cdot{\textbf{n}}\,\mathrm{d}s\,\mathrm{d}t\right)\bigg|_{Z=0}\nonumber\\ &=\left(H_T\int_{t=0}^{1/H_T}\iint_{\omega_f}{\textbf{u}}^{(0)}\cdot{\textbf{e}}_2\,\mathrm{d}\xi\mathrm{d}\eta\,\mathrm{d}t\right)\bigg|_{Z=0}=q|_{Z=0},\end{align}

which is precisely the (vertical) Darcy velocity. Substituting back into (3.21) (and using $\mathcal{F}=\phi$ ), we find the flux condition for the outer problem, namely

(3.28) \begin{equation} q=-(1-\delta)\phi H_T \qquad \text{ at }Z=0,\end{equation}

which provides the boundary condition for our outer, Darcy flow, problem.

We note that we do not need to find a macroscale boundary condition analogous to the no-tangential-slip condition (3.18d) that holds on the microscale, since the macroscale flow equations (3.13) and (3.17) only require one boundary condition for q on $Z=0$ (even in higher spatial dimensions), which is the normal flux condition (3.28). Since we do not need to explicitly compute the microscale velocity field $\bar{{\textbf{u}}}^{(0)}$ , the condition (3.18d) has not been used explicitly in the analysis, although it is formally needed for closure of the leading-order microscale system (3.18).

3.3. Advection–diffusion problem

Having homogenised the Stokes flow problem to Darcy’s law, with a flux boundary condition given in terms of the average interface velocity $H_T$ , we may now study the advection–diffusion problem for the vapour density $\rho$ . This analysis is similar to that in [16], although we now study an advection–diffusion process rather than simply a diffusion process. Indeed the system studied in [16] is the ‘no-flow’ limit, $\delta\rightarrow 1$ , of the present analysis. Our analysis in this section is also more general than [16], since we do not at this stage specify the chemistry interface condition, leaving this to be discussed in Section 4. Finally, we note that the matching process between the inner and intermediate regions given in [16] is not fully correct. This should instead be done as described in Section 3.3.3 below, resulting in a slightly different homogenised boundary condition to that derived in [16]. We provide a correction to [16] in Section 4.4.

As in [16], we will (in general) require a wider boundary layer of width $O(\epsilon^{1/2})$ at the interface $Z=0$ , in addition to the $O(\epsilon)$ inner region, in order to smooth out the oscillations in the vapour density due to the interface motion. This additional layer, which we refer to as the intermediate region, is necessary in order to match between the inner and outer regions. In this intermediate region, we change to the intermediate spatial variables, $\tilde{z}$ and $\tilde{x}$ , defined by:

(3.29) \begin{equation}Z=\epsilon^{1/2}\tilde{z}, \qquad X=\epsilon^{1/2}\tilde{x}.\end{equation}

In the intermediate region, we assume that the dependent variables are functions of both these intermediate variables $\tilde{z}$ and $\tilde{x}$ and of the microscale variables $\xi$ and $\eta$ . They also depend on both t and T but are independent of the spatial macroscale, Z. We use the notation $\tilde{\rho}$ for the vapour density in the intermediate region. With this change of variables, the advection–diffusion equation (2.9c) in the intermediate region becomes

(3.30a) \begin{align} \epsilon\delta&\left(\tilde{\rho}_t+\epsilon\tilde{\rho}_T-H_T\left(\tilde{\rho}_\eta+\epsilon^{1/2}\tilde{\rho}_{\tilde{z}}\right)\right)+\epsilon\tilde{{\textbf{u}}}\cdot\left(\nabla_\xi+\epsilon^{1/2}{\textbf{e}}_2\frac{\partial}{\partial \tilde{z}}\right)\tilde{\rho}\nonumber\\ &=\left(\nabla^2_\xi+\epsilon^{1/2}\left(\nabla_\xi\cdot{\textbf{e}}_2\frac{\partial}{\partial \tilde{z}}+{\textbf{e}}_2\frac{\partial}{\partial \tilde{z}}\cdot\nabla_\xi\right)+\epsilon\frac{\partial^2}{\partial \tilde{z}^2}\right)\tilde{\rho}.\end{align}

As in the outer region, $\tilde{\rho}$ is periodic with period 1, in $\xi,\eta$ , and with period $1/H_T$ , in t, and on the solid $\partial\omega_s(t)$ , the boundary condition (2.9d) becomes

(3.30b) \begin{align} \left(\nabla_\xi+\epsilon^{1/2}{\textbf{e}}_2\frac{\partial}{\partial \tilde{z}}\right)\tilde{\rho}\cdot{\textbf{n}}_s=0.\end{align}

We note that the intermediate problem (3.30) for $\tilde{\rho}$ depends on the velocity $\tilde{{\textbf{u}}}$ , but for the gas flow analysis in Section 3.2 we did not need to use information from this intermediate region. We will see that, for the analysis of the advection–diffusion problem in the intermediate region, we do not need to know anything about $\tilde{{\textbf{u}}}$ , except to assume it is well defined.

Since (3.30) involves powers of $\epsilon^{1/2}$ , we expand the dependent variables (in all of the inner, intermediate and outer regions) in powers of $\epsilon^{1/2}$ . We summarise the key facets of the analysis of the advection–diffusion problem, in the inner, intermediate and outer regions, as shown in Figure 5.

Figure 5. Summary of the results of the analysis of the advection–diffusion problem, at each order and in each region of the domain. The red arrows show the matchings leading to the effective interface conditions: the horizontal arrows at O(1) illustrate the leading-order vapour density matchings (for the chemistry interface condition that we will consider later in Section 4), and the diagonal arrows show the leading-order vapour flux matchings.

3.3.1. Outer region

We first analyse the vapour advection–diffusion problem (3.2c) with (3.2e), in the outer region of the domain. At leading order in $\epsilon$ , (3.2c) becomes

(3.31a) \begin{equation} \nabla_\xi^2\rho^{(0)}=0,\end{equation}

and the boundary condition (3.2e) becomes

(3.31b) \begin{equation} \nabla_\xi\rho^{(0)}\cdot{\textbf{n}}_s=0 \qquad \text{ on }\partial\omega_s(t),\end{equation}

which must be solved along with periodicity of $\rho^{(0)}$ in $\xi$ and $\eta$ . We multiply the equation (3.31a) by $\rho^{(0)}$ and integrate over the cell, using periodicity in $\xi$ and $\eta$ and the boundary condition to show that $\rho^{(0)}(t,Z,T)$ is independent of the microscale spatial variables $\xi$ and $\eta$ . At $O(\epsilon^{1/2})$ we obtain an identical problem for $\rho^{(1/2)}$ , and therefore find that $\rho^{(1/2)}(t,Z,T)$ is also independent of $\xi,\eta$ .

At $O(\epsilon)$ the equation (3.2c) reads

(3.32a) \begin{equation} \delta(\rho^{(0)}_t-H_T\rho^{(0)}_\eta)+{\textbf{u}}^{(0)}\cdot\nabla_\xi\rho^{(0)}=\nabla_\xi^2\rho^{(1)}+\left(\nabla_\xi\cdot{\textbf{e}}_2\frac{\partial}{\partial Z}+{\textbf{e}}_2\frac{\partial}{\partial Z}\cdot\nabla_\xi\right)\rho^{(0)},\end{equation}

subject to the periodicity of $\rho^{(1)}$ in $\xi$ , $\eta$ , and t and, from (3.2e), the boundary condition:

(3.32b) \begin{equation} (\nabla_\xi\rho^{(1)}+{\textbf{e}}_2\rho^{(0)}_Z)\cdot{\textbf{n}}_s=0 \qquad \text{ on }\partial\omega_s(t).\end{equation}

We note that the advection term at this order is zero since $\rho^{(0)}$ is independent of the micro-lengthscale and so, as in [16], we integrate (3.32a) over the cell to find that $\rho^{(0)}_t=0$ , and thus $\rho^{(0)}$ is a function of Z and T only.

The $O(\epsilon)$ problem therefore reduces to $\nabla_\xi^2\rho^{(1)}=0$ with the boundary condition (3.32b). We recall that $\rho^{(1)}$ is a function of $\xi,\eta,t,Z$ and T. Exploiting the linearity of the system, we write

(3.33) \begin{equation} \rho^{(1)}=W\rho^{(0)}_Z+F^{(0)},\end{equation}

where $F^{(0)}(t,Z,T)$ is independent of $\xi$ and $\eta$ , and $W(\xi,\eta,t)$ satisfies the cell problem:

(3.34a) \begin{align} \nabla_\xi^2W=0 \qquad &\text{ in }\omega_f(t),\end{align}

(3.34b) \begin{align} (\nabla_\xi W+{\textbf{e}}_2)\cdot{\textbf{n}}_s=0 \qquad &\text{ on }\partial\omega_s(t),\end{align}

(3.34c) \begin{align} W \text{ periodic in }\xi,\eta&\text{ and }t.\end{align}

We will not need any information from the order $O(\epsilon^{3/2})$ problem in the outer region, since the $O(\epsilon^2)$ problem (below) has no dependence on $\rho^{(3/2)}$ and neither will we need information from $O(\epsilon^{3/2})$ in order to derive the leading-order effective interface conditions. We therefore do not study the equations at this order. At $O(\epsilon^2)$ , we find that (3.2c) becomes

(3.35a) \begin{align} \delta\left(\rho^{(0)}_T-H_T\rho^{(0)}_Z+\rho^{(1)}_t-H_T\rho^{(1)}_\eta\right)&+\nabla_\xi\cdot\left(\rho^{(1)}{\textbf{u}}^{(0)}+\rho^{(0)}{\textbf{u}}^{(1)}-\left(\nabla_\xi\rho^{(2)}+{\textbf{e}}_2\rho^{(1)}_Z\right)\right)\nonumber\\[5pt] &+{\textbf{e}}_2\frac{\partial}{\partial Z}\cdot\left(\rho^{(0)}{\textbf{u}}^{(0)}-\left(\nabla_\xi\rho^{(1)}+{\textbf{e}}_2\rho^{(0)}_Z\right)\right)=0,\end{align}

and the boundary condition (3.2e) on $\partial\omega_s(t)$ becomes

(3.35b) \begin{equation} \left(\nabla_\xi\rho^{(2)}+{\textbf{e}}_2\rho_Z^{(1)}\right)\cdot{\textbf{n}}_s=0,\end{equation}

which must be solved along with the periodicity of $\rho^{(2)}$ . Integrating (3.35a) over the fluid-occupied region of the cell, $\omega_f(t)$ , using the divergence theorem, and applying the boundary and periodicity conditions for both $\rho^{(i)}$ and for ${\textbf{u}}^{(i)}$ , we find

(3.36) \begin{align} \delta\phi(\rho^{(0)}_T-H_T\rho^{(0)}_Z)&+\delta\frac{\mathrm{d}}{\mathrm{d}t}\left(\iint_{\omega_f(t)}\rho^{(1)}\,\mathrm{d}\xi\mathrm{d}\eta\right) \nonumber\\ & +\iint_{\omega_f(t)}{\textbf{e}}_2\frac{\partial}{\partial Z}\cdot\left(\rho^{(0)}{\textbf{u}}^{(0)}-(\nabla_\xi\rho^{(1)}+{\textbf{e}}_2\rho^{(0)}_Z)\right)\,\mathrm{d}\xi\mathrm{d}\eta=0,\end{align}

where we have used a transport theorem in t, which reads, for a general function F,

(3.37) \begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\left(\iint_{\omega_f(t)}F\,\mathrm{d}\xi\mathrm{d}\eta\right)&=\iint_{\omega_f(t)}F_t\,\mathrm{d}\xi\mathrm{d}\eta-\int_{\partial\omega_s(t)}FH_T{\textbf{e}}_2\cdot{\textbf{n}}_s\,\mathrm{d}s\nonumber\\ &=\iint_{\omega_f(t)}F_t\,\mathrm{d}\xi\mathrm{d}\eta-\iint_{\omega_f(t)}H_TF_\eta\,\mathrm{d}\xi\mathrm{d}\eta.\end{align}

This transport theorem encorporates the contributions to the total time derivative of the integral due to the time derivative of the integrand and the changing shape of the integration domain $\omega_f(t)$ . Since $\omega_f(t)$ changes by the solid inclusions moving up through the domain (in the ${\textbf{e}}_2$ direction), there is no ${\textbf{e}}_1$ term. The second line of (3.37) follows from the divergence theorem.

Using the form of $\rho^{(1)}$ given in (3.33) and the fact that $\rho^{(0)}$ is independent of $\xi$ and $\eta$ , equation (3.36) becomes

(3.38) \begin{align} \delta\phi(\rho^{(0)}_T-H_T\rho^{(0)}_Z)&+\delta\frac{\mathrm{d}}{\mathrm{d}t}\left(\iint_{\omega_f(t)}\rho^{(1)}\,\mathrm{d}\xi\mathrm{d}\eta\right) \nonumber\\ & ={\textbf{e}}_2\frac{\partial}{\partial Z}\cdot \left(\rho^{(0)}\iint_{\omega_f(t)}{\textbf{u}}^{(0)}\,\mathrm{d}\xi\mathrm{d}\eta-\rho^{(0)}_Z\iint_{\omega_f(t)}(\nabla_\xi W+{\textbf{e}}_2)\,\mathrm{d}\xi\mathrm{d}\eta\right).\end{align}

Averaging over t, the total-time-derivative term integrates to zero by the periodicity of $\rho^{(1)}$ in t, and so we arrive at a homogenised advection–diffusion equation describing the evolution of $\rho^{(0)}$ , namely

(3.39) \begin{align} \delta\phi&(\rho^{(0)}_T-H_T\rho^{(0)}_Z)+\frac{\partial}{\partial Z}\left(\rho^{(0)}q-\mathcal{D}\rho^{(0)}_Z\right)=0.\end{align}

We note that, although W and the cell structure $\omega_f$ vary with t (as the solid microstructure moves up through the domain), due to the periodicity of W over the cell and since the volume of the fluid-occupied region of the cell, $\omega_f$ , remains uniform, the effective diffusivity of the medium, taking into account the solid microstructure, is given by:

(3.40) \begin{equation} \mathcal{D}\;:\!=\;{\textbf{e}}_2\cdot\left(\iint_{\omega_f(t)}\left(\nabla_\xi W+{\textbf{e}}_2\right)\,\mathrm{d}\xi\mathrm{d}\eta\right)=\iint_{\omega_f(t)}\left(W_\eta+1\right)\,\mathrm{d}\xi\mathrm{d}\eta,\end{equation}

and is independent of t.

3.3.2. Intermediate region

We next consider the intermediate region, where the model is given by (3.30). The analysis of the intermediate region is closely linked with that of the inner region, which we will present in Section 3.3.3 below. We will require results from the inner region in order to complete the analysis in the intermediate region; these are highlighted in the discussion.

At O(1), equation (3.30a) becomes simply $\nabla_\xi^2\tilde{\rho}^{(0)}=0$ with (3.30b) giving homogeneous Neumann boundary conditions on the solid, along with the periodicity in t and at the cell boundary. Thus, we conclude, as in the outer region, that $\tilde{\rho}^{(0)}$ is independent of $\xi$ and $\eta$ .

At $O(\epsilon^{1/2})$ we find the problem (3.30a) in $\omega_f$ becomes

(3.41a) \begin{equation} \nabla_\xi^2\tilde{\rho}^{(1/2)}=0,\end{equation}

with periodicity and the boundary condition (3.30b) reading

(3.41b) \begin{equation} \left(\nabla_\xi\tilde{\rho}^{(1/2)}+{\textbf{e}}_2\tilde{\rho}^{(0)}_{\tilde{z}}\right)\cdot{\textbf{n}}_s=0 \qquad \text{ on }\partial\omega_s.\end{equation}

As in the outer problem (at $O(\epsilon)$ ), we can exploit the linearity of (3.41) to write $\tilde{\rho}^{(1/2)}$ , which we recall is allowed to be a function of $\xi,\eta,t,\tilde{z}$ and T, as:

(3.42) \begin{equation} \tilde{\rho}^{(1/2)}=W\tilde{\rho}^{(0)}_{\tilde{z}}+\tilde{f}^{(1/2)},\end{equation}

where $W(\xi,\eta,t)$ is the solution of the cell problem (3.34), and $\tilde{f}^{(1/2)}$ is independent of $\xi$ and $\eta$ .

The $O(\epsilon)$ version of (3.30a) is

(3.43a) \begin{equation} \delta \tilde{\rho}^{(0)}_t=\nabla_\xi^2\tilde{\rho}^{(1)}+\left(\nabla_\xi\cdot{\textbf{e}}_2\frac{\partial}{\partial \tilde{z}}+{\textbf{e}}_2\frac{\partial}{\partial \tilde{z}}\cdot\nabla_\xi\right)\tilde{\rho}^{(1/2)},\end{equation}

in $\omega_f$ , along with periodicity and the boundary condition (3.30b), which at $O(\epsilon)$ reads

(3.43b) \begin{equation} \left(\nabla_\xi\tilde{\rho}^{(1)}+{\textbf{e}}_2\tilde{\rho}^{(1/2)}_{\tilde{z}}\right)\cdot{\textbf{n}}_s=0 \qquad \text{ on }\partial\omega_s.\end{equation}

There are no advection terms in (3.43a) since $\nabla_\xi\tilde{\rho}^{(0)}=0$ . Similarly to the derivation of (3.39) in the outer problem (at $O(\epsilon^2)$ ), we now integrate (3.43a) over $\omega_f$ , using (i) the fact that $\tilde{\rho}^{(0)}$ is independent of $\xi$ and $\eta$ , (ii) the divergence theorem, (iii) the boundary condition (3.43b), (iv) the periodicity over the cell, and (v) the form (3.42) of $\tilde{\rho}^{(1/2)}$ . The outcome is that $\tilde{\rho}^{(0)}$ satisfies the diffusion equation:

(3.44) \begin{equation} \delta\phi\tilde{\rho}^{(0)}_t=\mathcal{D}\tilde{\rho}^{(0)}_{\tilde{z}\tilde{z}},\end{equation}

on the fast timescale t, and the intermediate lengthscale $\tilde{z}$ . Here $\mathcal{D}$ , the effective diffusivity, is given by (3.40). We note that the solution $\tilde{\rho}^{(0)}$ must be periodic in t and, by matching with the outer problem at leading order, we must have $\tilde{\rho}^{(0)}\rightarrow\rho^{(0)}|_{Z=0}$ as $\tilde{z}\rightarrow\infty$ . We may also match with the inner region at leading order. In Section 3.3.3, we will show that this matching gives the condition $\tilde{\rho}^{(0)}_{\tilde{z}}=0$ at $\tilde{z}=0$ . The solution of (3.44) is uniquely determined by these conditions (this uniqueness of solution may be shown by an energy method, in a similar way to in the Appendix of [16]), and is given by:

(3.45) \begin{equation} \tilde{\rho}^{(0)}=\rho^{(0)}|_{Z=0}.\end{equation}

We therefore conclude that $\tilde{\rho}^{(0)}$ is, in fact, independent both of $\tilde{z}$ and of t (since $\rho^{(0)}$ is independent of t). Furthermore, since $\tilde{\rho}^{(0)}_{\tilde{z}}$ is therefore everywhere 0, from (3.42) we see that $\tilde{\rho}^{(1/2)}$ is in fact independent of $\xi$ and $\eta$ .

With this knowledge we return to the $O(\epsilon)$ problem (3.43), and we see that (3.43a) reduces to simply $\nabla_\xi^2\tilde{\rho}^{(1)}=0$ . Arguing again by linearity, we see (3.43) therefore has solution:

(3.46) \begin{equation} \tilde{\rho}^{(1)}=W\tilde{\rho}^{(1/2)}_{\tilde{z}}+\tilde{f}^{1},\end{equation}

where $W(\xi,\eta,t)$ is the solution of the cell problem (3.34), and $\tilde{f}^{1}$ is independent of $\xi$ and $\eta$ .

The final order of interest is $O(\epsilon^{3/2})$ . Here, we find that (3.30a) becomes, in $\omega_f$ ,

(3.47a) \begin{align} \delta\left(\tilde{\rho}^{(1/2)}_t-H_T\tilde{\rho}^{(1/2)}_\eta\right)=\nabla_\xi\cdot\left(\nabla_\xi\tilde{\rho}^{(3/2)}+{\textbf{e}}_2\tilde{\rho}^{(1)}_{\tilde{z}}\right)+{\textbf{e}}_2\frac{\partial}{\partial\tilde{z}}\cdot\nabla_\xi\tilde{\rho}^{(1)}+\tilde{\rho}^{(1/2)}_{\tilde{z}\tilde{z}}\end{align}

(where there are no additional advection terms since both $\tilde{\rho}^{(0)}$ and $\tilde{\rho}^{(1/2)}$ are independent of the spatial microscale, and $\tilde{\rho}^{(0)}_{\tilde{z}}=0$ ). We also have the usual periodicity, while the boundary condition (3.30b) reads

(3.47b) \begin{equation} \left(\nabla_\xi\tilde{\rho}^{(3/2)}+{\textbf{e}}_2\tilde{\rho}^{(1)}_{\tilde{z}}\right)\cdot{\textbf{n}}_s=0 \qquad \text{ on }\partial\omega_s.\end{equation}

Again we argue as for the $O(\epsilon^2)$ problem in the outer region: we integrate the equation (3.47a) over the cell, using the transport theorem (3.37), the form (3.46) and apply both the boundary condition (3.47b) and periodicity, to find that

(3.48a) \begin{equation} \delta\phi\tilde{\rho}^{(1/2)}_t=\mathcal{D}\tilde{\rho}^{(1/2)}_{\tilde{z}\tilde{z}}.\end{equation}

We have therefore derived an effective equation for $\tilde{\rho}^{(1/2)}(t,\tilde{z},T)$ in the intermediate layer. The solution of (3.48a) must be periodic in t. Matching with the outer region, we see that we must have

(3.48b) \begin{equation} \tilde{\rho}^{(1/2)}_{\tilde{z}}\rightarrow\rho^{(0)}_Z|_{Z=0} \qquad\text{ as }\tilde{z}\rightarrow\infty.\end{equation}

We must also match with the inner region to close this problem for $\tilde{\rho}^{(1/2)}$ . As we will discuss in Section 3.3.3, there will be two matching conditions that relate $\tilde{\rho}^{(1/2)}$ and $\tilde{\rho}^{(1/2)}_{\tilde{z}}$ to the inner variables, namely

(3.48c) \begin{equation} \tilde{\rho}^{(1/2)}=\lim_{z\rightarrow\infty}\bar{\rho}^{(1/2)} \quad \text{ or } \quad \tilde{\rho}^{(1/2)}_{\tilde{z}}=\lim_{z\rightarrow\infty}\bar{\rho}^{(1)}_z \qquad \text{ at }\tilde{z}=0.\end{equation}

(These are the matching conditions (3.50c) and (3.50e) given below.) Regardless of which of these two matching conditions (3.48c), we view as the boundary condition closing (3.48), the time-periodic solution of (3.48) is unique. (The other condition (3.48c) then provides a far-field condition for the inner problem.) Unlike for $\tilde{\rho}^{(0)}$ , where the matching with the inner region required the solution of (3.34) to be independent of t and $\tilde{z}$ , we will see in Section 4 that the $\tilde{\rho}^{(1/2)}$ problem (3.48) has a solution with non-trivial dependence on t and $\tilde{z}$ , at least for certain chemistry interface conditions.

In summary, we have seen that the leading-order solution $\tilde{\rho}^{(0)}$ is constant through the intermediate region, but that the $O(\epsilon^{1/2})$ correction, $\tilde{\rho}^{(1/2)}$ , satisfies the effective diffusion equation (3.48a). By the matching (3.48b), we see that this problem (3.48) describes the evolution of the leading-order vapour flux through the intermediate region. It is precisely on this intermediate lengthscale $\tilde{z}$ that temporal oscillations in the leading-order vapour flux coming out of the inner region (off the evaporating interface) are smoothed out by diffusion.

3.3.3. Inner region

We now study the advection–diffusion problem for $\bar{\rho}$ and ${\mathcal{H}}$ in the inner region, namely (3.5c) with (3.5e) and (3.5g), with periodicity in x and t, and the appropriate matching conditions (given by (3.50), discussed below). The problem is closed using the chemistry condition (3.6) on $z={{\mathcal{H}}}^{(0)}$ , although we will not need to use this in the analysis of this section.

We will also require several results from the intermediate region, derived in Section 3.3.2, which we will incorporate by matching between the inner and intermediate regions. We note that this matching between the inner region, where there is a single variable z in the ${\textbf{e}}_2$ direction, and the intermediate region, where we have both $\eta$ and $\tilde{z}$ , is a subtle process. As z grows and we leave the inner region, we still expect fine scale oscillations in $\bar{\rho}$ due to the solid microstructure. In the limit $z\rightarrow\infty$ , we must therefore introduce multiple lengthscales by introducing the additional variable $\eta$ , and assuming that $\bar{\rho}$ depends on both $\eta$ and z independently in this limit. Thus, our derivatives in the limit $z\rightarrow \infty$ must be replaced with the multiple scales ansatz:

(3.49) \begin{equation} \nabla_x\rightarrow \nabla_\xi+{\textbf{e}}_2\frac{\partial}{\partial z} \qquad \text{ as }z\rightarrow\infty.\end{equation}

We derive the Van Dyke matching conditions between the inner and intermediate regions, following [4], and the conditions we will need for our analysis are as follows:

(3.50a) \begin{align} \lim_{z\rightarrow\infty}(\bar{\rho}^{(0)})&=\tilde{\rho}^{(0)}|_{\tilde{z}=0},\end{align}

(3.50b) \begin{align} \lim_{z\rightarrow\infty}(\bar{\rho}^{(1/2)}_z)&=\tilde{\rho}^{(0)}_{\tilde{z}}|_{\tilde{z}=0},\end{align}

(3.50c) \begin{align} \lim_{z\rightarrow\infty}(\bar{\rho}^{(1/2)}-z\tilde{\rho}^{(0)}_{\tilde{z}}|_{\tilde{z}=0})&=\tilde{\rho}^{(1/2)}|_{\tilde{z}=0},\end{align}

(3.50d) \begin{align} \lim_{z\rightarrow\infty}(\bar{\rho}^{(1)}_{zz})&=\tilde{\rho}^{(0)}_{\tilde{z}\tilde{z}}|_{\tilde{z}=0},\end{align}

(3.50e) \begin{align} \lim_{z\rightarrow\infty}\left(\bar{\rho}^{(1)}_z-\frac{z}{2}\tilde{\rho}^{(0)}_{\tilde{z}\tilde{z}}|_{\tilde{z}=0}\right)&=\tilde{\rho}^{(1/2)}_{\tilde{z}}|_{\tilde{z}=0},\end{align}

(3.50f) \begin{align} \lim_{z\rightarrow\infty}\left(\bar{\rho}^{(1)}-z\tilde{\rho}^{(1/2)}_{\tilde{z}}|_{\tilde{z}=0}-\frac{z^2}{2}\tilde{\rho}^{(0)}_{\tilde{z}\tilde{z}}|_{\tilde{z}=0}\right)&=\tilde{\rho}^{(1)}|_{\tilde{z}=0}.\end{align}

We note that the factors of z and $z^2$ in these matching conditions (3.50) are the growing lengthscale z, not $\eta$ , and we emphasise that $ \bar{\rho}_z^i=\bar{\rho}_z^i(\eta,z)$ in this large-z limit.

At leading order, the inner advection–diffusion problem (3.5c), in $\omega_{{{\mathcal{H}}}^{(0)}}$ , becomes

(3.51a) \begin{equation} \nabla_x^2\bar{\rho}^{(0)}=0,\end{equation}

which must be solved with periodicity in x and t, while the boundary conditions (3.5e) and (3.5g) become

(3.51b) \begin{align} \nabla_x\bar{\rho}^{(0)}\cdot{\textbf{n}}_s=0 \qquad &\text{ on }\partial\omega_s, \end{align}

(3.51c) \begin{align} \nabla_x\bar{\rho}^{(0)}\cdot{\textbf{m}}^{(0)}=0, \qquad &\text{ on }z={{\mathcal{H}}}^{(0)}. \end{align}

Multiplying (3.51a) by $\bar{\rho}^{(0)}$ , integrating over the cell (using divergence theorem), and using the no-flux boundary conditions (3.51b)–(3.51c) and periodicity, we find that

(3.52) \begin{equation} \iint_{\omega_{{{\mathcal{H}}}^{(0)}}}|\nabla_x\bar{\rho}^{(0)}|^2\,\mathrm{d}x\mathrm{d}z=\lim_{z\rightarrow\infty}\int_{\Gamma}\bar{\rho}^{(0)}\nabla_x\bar{\rho}^{(0)}\cdot{\textbf{n}}\,\mathrm{d}s,\end{equation}

where, as in Section 3.2.2, $\Gamma$ is any curve through $\omega_{{{\mathcal{H}}}^{(0)}}$ from $x=-1/2$ to $x=1/2$ , at height z, with normal ${\textbf{n}}$ . As we take the limit $z\rightarrow\infty$ , we introduce the multiple scales via (3.49) and use the matching (3.50a) so that (3.52) becomes

(3.53) \begin{equation} \iint_{\omega_{{{\mathcal{H}}}^{(0)}}}|\nabla_x\bar{\rho}^{(0)}|^2\,\mathrm{d}x\mathrm{d}z=\int_{\Gamma}\tilde{\rho}^{(0)}(\nabla_\xi\tilde{\rho}^{(0)}+{\textbf{e}}_2\tilde{\rho}^{(0)}_z)\cdot{\textbf{n}}\,\mathrm{d}s.\end{equation}

In the intermediate layer, we found that the leading-order density $\tilde{\rho}^{(0)}$ is a function only of t and T. Thus, the right-hand side of (3.53) is zero, and we conclude from (3.53) that $\bar{\rho}^{(0)}$ is independent of x and z. Furthermore, since in the intermediate layer $\tilde{\rho}^{(0)}$ is independent of $\tilde{z}$ , we may match directly through the intermediate layer at this order, to see that

(3.54) \begin{equation} \bar{\rho}^{(0)}=\tilde{\rho}^{(0)}|_{\tilde{z}=0}=\lim_{\tilde{z}\rightarrow\infty}\tilde{\rho}^{(0)}=\rho^{(0)}|_{Z=0}.\end{equation}

Thus, since $\rho^{(0)}$ in the outer region is independent of t, we also know that $\bar{\rho}^{(0)}$ is independent of t.

At the next order, $O(\epsilon^{1/2})$ , of the inner problem (3.5c) we find the same problem for $\bar{\rho}^{(1/2)}$ as we had for $\bar{\rho}^{(0)}$ at leading order, namely, in $\omega_{{{\mathcal{H}}}^{(0)}}$ ,

(3.55a) \begin{equation} \nabla_x^2\bar{\rho}^{(1/2)}=0,\end{equation}

which we solve with periodicity in x and t, and the boundary conditions (3.5e) and (3.5g) become

(3.55b) \begin{align} \nabla_x\bar{\rho}^{(1/2)}\cdot{\textbf{n}}_s=0, \qquad &\text{ on }\partial\omega_s, \end{align}

(3.55c) \begin{align} \nabla_x\bar{\rho}^{(1/2)}\cdot{\textbf{m}}^{(0)}=0, \qquad &\text{ on }z={{\mathcal{H}}}^{(0)}.\end{align}

We might expect there to be additional terms in the $O(\epsilon^{1/2})$ problem, from expanding ${{\mathcal{H}}}={{\mathcal{H}}}^{(0)}+\epsilon^{1/2}{{\mathcal{H}}}^{(1/2)}+...$ and from evaluating the boundary conditions at ${{\mathcal{H}}}^{(0)}$ rather than ${{\mathcal{H}}}$ . However, these terms are zero, since $\bar{\rho}^{(0)}$ is independent of the microscale. Integrating (3.55a) over the domain $\omega_{{{\mathcal{H}}}^{(0)}}$ , and using the divergence theorem, periodicity and the boundary conditions (3.55b)–(3.55c), we find that

(3.56) \begin{equation} 0=\lim_{z\rightarrow\infty}\int_{\Gamma}\nabla_x\bar{\rho}^{(1/2)}\cdot{\textbf{n}}\,\mathrm{d}s=\lim_{z\rightarrow\infty}\int_{\Gamma}\left(\nabla_\xi+{\textbf{e}}_2\frac{\partial}{\partial z}\right)\bar{\rho}^{(1/2)}\cdot{\textbf{n}}\,\mathrm{d}s,\end{equation}

where we have made use of the multiple scales ansatz (3.49). Using the matching condition (3.50c), we see that

(3.57) \begin{equation} \bar{\rho}^{(1/2)}\sim z\tilde{\rho}^{(0)}_{\tilde{z}}|_{\tilde{z}=0}+\tilde{\rho}^{(1/2)}|_{\tilde{z}=0}=\tilde{\rho}^{(0)}_{\tilde{z}}|_{\tilde{z}=0}\left(z+W(\xi,\eta,t)\right) \qquad \text{ as }z\rightarrow\infty,\end{equation}

where for the equality we have used the form (3.42) of $\tilde{\rho}^{(1/2)}$ . (We note that the form of $\bar{\rho}^{(1/2)}$ given by (3.57) is compatible with the matching condition (3.50b) precisely since the z derivative in (3.50b) is only with respect to the growing z variable, and not $\eta$ .) Substituting (3.57) into (3.56), and since $\tilde{\rho}^{(0)}$ is independent of $\xi$ and $\eta$ , we have

(3.58) \begin{equation} 0=\tilde{\rho}^{(0)}_{\tilde{z}}|_{\tilde{z}=0}\int_{\Gamma}\left(\nabla_\xi W+{\textbf{e}}_2\right)\cdot{\textbf{n}}\,\mathrm{d}s.\end{equation}

The integral on the right-hand side is independent of the choice of curve $\Gamma$ by definition of W, and is non-zero, and thus we conclude that $\tilde{\rho}^{(0)}_{\tilde{z}}|_{\tilde{z}=0}=0$ . (This fact was used in the analysis of the intermediate region, to show that $\tilde{\rho}^{(0)}$ is in fact everywhere independent of $\tilde{z}$ , and that $\tilde{\rho}^{(1/2)}$ is independent of $\xi$ and $\eta$ .)

Returning to the $O(\epsilon^{1/2})$ inner problem (3.55) once more, we multiply (3.55a) through by $\bar{\rho}^{(1/2)}$ and integrate over $\omega_{{{\mathcal{H}}}^{(0)}}$ (much as for the O(1) inner problem), to find

(3.59) \begin{align} \iint_{\omega_{{{\mathcal{H}}}^{(0)}}}|\nabla_x\bar{\rho}^{(1/2)}|^2\,\mathrm{d}x\mathrm{d}z&=\lim_{z\rightarrow\infty}\int_{\Gamma}\bar{\rho}^{(1/2)}\nabla_x\bar{\rho}^{(1/2)}\cdot{\textbf{n}}\,\mathrm{d}s\nonumber\\ &=\lim_{z\rightarrow\infty}\int_{\Gamma}\bar{\rho}^{(1/2)}\left(\nabla_\xi+{\textbf{e}}_2\frac{\partial}{\partial z}\right)\bar{\rho}^{(1/2)}\cdot{\textbf{n}}\,\mathrm{d}s\nonumber\\ &=\lim_{z\rightarrow\infty}\int_{\Gamma}(\tilde{\rho}^{0}_{\tilde{z}}|_{\tilde{z}=0})^2\left(W+z\right)\left(\nabla_\xi W+{\textbf{e}}_2\right)\cdot{\textbf{n}}\,\mathrm{d}s,\end{align}

again using (3.49), and (3.57). Since $\tilde{\rho}^{(0)}_{\tilde{z}}|_{\tilde{z}=0}=0$ , the right-hand side is 0 and thus $\bar{\rho}^{(1/2)}$ is also independent of $\xi$ and $\eta$ .

The final order of interest in the inner region is $O(\epsilon)$ . At this order, (3.5c) becomes, in $\omega_{{{\mathcal{H}}}^{(0)}}$ ,

(3.60a) \begin{equation}\nabla_x^2\bar{\rho}^{(1)}=0,\end{equation}

which must be solved with periodicity in x and t and the boundary conditions (3.5e) and (3.5g), which become

(3.60b) \begin{align} \nabla_x\bar{\rho}^{(1)}\cdot{\textbf{n}}_s=0, \qquad &\text{ on }\partial\omega_s, \end{align}

(3.60c) \begin{align} \nabla_x\bar{\rho}^{(1)}\cdot{\textbf{m}}^{(0)}=(1-\nu\bar{\rho}^{(0)})\frac{H_T+{{\mathcal{H}}}^{(0)}_t}{\sqrt{1+\left({{\mathcal{H}}}^{(0)}_x\right)^2}}, \qquad &\text{ on }z={{\mathcal{H}}}^{(0)}.\end{align}

Again we have used the fact that $\bar{\rho}^{(0)}$ is independent of x, z and t, as well as the fact that $\bar{{\textbf{u}}}^{(0)}$ is divergence free in the microscale variables (as shown in Section 3.2.2). Since $\bar{\rho}^{(1/2)}$ is independent of the spatial microscale, we have no correction terms due to using the domain $\omega_{{{\mathcal{H}}}^{(0)}}$ rather than $\omega_{{\mathcal{H}}}$ . The problem is closed using the matching condition (3.50f) as $z\rightarrow\infty$ , which, since we know $\tilde{\rho}^{(1)}=W\tilde{\rho}^{(1/2)}+\tilde{f}^{(1)}$ in the intermediate region, is given by:

(3.60d) \begin{equation} \lim_{z\rightarrow\infty}(\bar{\rho}^{(1)}-z\tilde{\rho}^{(1/2)}_{\tilde{z}}|_{\tilde{z}=0})=\tilde{\rho}^{(1/2)}_{\tilde{z}}|_{\tilde{z}=0}W(\xi,\eta,t)+\tilde{f}^{(1)}_{\tilde{z}}|_{\tilde{z}=0}.\end{equation}

The problem (3.60) along with the appropriate chemistry boundary condition forms a closed system for $\bar{\rho}^{(1)}$ and ${{\mathcal{H}}}^{(0)}$ , given in terms of $H_T$ and $\tilde{\rho}^{(1/2)}_{\tilde{z}}|_{\tilde{z}=0}$ .

However, without computing $\bar{\rho}^{(1)}$ and ${{\mathcal{H}}}^{(0)}$ we may derive the effective vapour–flux interface condition for the macroscale problem. Integrating (3.60a) over $\omega_{{{\mathcal{H}}}^{(0)}}$ , using the divergence theorem, and the periodicity and boundary conditions (3.60b)–(3.60c), we find that the total diffusive flux through the evaporating interface equals the flux out of the inner cell as $z\rightarrow\infty$ , so that

(3.61) \begin{align} \int_{z={{\mathcal{H}}}^{(0)}}(1-\nu\bar{\rho}^{(0)})\frac{H_T+{{\mathcal{H}}}^{(0)}_t}{\sqrt{1+({{\mathcal{H}}}^{(0)}_x)^2}}\,\mathrm{d}s&=\lim_{z\rightarrow\infty}\int_{\Gamma}\nabla_x\bar{\rho}^{(1)}\cdot{\textbf{n}}\,\mathrm{d}s\nonumber\\ &=\lim_{z\rightarrow\infty}\int_{\Gamma}\left(\nabla_\xi+{\textbf{e}}_2\frac{\partial}{\partial z}\right)\bar{\rho}^{(1)}\cdot{\textbf{n}}\,\mathrm{d}s\nonumber\\ &=\tilde{\rho}^{(1/2)}_{\tilde{z}}|_{\tilde{z}=0}\int_{\Gamma}\left(\nabla_\xi W+{\textbf{e}}_2\right)\cdot{\textbf{n}}\,\mathrm{d}s,\end{align}

where for the final equality we have used the matching condition (3.60d), and the fact that $\tilde{\rho}^{(1/2)}$ is independent of $\xi$ and $\eta$ .

Since W satisfies the cell problem (3.34), the function ${\textbf{F}}\;:\!=\;\nabla_\xi W+{\textbf{e}}_2$ is divergence free on the cell $\omega_f$ , and ${\textbf{F}}\cdot{\textbf{n}}_s=0$ on the solid $\partial\omega_s$ . The integral over the curve $\Gamma$ in (3.61) is therefore independent of the choice of $\Gamma$ (by application of the divergence theorem). Similarly to the analysis of Section 3.2.2, by choosing $\Gamma$ to be the flat curve $\eta=\eta^*$ where this lies within $\omega_f$ , and on the surface $\partial\omega_s$ wherever the solid obstructs the line $\eta=\eta^*$ , we see that the integral in (3.61),

(3.62) \begin{equation} \int_{\Gamma}\left(\nabla_\xi W+{\textbf{e}}_2\right)\cdot{\textbf{n}}\,\mathrm{d}s=\int_{\{\eta=\eta^*\}\cap\omega_f}W_\eta+1\,\mathrm{d}\xi,\end{equation}

is independent of the value of $\eta^*$ . Thus, this integral is equal to its average over $\eta^*$ ,

(3.63) \begin{equation} \int_{\Gamma}\left(\nabla_\xi W+{\textbf{e}}_2\right)\cdot{\textbf{n}}\,\mathrm{d}s=\iint_{\omega_f}W_\eta+1\,\mathrm{d}\xi\mathrm{d}\eta=\mathcal{D},\end{equation}

and so is precisely the effective diffusivity $\mathcal{D}$ . Substituting this into (3.61), we find

(3.64) \begin{equation} (1-\nu\bar{\rho}^{(0)})\int_{z={{\mathcal{H}}}^{(0)}}\frac{H_T+{{\mathcal{H}}}^{(0)}_t}{\sqrt{1+({{\mathcal{H}}}^{(0)}_x)^2}}\,\mathrm{d}s=\mathcal{D}\tilde{\rho}^{(1/2)}_{\tilde{z}}|_{\tilde{z}=0},\end{equation}

where we have used the fact that $\bar{\rho}^{(0)}$ is independent of x and z to rearrange the left-hand side. This equation (3.64) relates the instantaneous leading-order diffusive flux of vapour entering the intermediate layer, to the flux of vapour off the interface. We note that the integral on the left of (3.64) is the instantaneous total interface velocity, and that this depends on t, since the shape and velocity of the microscale interface, $z={{\mathcal{H}}}^{(0)}$ , in general vary with t. This is the reason why we need an intermediate layer: we could not match directly between the inner and outer regions because the leading-order outer density flux $\rho^{(0)}_Z$ is independent of the micro-timescale t. The role of the intermediate region is that it allows these temporal oscillations to be smoothed out.

We now average (3.64) in t, using the fact that $\bar{\rho}^{(0)}$ is independent of t and the definition (3.22) of $\mathcal{F}$ (which from our discussion in Section 3.2.2 we know to be given by $\mathcal{F}=\phi$ ), to simplify the left-hand side, to give

(3.65) \begin{equation} (1-\nu\bar{\rho}^{(0)})\phi H_T=\mathcal{D}H_T\int_{t=0}^{1/H_T}\tilde{\rho}^{(1/2)}_{\tilde{z}}|_{\tilde{z}=0}\,\mathrm{d}t.\end{equation}

We note that by integrating (3.48a) over the period of t, and using the periodicity of $\tilde{\rho}^{(1/2)}$ , we must have

(3.66) \begin{equation} 0=H_T\int_{t=0}^{1/H_T}\mathcal{D}\tilde{\rho}^{(1/2)}_{\tilde{z}\tilde{z}}\,\mathrm{d}t=\left(H_T\int_{t=0}^{1/H_T}\mathcal{D}\tilde{\rho}^{(1/2)}_{\tilde{z}}\,\mathrm{d}t\right)_{\tilde{z}},\end{equation}

so that the time-averaged diffusive flux of vapour through the intermediate region is independent of $\tilde{z}$ . In particular, we may therefore replace $\tilde{z}=0$ with the limit $\tilde{z}\rightarrow\infty$ in (3.65). We then apply the matching condition:

(3.67) \begin{equation} \lim_{\tilde{z}\rightarrow\infty}\tilde{\rho}^{(1/2)}_{\tilde{z}}=\rho^{(0)}_Z|_{Z=0},\end{equation}

and the fact that $\rho^{(0)}_Z$ is independent of t to find that

(3.68) \begin{align} (1-\nu\bar{\rho}^{(0)})\phi H_T =\mathcal{D}\rho^{(0)}_Z|_{Z=0}.\end{align}

Finally, using (3.54), we see that (3.68) is an effective boundary condition for the macroscale leading-order vapour density $\rho^{(0)}$ , namely

(3.69) \begin{equation} (1-\nu\rho^{(0)})\phi H_T=\mathcal{D}\rho^{(0)}_Z \qquad \text{ at }Z=0.\end{equation}

3.4. Summary so far

In summary, we have so far derived the following effective model in $Z>0$ ,

(3.70a) \begin{align} \frac{\partial q}{\partial Z}&=0,\end{align}

(3.70b) \begin{align} q&=-kP_Z,\end{align}

(3.70c) \begin{align} \delta\phi(\rho^{(0)}_T-H_T\rho^{(0)}_Z)&=\frac{\partial}{\partial Z}\left(\mathcal{D}\rho^{(0)}_Z-\nu\rho^{(0)}q\right),\end{align}

with the interface conditions at $Z=0$ :

(3.70d) \begin{align} q&=-(1-\delta)\phi H_T, \end{align}

(3.70e) \begin{align} \mathcal{D}\rho^{(0)}_Z&=(1-\nu\rho^{(0)})\phi H_T, \end{align}

where the leading-order macroscale-dependent variables, q, P and $\rho^{(0)}$ , are functions of the macroscale-independent variables, Z and T, only, and the interface position H is a function of T only. The effective permeability k and diffusivity $\mathcal{D}$ , given by (3.14) and (3.40), are found by solution of the cell problems (3.10) and (3.34), respectively, for the specific pore space geometry.

We note that, by combining (3.70d) and (3.70e), we find that

(3.71) \begin{equation} \left(\rho^{(0)}\nu(q-\delta\phi H_T)-\mathcal{D}\rho^{(0)}_Z\right)|_{Z=0}=-\phi H_T,\end{equation}

and thus the total flux of vapour off the interface is equal to the flux of liquid into the interface. Therefore, by assuming conservation of liquid mass on the microscale, we have ensured that our macroscale model also conserves liquid mass.

We also note that the shape and velocity variations of the microscale interface do not affect these flux boundary conditions (3.70d) and (3.70e) at all: we might have written these down with no thought about the microscale behaviour. The chemistry interface condition will be analysed in the following section and will provide one further effective interface condition at the macroscale evaporating interface $Z=0$ . We will find that, for certain ranges of $\alpha$ , our detailed multiscale analysis is necessary in order to derive this final interface condition and that here the microscale interface geometry does play a role.

4. Chemistry interface condition

We now close the homogenised model (3.70), derived in Section 3, by deriving the effective chemistry interface condition on the evaporating interface. We recall that the microscale chemistry condition, given by (3.6), takes the form:

(4.1) \begin{equation}1-\bar{\rho}=-\alpha \frac{H_T+{{\mathcal{H}}}_t}{\sqrt{1+{{\mathcal{H}}}_x^2}} \qquad \text{ on }z={{\mathcal{H}}}(x,t).\end{equation}

In this section, we will look at several different distinguished limits, namely $\alpha=O(\epsilon)$ , $\alpha=O(\epsilon^{1/2})$ and $\alpha=O(1)$ . As well as leading to different forms of the effective interface condition, we will also find quite dramatically different microscale behaviours for the different ranges of $\alpha$ . The main findings in each of these cases are summarised in Table 1.

Table 1. Summary and comparison of the main results of the chemistry interface condition, for each range of $\alpha$

4.1. Case $\alpha\ll1$

We begin with the case of small $\alpha\ll1$ so that the vapour density equilibrates to the saturation density at much faster rate than the rate of interface motion. For $\alpha\ll1$ , to leading order, the chemistry interface condition (4.1) is

(4.2) \begin{equation} \bar{\rho}^{(0)}=1 \qquad \text{ on }z={{\mathcal{H}}}^{(0)},\end{equation}

so that $\bar{\rho}^{(0)}\equiv 1$ throughout the inner region, as it is independent of x, z and t. By matching at leading order (through the intermediate region as before according to (3.54)), we derive the effective interface condition:

(4.3) \begin{equation} \rho^{(0)}=1 \qquad \text{ on }Z=0,\end{equation}

which closes the homogenised model (3.70) in this case. We note that in this case, the chemistry condition becomes a Dirichlet condition on $\rho^{(0)}$ , while the diffusive flux condition (3.70e) determines the interface velocity $H_T$ .

For any $\alpha\ll1$ the Dirichlet boundary condition (4.3) is the correct macroscale interface condition. However, we next explore the different microscale behaviours of the distinguished limits $\alpha=O(\epsilon)$ (which contains the behaviour for $\alpha\ll\epsilon$ ) and $\alpha=O(\epsilon^{1/2})$ .

4.1.1. Case $\alpha=O(\epsilon)$ or smaller

If $\alpha\ll\epsilon^{1/2}$ then at $O(\epsilon^{1/2})$ , the chemistry interface condition (4.1) becomes

(4.4) \begin{equation} \bar{\rho}^{(1/2)}=0 \qquad \text{ on }z={{\mathcal{H}}}^{(0)}.\end{equation}

Since we know that $\bar{\rho}^{(1/2)}$ is independent of the spatial variables x and z, we conclude that $\bar{\rho}^{(1/2)}\equiv 0$ is always zero throughout the inner region. This is surprising: we have no $O(\epsilon^{1/2})$ correction term in the inner region.

Furthermore, we can close the problem (3.48a) for $\tilde{\rho}^{(1/2)}$ in the intermediate region with the Dirichlet boundary condition (3.48c, left), which is

(4.5) \begin{equation} \tilde{\rho}^{(1/2)}|_{\tilde{z}=0}=\lim_{z\rightarrow\infty}\bar{\rho}^{(1/2)}=0.\end{equation}

The unique time-periodic solution of (3.48) is then

(4.6) \begin{equation} \tilde{\rho}^{(1/2)}=\tilde{z}\rho^{(0)}_Z|_{Z=0},\end{equation}

which is, in fact, independent of t. The intermediate region in this case is not performing its expected function of smoothing temporal oscillations in the flux off the interface. Indeed, the flux of vapour through the intermediate region, $\rho^{(0)}_Z|_{Z=0}$ , is constant: the intermediate region does nothing at this order. Since (4.6) holds, the far-field condition (3.60d) for the inner region flux problem (3.60) becomes simply

(4.7) \begin{equation} \bar{\rho}^{(1)}\sim\rho^{(0)}_Z|_{Z=0}(z+W) \qquad \text{ as }z\rightarrow\infty.\end{equation}

Thus, in the case when $\alpha\ll\epsilon^{1/2}$ , we could have greatly simplified the analysis of Section 3.3, as we do not need an intermediate layer or an expansion in terms of $O(\epsilon^{1/2})$ at all.

Substituting (4.6) into (3.64) (and using the fact that $\bar{\rho}^{(0)}=1$ ), we find that the instantaneous total flux of vapour coming off the interface on the microscale, namely

(4.8) \begin{equation} (1-\nu)\int_{z={{\mathcal{H}}}^{(0)}}\frac{H_T+{{\mathcal{H}}}^{(0)}_t}{\sqrt{1+\left({{\mathcal{H}}}^{(0)}_x\right)^2}}\,\mathrm{d}s=\mathcal{D}\rho^{(0)}_Z|_{Z=0},\end{equation}

is independent of t (since the right-hand side is independent of t). The shape and/or speed of the microscale interface $z={{\mathcal{H}}}^{(0)}(x,t)$ must therefore vary in time as the solid microstructure moves up through the cell in order to maintain this constant total flux off the interface. For instance, at times when the fluid-occupied fraction of the curve $z=0$ is smaller than average, the interface $z={{\mathcal{H}}}^{(0)}(x,t)$ must move faster than average in order to maintain the same total evaporation flux. The shape and velocity of the microscale interface $z={{\mathcal{H}}}^{(0)}$ may be found, along with $\bar{\rho}^{(1)}$ , by solving the problem (3.60). However, we note that there is no need to compute this solution in order to understand the macroscale evaporation process, since the homogenised model is not affected in any way by the details of the microscale interface motion. Indeed, the macroscale, outer process informs the inner: the microscale interface shape and speed adapt to maintain the t-independent vapour flux drawn off the interface by the macroscale density gradient.

4.1.2. Case $\alpha=O(\epsilon^{1/2})$

We now suppose that $\alpha=O(\epsilon^{1/2})$ and write $\alpha=\beta\epsilon^{1/2}$ with $\beta=O(1)$ . While $\bar{\rho}^{(0)}$ satisfies the Dirichlet condition (4.2) at leading order, at $O(\epsilon^{1/2})$ (4.1) becomes

(4.9) \begin{equation} \beta\frac{H_T+{{\mathcal{H}}}^{(0)}_t}{\sqrt{1+\left({{\mathcal{H}}}^{(0)}_x\right)^2}}=\bar{\rho}^{(1/2)} \qquad \text{ on }z={{\mathcal{H}}}^{(0)}.\end{equation}

We recall that, like $\bar{\rho}^{(0)}$ , $\bar{\rho}^{(1/2)}$ is independent of x and z (although unlike $\bar{\rho}^{(0)}$ it may depend on t).

The value of $\bar{\rho}^{(1/2)}(t)=\tilde{\rho}^{(1/2)}|_{\tilde{z}=0}$ is found by solving the diffusion equation for $\tilde{\rho}^{(1/2)}$ in the intermediate region (3.48) (with the boundary condition (3.48c, right)). The two problems (3.60) and (3.48) are therefore fully coupled by both (4.9) and the matching condition (3.60d) and must be solved simultaneously in order to understand the microscale interface motion. We note that we can simplify this problem by deriving a new interface condition for $\tilde{\rho}^{(1/2)}$ that holds at $\tilde{z}=0$ , rather than matching with $\bar{\rho}^{(1)}$ in the inner region. Specifically, we note from (3.64) that the instantaneous flux of vapour out of the inner region is

(4.10) \begin{equation} \mathcal{D}\tilde{\rho}^{(1/2)}_{\tilde{z}}|_{\tilde{z}=0}=(1-\nu)\int_{z={{\mathcal{H}}}^{(0)}}\frac{H_T+{{\mathcal{H}}}^{(0)}_t}{\sqrt{1+({{\mathcal{H}}}^{(0)}_x)^2}}\,\mathrm{d}s=\frac{(1-\nu)}{\beta}L\bar{\rho}^{(1/2)}=\frac{(1-\nu)}{\beta}L\tilde{\rho}^{(1/2)}_{\tilde{z}=0},\end{equation}

where $L=\int_{z={{\mathcal{H}}}^{(0)}}\,\mathrm{d}s$ is the instantaneous length of the microscale interface. Equation (4.10) provides a Robin boundary condition for the problem (3.48a), replacing the matching conditions (3.48c).

In particular, we can use this expression for the instantaneous vapour flux (4.10) to argue that the intermediate region is needed in this case of $\alpha=O(\epsilon^{1/2})$ , unlike in the case $\alpha=O(\epsilon)$ above. We demonstrate that the intermediate region is needed using a contradiction, and so suppose that the intermediate region is not required, that is, that the flux through the intermediate region is constant in t and equal to the far-field flux $\rho^{(0)}_Z|_{Z=0}$ . The solution $\tilde{\rho}^{1/2}$ of (3.48) then takes form:

(4.11) \begin{equation} \tilde{\rho}^{(1/2)}=\rho^{(0)}_Z|_{Z=0}\tilde{z}+\bar{\rho}^{(1/2)}(t).\end{equation}

From (4.10), we see that $L\bar{\rho}^{(1/2)}$ must be independent of t, and since for a general pore-scale geometry L varies in t, we therefore require that $\bar{\rho}^{(1/2)}$ is non-constant in t (noting from (4.9) that we must have $\bar{\rho}^{(1/2)}$ non-zero in order for there to be an interface motion at all). This leads to the required contradiction, as $\tilde{\rho}^{(1/2)}$ given by (4.11), with $\bar{\rho}^{(1/2)}$ non-constant in t does not satisfy the diffusion equation (3.48a).

Thus, the instantaneous flux of vapour out of the inner region cannot be independent of t, and so the intermediate region is needed to smooth the temporal oscillations in the vapour flux in this case $\alpha=O(\epsilon^{1/2})$ .

4.2. Case $\alpha=O(1)$

With $\alpha=O(1)$ , the chemistry condition (4.1) relates the leading-order interface velocity to the leading-order vapour density. At O(1), (4.1) reads

(4.12) \begin{equation} 1-\bar{\rho}^{(0)}=-\alpha\frac{H_T+{{\mathcal{H}}}_t^{(0)}}{\sqrt{1+\left({{\mathcal{H}}}^{(0)}_x\right)^2}} \qquad \text{on }z={{\mathcal{H}}}^{(0)}.\end{equation}

We know from our analysis of Section 3.3.3 that $\bar{\rho}^{(0)}$ is independent of x, z and t. Integrating (4.12) over the interface $z={{\mathcal{H}}}^{(0)}$ , and averaging in time, we see

(4.13) \begin{align} (1-\bar{\rho}^{(0)})H_T\int_{t=0}^{1/H_T}\int_{z={{\mathcal{H}}}^{(0)}}\,\mathrm{d}s\,\mathrm{d}t&=-\alpha H_T\int_{t=0}^{1/H_T}\int_{z={{\mathcal{H}}}^{(0)}}\frac{H_T+{{\mathcal{H}}}_t^{(0)}}{\sqrt{1+\left({{\mathcal{H}}}^{(0)}_x\right)^2}}\,\mathrm{d}s\,\mathrm{d}t\nonumber\\ & =-\alpha \phi H_T,\end{align}

using the definition (3.22) of $\mathcal{F}$ and the fact that $\mathcal{F}=\phi$ as shown in Section 3.2.2.

We now define a new effective parameter, $\mathcal{L}$ , by

(4.14) \begin{equation} \mathcal{L}=H_T\int_{t=0}^{1/H_T}\int_{z={{\mathcal{H}}}^{(0)}}\,\mathrm{d}s\,\mathrm{d}t.\end{equation}

We note in particular that although the definition of $\mathcal{L}$ is similar to that of $\mathcal{F}$ (given by (3.22)), the two are not equal: $\mathcal{L}$ is the time-averaged length of the microscale interface (integrating with respect to $\mathrm{d}s$ ), whereas $\mathcal{F}$ was the average flux through the interface (integrating with respect to $\mathrm{d}x$ ). Using the fact that $\bar{\rho}^{(0)}=\rho^{(0)}|_{Z=0}$ from (3.54), we see that (4.13) becomes the effective interface condition:

(4.15) \begin{align} (1-\rho^{(0)})\mathcal{L}=-\alpha\phi H_T \qquad \text{ at }Z=0.\end{align}

This is the chemistry interface condition that closes the homogenised model in this case $\alpha=O(1)$ . We note that the presence of $\mathcal{L}$ in (4.15) means that the effective evaporation rate at the macroscale interface is different (by a factor of $\mathcal{L}/\phi$ ) to that at the microscale interface. For longer, more tortuous interfaces on the microscale, we expect the overall evaporation rate to be increased. This is intuitive, as the greater the interface length, the more liquid can evaporate simultaneously.

We emphasise that, unlike the case $\alpha\ll1$ where the variations in shape and speed of the microscale interface had no effect on the leading-order macroscale model, when $\alpha=O(1)$ we need to understand the motion of the microscale interface $z={{\mathcal{H}}}^{(0)}$ in order to compute the effective parameter $\mathcal{L}$ and thereby close the homogenised model (3.70).

Figure 6. Microscale interface motion for different pore-scale geometries. The interface shape $t=f(x,y)$ , or level sets of the f, are shown in blue at a number of times. The characteristics or rays are straight lines, shown in green. Left: for a diamond-shaped solid obstacle, the solution is given explicitly by (4.19). Centre and right: sketch of two different geometries with the same porosity, but different time period and average interface length $\mathcal{L}$ .

The motion of the microscale interface $z={{\mathcal{H}}}^{(0)}$ is governed by (4.12). Since $\bar{\rho}^{(0)}$ is independent of x, z and t, we may consider $\bar{\rho}^{(0)}$ to be a fixed, externally prescribed constant for the purposes of determining the microscale interface motion. We note that in this case $\alpha=O(1)$ , we may then compute the density correction $\bar{\rho}^{(1)}$ from (3.60), once the interface position ${{\mathcal{H}}}^0(x,t)$ has been determined from (4.12), although we do not need to do this to close the homogenised model.

The problem (4.12) is a non-linear first-order hyperbolic equation for ${{\mathcal{H}}}^0(x,t)$ . The problem is closed by the requirement that ${{\mathcal{H}}}^0$ be periodic in both t and x and is constrained by the solid geometry $\omega_s(t)$ . It is instructive to reformulate this problem for the interface position $z={{\mathcal{H}}}^{(0)}(x,t)$ as a problem for $t=f(x,y)$ , in the stationary coordinate system. The equation (4.12) states that the normal interface velocity is constant and equal to $v\;:\!=\;(1-\bar{\rho}^{(0)})/\alpha$ . In terms of f, (4.12) is therefore equivalent to

(4.16) \begin{equation} v=\frac{1}{\sqrt{f_x^2+f_y^2}},\end{equation}

or

(4.17) \begin{equation} f_x^2+f_y^2=\frac{1}{v^2},\end{equation}

which is the well-known eikonal equation (see, for instance, [19]). In this formulation, we interpret the interface at time t to be the level set $t=f(x,y)$ of the surface f. There is a large literature concerning solutions of the eikonal equation (4.17) motivated, in particular, by problems in optics. Our problem is unusual in that we require our solution f to be periodic in x, and also that the shape of the interface is periodic in time t, which is equivalent to the requirement that there exists a constant $H_T$ such that

(4.18) \begin{equation}f(x,y)+1/H_T=f(x,y+1)\end{equation}

for all x and y. These periodicity constraints in particular make this problem difficult to solve for general pore-scale geometries. However, we note that for special geometries in which there are regions of pore space with no obstacles in the path of the moving interface, we make progress analytically, using Charpit’s method.

For instance, for the diamond geometry shown in Figure 6 (left), the interface moves uniformly downwards unimpeded by any obstacles, as a flat surface, in the regions $|x|>|b|$ . The periodicity condition (4.18) in this case can be replaced by the boundary condition that $f(x,y)=0$ on $y=0$ , $x\in[b,1/2]$ . The solution in this case is

(4.19) \begin{equation} f(x,y)=\begin{cases} -\dfrac{y}{v} & \text{ for }x\in\left[b,\dfrac{1}{2}\right] \text{ and }x\in\left[-\dfrac{1}{2},-b\right]\\[8pt] \dfrac{1}{v}\sqrt{(b-x)^2+y^2} & \text{ for }x\in[-b,b]\setminus\omega_s. \end{cases}\end{equation}

The characteristics form an expansion fan from the points (b, 0) and $(-b,0)$ , and the interface shape in the regions spanned by these expansion fans are arcs of circles. We note there is a shock or ridge along the y-axis where characteristics meet; across this shock, the solution f is continuous. We note that, in this case, the time period of the interface motion, $1/H_T$ , is simply $1/v$ .

Given the solution f of (4.17), the effective interface length $\mathcal{L}$ is given by:

(4.20) \begin{equation} \mathcal{L}=H_T\int_{t=0}^{1/H_T}\int_{t=f(x,y)}\,\mathrm{d}s\,\mathrm{d}t.\end{equation}

For the example with diamond-shaped solid inclusions, we find from (4.15), along with the fact that $H_T=v$ , that $\mathcal{L}=\phi$ (this is true for any geometry where part of the interface can move vertically downwards through the entire cell, unimpeded by solid obstacles). To demonstrate that $\mathcal{L}$ is not always equal to $\phi$ , we consider two pore-scale geometries, illustrated in Figure 6 (centre and right), with the same porosity $\phi$ but where the pore space of one consists simply of a straight vertical channel, while the other is a thinner but more tortuous path. The average interface length $\mathcal{L}$ is larger for the first geometry than for the second, so $\mathcal{L}\neq\phi$ for the second geometry. Furthermore, since the normal velocity of the microscale interface is the same (v) for each case, the time taken to traverse the second, more tortuous cell geometry is longer than that for the first, and so the average vertical velocity $H_T$ is greater for the first cell than for the second.

Finally in this section, we note that the instantaneous vapour flux coming out of the inner region, from (3.64), is

(4.21) \begin{equation} \mathcal{D}\tilde{\rho}^{(1/2)}_{\tilde{z}}|_{\tilde{z}=0}=-\frac{(1-\nu\bar{\rho}^{(0)})(1-\bar{\rho}^{(0)})}{\alpha}\int_{z={{\mathcal{H}}}^{(0)}}\,\mathrm{d}s,\end{equation}

where we have used (4.12) and the fact that $\bar{\rho}^{(0)}$ is independent of the microscale. We note that this instantaneous vapour flux leaving the inner region does indeed vary on the fast-timescale, t, in general, since we have seen that the instantaneous interface length, $\int_{z={{\mathcal{H}}}^{(0)}}\,\mathrm{d}s$ varies in time. Thus, while the intermediate region was not required in the case $\alpha=O(\epsilon)$ , for this case $\alpha=O(1)$ (like the case $\alpha=O(\epsilon^{1/2})$ ) we really do require the intermediate smoothing region in our asymptotic structure.

4.3. Summary of the homogenised model

In summary, we have derived the homogenised model (3.70), along with a chemistry interface condition, which, depending on the size of $\alpha$ , manifests as either a Dirichlet boundary condition for the vapour density $\rho^{(0)}$ , or a relation between the interface velocity $H_T$ and the local vapour density at the interface.

Returning to variables $Y=Z+H$ fixed with respect to the porous medium (rather than moving with the macroscale evaporating interface), our homogenised model is

(4.22a) \begin{align} \frac{\partial q}{\partial Y}&=0,\end{align}

(4.22b) \begin{align} q&=-kP_Y,\end{align}

(4.22c) \begin{align} \delta\phi\rho^{(0)}_T+\nu q\rho^{(0)}_Y&=\frac{\partial}{\partial Y}\left(\mathcal{D}\rho^{(0)}_Y\right),\end{align}

with the interface conditions at $Y=H(T)$ ,

(4.22d) \begin{align} &q=-(1-\delta)\phi H_T,\end{align}

(4.22e) \begin{align} &\mathcal{D}\rho^{(0)}_Y=(1-\nu\rho^{(0)})\phi H_T,\end{align}

(4.22f) \begin{align} &\begin{cases} -\alpha \phi H_T=\mathcal{L}(1-\rho^{(0)}) & \qquad \text{ if }\alpha=O(1),\\[5pt] \rho^{(0)}=1 & \qquad \text{ if }\alpha\leq O(\epsilon^{1/2}). \end{cases}\end{align}

Initial conditions for q, P and $\rho^{(0)}$ , and boundary conditions at the surface of the porous medium close the model.

We note that, for sufficiently small $\alpha\lesssim \epsilon^{1/2}$ , the homogenised model depends on the microstructure only through the bulk effective parameters $\phi$ , $\mathcal{D}$ and k, whereas for $\alpha=O(1)$ there is additional dependence on the average length $\mathcal{L}$ of the microscale evaporating interface. The analysis of this section has also shown that the fine-scale behaviour at, and near, the macroscale interface is fundamentally different for the different ranges of $\alpha$ . For $\alpha=O(\epsilon)$ or smaller, it is the leading-order density gradient in the outer problem that drives the system. The far-field gradient draws vapour off the interface, which must move at an oscillating rate on the microscale in order to maintain the locally constant flux of vapour off the interface. In this way the outer problem informs the inner, and so an intermediate smoothing region is not needed. By contrast, for $\alpha=O(1)$ , the local density value determines the evaporation rate, so that the microscale interface moves with constant normal velocity independent of t. Since the vapour flux leaving the interface therefore oscillates in time, the intermediate region is required to smooth out these oscillations in order to match between the inner and outer regions.

4.4. Correction to previous article

Finally, we give a correction to [16]. The microscale model analysed in [16] has the same form as the evaporation model (2.9), in the special case that there is no fluid flow ( $\delta=1$ , so that ${\textbf{u}}={\textbf{0}}$ ), with a chemical-reaction-rate chemistry condition at the moving interface, directly analogous to the $\alpha=O(1)$ case for our current evaporation setting. The structure of the homogenisation argument in [16] is correct. However, the matching between the inner and intermediate regions is not fully correct; this should be done as in Section 3.3.3 above. Furthermore, the conclusion of [16] that the microscale interface is flat and moving at constant velocity is incorrect; as discussed in Section 4.2, the microscale interface should in fact have a constant normal velocity, which does not result in a flat interface shape. From our analysis above, we see that the correct homogenised model for the reactive decontamination problem studied in [16], replacing the equations (4.74)–(4.75) in that paper, is, for $Y>H(T)$

(4.23a) \begin{align} C_{T}=\frac{\mathcal{D}}{\phi}C_{YY},\end{align}

with, on $Y=H(T)$

(4.23b) \begin{align} CH_T+\frac{\mathcal{D}}{\phi}C_Y&=\beta^*\frac{\mathcal{L}}{\phi}C,\end{align}

(4.23c) \begin{align} H_T&=-\gamma\beta^*\frac{\mathcal{L}}{\phi}C \end{align}

where $\mathcal{D}$ is given by (3.40), $\mathcal{L}$ is given by (4.14) and $\phi$ , the porosity of the medium, is the liquid-occupied area of the cell. The parameters $\beta^*$ and $\gamma$ are as defined in [16]. The form of the boundary conditions given in [16] is therefore accurate, but the effective parameter factors in those conditions are incorrect.

5. Discussion and conclusion

We have derived a homogenised model (4.22) for the motion of an evaporation front through porous media, and describing the flow and advective–diffusive transport of vapour away from the evaporating interface. In particular, we have used a boundary layer analysis combined with homogenisation in both time and space to derive the effective boundary conditions at the evaporating interface to apply to the homogenised model. As is standard for homogenisation problems, the effect of the microstructure enters the bulk equations through the effective parameters $\phi$ , k and $\mathcal{D}$ , which are the porosity, effective permeability and effective diffusivity, respectively. We have shown that the interface conditions describing the conservation of mass across the interface also depend on the microstructure via the effective parameters $\phi$ and $\mathcal{D}$ . In particular, we found that the homogenised interface conditions do not depend on the microscale oscillations in the shape and speed of the evaporating interface. The homogenised model was closed by deriving a homogenised chemistry interface condition. The form of the microscale chemistry interface condition, (characterised by $\alpha$ , ratio of the interface velocity to the rate at which the vapour density reaches its equilibrium) was seen to affect the form of the corresponding effective interface condition on the macroscale. For $\alpha\ll1$ so that the rate at which the vapour reaches its saturation point at the interface is faster than the speed of the evaporating interface, we found a Dirichlet boundary condition on the macroscale for the leading-order vapour density. For the case $\alpha=O(1)$ , the homogenised boundary condition gives the interface velocity in terms of an effective evaporation rate. This effective evaporation rate for the homogenised model was shown to be different in general from the microscale evaporation rate, altered by a factor of $\mathcal{L}/\phi$ , where $\mathcal{L}$ is the average length of the microscale interface. The value of $\mathcal{L}$ may be found for a given pore-scale geometry by solving an eikonal equation for the motion of the microscale interface. From the discussion in Section 4.2, it seems likely that $\phi$ is an upper bound on $\mathcal{L}$ , so that the effective evaporation rate is bounded above by the microscale evaporation rate. Further investigation is needed to establish whether or not this is true. If it is indeed true, then this upper bound on the effective evaporation rate would correspond to a lower bound on the overall drying time of a porous material, or equivalently an upper bound on the rate of offgassing, independent of the microstructure, which is likely to be useful in industrial applications where drying is a production-limiting process, such as pharmaceuticals and ceramics.

In addition to resulting in different effective interface conditions, we also saw different fine-scale solution structures near the interface for different values of $\alpha$ . For $\alpha=O(1)$ , the local vapour density determines the rate of evaporation, forcing the microscale interface to move with a constant normal velocity. Due to the solid inclusions, the flux of vapour from the microscale interface oscillates in time, and we required a wider, intermediate boundary layer of $O(\epsilon^{1/2})$ about the evaporation front, over which these temporal oscillations are smoothed out. By contrast when $\alpha\ll \epsilon^{1/2}$ , we found that the overall gradient in the vapour density determines the rate of evaporation, with vapour being drawn off the evaporating interface at a rate that cannot oscillate in time. Instead the microscale interface must adapt its shape or speed in order maintain the locally constant vapour flux. Since there are no temporal oscillations in this case, no intermediate smoothing layer is needed, and the structure of the asymptotic analysis is greatly simplified. Between these two extremes, we found some elements of each case: for all $\alpha\ll 1$ the effective macroscale interface condition is a Dirichlet condition, but an intermediate smoothing layer is required for $\alpha=O(\epsilon^{1/2})$ .

The parameters $\delta$ and $\nu$ were taken to be order one in the homogenisation analysis. As discussed in Section 2.1, in many circumstances it may be reasonable to assume one or both of these is small. Taking the limit(s) $\delta\rightarrow0$ and or $\nu\rightarrow0$ is commutable with the homogenisation analysis in this paper. For simplicity, we assumed a two-dimensional porous medium, and we also assumed that the macroscale problem was one-dimensional, with unidirectional flow and density gradients on the macroscale, and a flat macroscale interface. We anticipate that our analysis will extend fairly straightforwardly to the case of a flat evaporating interface (surface) within a three-dimensional porous medium. It is less clear how our analysis might extend to a non-flat macroscale interface, since the macroscale interface shape may impact the local shape and velocity of the microscale interface.

By assuming that the liquid was stationary, our analysis precluded effects due to capillarity, surface tension and viscous forces. We determine the average interface length $\mathcal{L}$ purely by the removal of liquid to become vapour, like an ‘etching’ problem, and this approximation may be unrealistic in many practical situations. An important area for future work is to include the flow of liquid, and the associated viscous and capillary effects, into this analysis. While we expect these processes to affect the microscale interface shape and motion, we note that we expect the homogenised model within the gas phase to be the same as in this paper: we anticipate that the fluid flow effects will affect $\mathcal{L}$ , but not the structure of the homogenised evaporation model.

We note that our homogenisation of the advection–diffusion equation was very similar to the analysis of [16], for a purely diffusive system. This was largely because, since the liquid density decreases as it vaporises ( $\delta<1$ ) the flow generated by the evaporation is away from the interface, so that (as in Section 2.1) the lengthscale over which the vapour density varies is that of the entire porous medium. In processes with $\delta>1$ (for instance, a reactive decontamination system such as that in [7, 16], with a denser reaction product than reactant), fluid would flow in towards the interface, and we might anticipate diffusive boundary layer structures at the interface.

Our homogenisation analysis shows how we may incorporate the effect of the pore microstructure into the interface conditions for evaporation in porous media, and in particular quantifies the effect of the microstructure on the evaporation rate, via the effective parameter $\mathcal{L}$ . Although in this paper we modelled an evaporation process, our results are also pertinent to chemical reaction fronts in porous media, where (as in Section 4.4) we likewise see the average interface length $\mathcal{L}$ impacting the effective reaction rate. Crucial future work will be to compute the quantity $\mathcal{L}/\phi$ for a variety of microstructures, in order to quantify the effect that microscale geometry has on the macroscale rate of evaporation and to solve our homogenised model in scenarios of practical interest to those working in both filtration and decontamination.