3.1. Model statement

As derived in [Reference Luckins, Breward, Griffiths and Wilmott10] (and including the minor correction given in [Reference Luckins, Breward, Griffiths and Please11]), the dimensional model for the cleanser concentration $c$ in the sharp-interface case is, for $0\lt y\lt h\lt L_{\text{SI}}$ ,

(3.2a) \begin{align} c_t=D_0c_{yy}, \end{align}

where subscripts $t$ and $y$ denote partial derivatives with respect to time and space, respectively, along with boundary conditions at the reacting interface at $y=h(t)$ :

(3.2b) \begin{align} ch_t+D_0 c_y = -\beta kc, \end{align}

(3.2c) \begin{align} h_t =\chi \beta kc, \end{align}

while at $y=0$ ,

(3.2d) \begin{align} c=c_0. \end{align}

Here, the diffusivity $D_0$ is the effective diffusivity of the cleanser through the (constant porosity) porous material, given by solving a ‘cell problem’ (see [Reference Luckins, Breward, Griffiths and Wilmott10]), $k$ [ms $^{-1}$ ] is the surface reaction rate, $\chi$ [m $^3\,$ mol $^{-1}$ ] is the molar volume of the neat agent and $\beta$ is a dimensionless geometric parameter which is the ratio of the average length of the microscale reacting interface to the porosity of the porous medium that arises from averaging the microscale problem (see [Reference Luckins, Breward, Griffiths and Please11]). The dimensional parameters in our models are summarised in Table 1. We assume that, initially, the agent fills the porous medium, so that $h=0$ at $t=0$ . An initial condition for $c$ is not required since there is no cleanser-occupied region at $t=0$ .

Table 1. Dimensional parameters in the decontamination models (SI is sharp-interface model, AW is agent-on-walls model)

3.2. Nondimensionalisation

We nondimensionalise the model using the scalings

(3.3) \begin{align} c=c_0c^{\prime}, \qquad y=L_{\text{SI}}y^{\prime}, \qquad h=L_{\text{SI}}h^{\prime}, \qquad t=\frac{L_{\text{SI}}}{\chi c_0 \beta k}t^{\prime}, \end{align}

where we have fixed the cleanser concentration scaling, $c_0$ , by the boundary condition at $y=0$ , the length scaling to be the initial depth of the agent spill, and we have chosen the timescale associated with the motion of the interface due to the chemical reaction by balancing (3.2c). This is the timescale of a reaction-limited decontamination, which we will find is the shortest timescale the overall decontamination can take.

Substituting these scalings (3.3) into the model (3.2) and dropping the prime notation, we obtain the dimensionless problem, for $0\lt y\lt h(t)\leq 1$

(3.4a) \begin{align} \alpha \gamma c_t=c_{yy}, \end{align}

with boundary conditions at $y=h(t)$ ,

(3.4b) \begin{align} \alpha \gamma ch_t+c_y=-\gamma c, \end{align}

(3.4c) \begin{align} h_t=c, \end{align}

and on $y=0$ ,

(3.4d) \begin{align} c=1. \end{align}

We have introduced two dimensionless parameters, $\alpha$ and $\gamma$ , defined in Table 2. Since $\chi$ is the molar volume of the neat agent, $\alpha$ characterises the ratio of cleanser concentration in solution to the agent concentration (in its neat form). Since the cleanser is dilute while the agent is neat, it is likely that $\alpha$ is less than one. This justifies our choice of the reaction timescale over the diffusion timescale, since for relatively dilute cleanser ( $\alpha \lt 1$ ) the reaction will often be the controlling process. We note that $\alpha ^{-1}$ may be interpreted as the Stefan number, as defined in [Reference Evans and King6]. Meanwhile, $\gamma$ is a ratio of the lengthscale of the spill to be decontaminated, $L_{\text{SI}}$ , to the lengthscale of diffusion over the reaction timescale. We may therefore consider $\gamma$ to be a proxy for the depth of the spill: $\gamma$ is small for shallow spills and large for deep spills.

Table 2. Dimensionless parameters in the sharp-interface model (3.4)

This problem (3.4) is a simplified version of the decontamination model studied in [Reference Dalwadi, O’Kiely, Thomson, Khaleque and Hall4], in the case that the reaction product is only soluble in the cleanser fluid (or $\mathcal{K}\rightarrow 0$ in the notation of [Reference Dalwadi, O’Kiely, Thomson, Khaleque and Hall4]). As noted in [Reference Dalwadi, O’Kiely, Thomson, Khaleque and Hall4], the problem (3.4) is a one-phase Stefan problem with kinetic undercooling and has been studied in detail, for instance in [Reference Evans and King6]. In the remainder of this section, we solve (3.4) numerically and discuss some relevant limits, which will be helpful when we compare with the agent-on-walls model in Section 4.

3.3. Numerical method and early-time behaviour

To solve (3.4) numerically, we first transform to a fixed domain by making the change of variables $y=hY$ , for $Y\in [0,1]$ , and we then use the method of lines on the transformed problem with 100 spatial meshpoints, which we find to be sufficient for convergence of the solution. For the spatial discretisation, we use a central-difference scheme for the diffusive term and a first-order upwind scheme for the artificial advection due to change of variables to a fixed domain, so that the overall scheme is first order. We then use the inbuilt ODE solver ode15s in Matlab for the timestepping, noting that the system is stiff for $\alpha \gamma \ll 1$ . This is a multistep solver, using numerical differentiation formulas of order 1–5, specifically designed for stiff systems [Reference Shampine and Reichelt14].

The numerical method is not valid at $t=0$ since initially we have $h=0$ . In order to initialise the numerical simulations, we therefore use asymptotic expansions for early time. Specifically, we consider $t\ll 1$ to be small. As in, for instance, [Reference Dalwadi, O’Kiely, Thomson, Khaleque and Hall4, Reference Evans and King6], we see that in order to balance (3.4c) for small $t$ we must have $h$ and $y$ of order $t$ . In (3.4a) we see that the transport of cleanser is fast at early times, with $c_{yy}=0$ . Further, the boundary condition (3.4b) requires that $c_y=0$ at $y=h$ , and thus $c$ is uniform at early times, and equal to 1 using the boundary condition at $y=0$ . At next order in $t$ ( $\ll 1$ ), we find linear equations for the first-order correction terms, which we solve to find

(3.5) \begin{align} c= 1-\gamma (1+\alpha )y+O(t^2), \qquad h=t-\frac{\gamma (1+\alpha )}{2}t^2+O(t^3). \end{align}

These early-time expansions (3.5) are valid so long as $t\ll \max (1,\gamma ^{-1},(\alpha \gamma )^{-1})$ . We use these early-time expansions with appropriately chosen $t$ to initialise numerical simulations of (3.4).

In the remainder of this section, we solve the model both numerically and asymptotically in relevant limits of small and order-one $\alpha$ . For ease of later comparison with the agent-on-walls case, we summarise these asymptotic limits in Figure 2 and refer to them as case i (defined by $\alpha \ll 1$ and with the sublimits cases ia and ib for which $\gamma \ll 1$ and $\gamma \gg 1$ , respectively), case ii ( $\alpha =O(1)$ and $\gamma \ll 1$ ) and case iii ( $\alpha =O(1)$ and $\gamma \gg 1$ ), respectively. We will compare the analytic solutions in these asymptotic limits with full numerical solutions throughout.

Figure 2. Schematic of the asymptotic limits considered for the sharp-interface model. Case i is $\alpha \ll 1$ , with sublimits ia for $\gamma \ll 1$ and ib for $\gamma \gg 1$ . Case ii is $\alpha =O(1)$ and $\gamma \ll 1$ , and case iii is $\alpha =O(1)$ and $\gamma \gg 1$ .

3.4. Case i: dilute cleanser, $\alpha \ll 1$

As discussed in Section 3.2, we expect that a physically relevant limit is $\alpha \ll 1$ , in which the cleanser chemical is dilute in relation to the neat agent. In this section, we therefore study asymptotic solutions of (3.4) in the limit $\alpha \ll 1$ , which we refer to as case i in the schematic of Figure 2. We assume first that $\gamma =O(1)$ , before summarising the behaviour in the sublimits $\gamma \ll 1$ (case ia) and $\gamma \gg 1$ (case ib).

We assume first that $\gamma =O(1)$ . To leading order in $\alpha \ll 1$ , the system (3.4) is simply

(3.6a) \begin{align} c_{yy}=0, \end{align}

with boundary conditions at $y=h$ ,

(3.6b) \begin{align} c_y=-\gamma c, \qquad h_t=c, \end{align}

along with $c=1$ on $y=0$ , at $h=0$ at $t=0$ . This reduced problem has solution

(3.7) \begin{align} c(y,t)=1-\frac{\gamma y}{\sqrt{1+2\gamma t}}, \qquad h=\frac{1}{\gamma }\left ({-}1+\sqrt{1+2\gamma t}\right ), \end{align}

for $y\in [0,h]$ and $h\in [0,1]$ . In particular, the leading-order time at which the decontamination is completed, defined to be the time at which $h=1$ and the reacting interface has moved through the entire spill, is given by

(3.8) \begin{align} t_D=1+\frac{\gamma }{2}. \end{align}

We note that the diffusion of cleanser is quasi-steady in this limit.

This approximate solution (3.7) is valid for the sublimit $\gamma \ll 1$ , case ia, corresponding to a thin spill depth, and we note that in this case $c\approx 1$ at all times, as the cleanser consumed in the reaction is quickly replenished from the surface of the porous medium by diffusion. Additionally, $h\approx t$ in this sublimit, that is, the front moves linearly with time through the material, rather than like $\sqrt{t}$ . The decontamination is limited by the rate of reaction in this case.

We now consider large spill depths, $\gamma \gg 1$ , referred to as case ib in Figure 2. It is not immediately clear that the case-i solution is valid in this limit, since for $\gamma \gg 1/\alpha$ , we might expect the timescale of diffusion to re-enter the problem. We will, however, show that the case-i solution (3.7) is indeed valid for all $\gamma \gg 1$ .

For large $\gamma$ , the case-i solution is valid for early $O(1)$ time. However, the timescale of decontamination must be longer. Making the change of variables $t=\gamma T$ , we see that the model (3.4) becomes

(3.9a) \begin{align} \alpha c_T=c_{yy}, \end{align}

with, on $y=h$ ,

(3.9b) \begin{align} h_T=-\frac{c_y}{1+\alpha c}, \qquad \frac{1}{\gamma }h_T=c, \end{align}

(where we have substituted (3.4c) into (3.4b)), and as usual $c=1$ on $y=0$ . For $\alpha \ll 1$ and $\gamma \gg 1$ (for the leading-order behaviour it is not relevant whether $\alpha \gamma$ is large or small), we see that, at leading order in $\gamma ^{-1}\ll 1$ and $\alpha \ll 1$ , $c$ is linear and equal to zero at $y=h$ . The leading-order interface speed $h_T=-c_y|_{y=h}$ is determined by the flux of cleanser into the interface. The leading-order solution is therefore

(3.10) \begin{align} c=1-\frac{y}{\sqrt{2T}}, \qquad h=\sqrt{2T}, \end{align}

which we note is the sublimit $\gamma \gg 1$ of (3.7). The fact that $\alpha \ll 1$ means that the cleanser concentration profile evolves instantaneously with the interface motion. Indeed both diffusion and reaction are quasi-steady processes in the $\gamma \gg 1$ sublimit; the decontamination is limited by the flux of cleanser into the reacting interface.

The position of the reaction front $h$ against time and the concentration $c$ of cleanser against $y$ for various times through the decontamination process, in case i, given by (3.7), are plotted in Figure 3, along with numerical solutions of (3.4) for various $\alpha$ and $\gamma$ . We see that the cleanser diffuses into the porous medium, with a monotonic decreasing concentration front the surface $y=0$ to the reaction front $y=h(t)$ . The reaction front $h(t)$ progresses into the porous material, with a speed that is non-increasing in time. We observe excellent agreement between the numerical and analytic solutions for the relatively small value of $\alpha =0.1$ . For larger $\alpha =0.5$ , corresponding to more-concentrated cleanser solutions, we see that the dimensionless decontamination time increases, since the diffusive transport of cleanser to the reaction front is no longer quasi-steady. (We note that our chosen timescale, given in (3.3), is inversely proportional to $\alpha$ ; the dimensional decontamination time is considered in Section 3.6 below). For non-zero $\alpha \gamma$ , we also have a nonlinear interface condition (3.4b), which slows the interface motion. We observe in Figure 3b, d and f that $c$ is nearly linear in $y$ , although for larger $\alpha$ and $\gamma$ the profile becomes slightly convex. We note that for larger $\gamma$ there is $O(1)$ variation in $c$ over the domain (for instance in Figure 3b), while for smaller $\gamma$ (smaller spill depths) $c$ remains close to 1 throughout the domain (as in Figure 3f).

Figure 3. Numerical and analytic solutions of the sharp-interface model (3.4), for small $\alpha$ , with $\gamma =5$ (top), $\gamma =1$ (middle) and $\gamma =0.2$ (bottom). On the left we show the position, $h(t)$ , of the reacting interface against time $t$ for various $\gamma$ and for various small values of $\alpha$ , along with the leading-order asymptotic solution (3.7) in the limit of small $\alpha$ . On the right, we show the cleanser concentration profiles $c(y)$ at times $t_D/4$ , $t_D/2$ , $3t_D/4$ and $t_D$ (the decontamination time, at which $h=1$ ), with the arrows showing the direction of increasing time, for the corresponding $\gamma$ and $\alpha$ , and again showing the asymptotic small- $\alpha$ solution (3.7).

3.5. Strong cleanser, $\alpha =O(1)$

With $\alpha =O(1)$ the cleanser applied at the surface of the porous medium has a similar concentration to the (neat) agent. We can make analytical progress in both the limit of shallow ( $\gamma \ll 1$ ) and deep ( $\gamma \gg 1$ ) agent spills, corresponding to cases ii and iii, respectively, in Figure 2. The case $\alpha,\gamma =O(1)$ is intractable analytically, but we will also show numerical solutions in this regime in Figure 5.

3.5.1. Case ii: shallow spills, $\gamma \ll 1$

For shallow spills, with $\gamma \ll 1$ , we look for a solution of (3.4) of the form

(3.11) \begin{align} c=c_0+\gamma c_1+\ldots, \qquad h=h_0+\gamma h_1+\ldots . \end{align}

Substituting these expansions (3.11) into (3.4), we find that, at leading order, $c_0$ and $h_0$ are given by

(3.12) \begin{align} c_0=1, \qquad h_0=t. \end{align}

Continuing to next order in $\gamma$ , we find that

(3.13) \begin{align} c= 1-\gamma (1+\alpha )y+O\left(\gamma ^2\right), \qquad h=t-\frac{\gamma (1+\alpha )}{2}t^2+O\left(\gamma ^2\right), \end{align}

which are in fact the early-time expansions (3.5). Thus, for sufficiently shallow spills, the linear early-time approximation is valid for the entirety of the decontamination process.

We therefore see that, in this shallow-spill case, the reaction front moves with constant velocity to leading order (in agreement with case ia), and that the concentration profile is independent of time at both leading and first order in $\gamma$ . The leading-order solution is shown in Figure 5a and b and discussed further in Section 3.5.2 below.

The decontamination time, at which $h=1$ , is given by

(3.14) \begin{align} t_D=1+\frac{\gamma (1+\alpha )}{2}+O\left(\gamma ^2\right). \end{align}

We note that we regain the dilute-cleanser solution (3.8) in the limit $\alpha \rightarrow 0$ .

3.5.2. Case iii: deep spills $\gamma \gg 1$

In the case of deep agent spills, so that $\gamma \gg 1$ , we note from (3.4a) that the diffusive flux of cleanser through the domain is a limiting process over the long lengthscale of the spill. The decontamination therefore occurs over a longer timescale, and so we make the change of variables

(3.15) \begin{align} t=\gamma T, \end{align}

for $T=O(1)$ , and we write $h(t)=H(T)$ . On this new timescale, the model (3.4) becomes

(3.16a) \begin{align} \alpha c_T&=c_{yy} & &\text{ for }0\lt y\lt H(T), \end{align}

(3.16b) \begin{align} c&=1 & &\text{ at }y=0, \end{align}

(3.16c) \begin{align} H_T&=-\frac{c_y}{1+\alpha c} \text{ and } H_{T}=\gamma c & &\text{ at }y=H(T), \end{align}

where we have combined the boundary conditions at $H(T)$ .

We look for an expansion of the form

(3.17) \begin{align} c=c_0+\frac{1}{\gamma }c_1+O\left (\frac{1}{\gamma ^2}\right ), \qquad H=H_0+\frac{1}{\gamma }H_1+O\left (\frac{1}{\gamma ^2}\right ). \end{align}

Substituting (3.17) into (3.16), we find that, at leading order in $\gamma ^{-1}$ ,

(3.18a) \begin{align} \alpha c_{0T}&=c_{0yy} & &\text{ for }0\lt y\lt H_0(T), \end{align}

(3.18b) \begin{align} c_0&=1 & &\text{ at }y=0, \end{align}

(3.18c) \begin{align} c_{0y}&=-H_{0T} \text{ and } c_0=0 & &\text{ at }y=H_0(T). \end{align}

This leading-order model admits a similarity solution, of the form

(3.19) \begin{align} c_0=1-\frac{\textrm{erf}\left (\frac{\sqrt{\alpha }y}{2\sqrt{T}}\right )}{\textrm{erf}\left (\frac{\sqrt{\alpha }\eta _*}{2}\right )}, \qquad H_0=\eta _*\sqrt{T}, \end{align}

where the constant $\eta _*$ satisfies the algebraic equation

(3.20) \begin{align} 2\sqrt{\frac{\alpha }{\pi }}=\eta _*\exp \left (\frac{\alpha \eta _*^2}{4}\right )\textrm{erf}\left (\frac{\sqrt{\alpha }\eta _*}{2}\right ). \end{align}

Equation (3.20) may be solved numerically for $\eta _*$ for any given $\alpha$ . The solution for a range of $\alpha$ is shown in Figure 4. In particular, we note that as $\alpha \rightarrow 0$ , $\eta _*\rightarrow \sqrt{2}$ , (indeed for $\alpha \lt 1$ we find $\eta _*$ very close to $\sqrt{2}$ ) and so we recover the solution (3.10) for the dilute-cleanser, deep-spill case, in this limit.

Figure 4. The numerically computed solution $\eta _*$ of (3.20), as a function of $\alpha$ .

The leading-order decontamination time is

(3.21) \begin{align} T_D=\frac{1}{\eta _*^2}+O\left (\frac{1}{\gamma }\right ), \end{align}

so that, in terms of our original time variable,

(3.22) \begin{align} t_D=\frac{\gamma }{\eta _*^2}+O(1). \end{align}

In Figure 5, we show solutions of the sharp-interface model (3.4) for $\alpha =1$ , and various values of $\gamma$ . We see excellent agreement between the numerical solutions and both the analytic solution (3.13) in the small- $\gamma$ limit (in Figure 5a and b) and the analytic solution (3.19) in the large- $\gamma$ limit (in Figure 5c and d).

Figure 5. Numerical and analytic solutions of the sharp-interface model (3.4), examining the effect of $\gamma$ , with $\alpha =1$ throughout. In (a) and (b), we show the small- $\gamma$ numerical solutions of (3.4), with the leading-order small- $\gamma$ asymptotic solution (3.12). In (c) and (d), we show large- $\gamma$ solutions with the leading-order large- $\gamma$ asymptotic solution (3.19). In (a) and (c), we show the motion of the reaction front $h(t)$ . In (b) and (d), we show the cleanser concentration profiles at times $t_D/4$ , $t_D/2$ , $3t_D/4$ and $t_D$ (the decontamination time, at which $h=1$ ), with arrows indicating increasing time. In (b), we mark the end of each line with a circle for the cases $\gamma =0.1$ and $\gamma =0.01$ , as otherwise the curves are difficult to distinguish.

3.6. Decontamination times and decontamination efficiency

There are two practical metrics we can evaluate to assess the quality of the decontamination process. These are:

(i) the decontamination time, which is the time $t_D$ at which all the agent has been reacted away, that is, when $h(t_D)=1$ . In the various asymptotic limits above, we found the analytical expressions (3.8), (3.14) and (3.22) for the leading-order decontamination time in cases i, ii and iii, respectively.

(ii) the efficiency of the decontamination process, in terms of the amount of cleanser chemical used. We would ideally like to minimise the excess cleanser introduced to the system, since (although less toxic than the agent) the cleanser can also be harmful to people or to the environment. Since we have assumed the reaction has stoichiometry 1 (one mole of agent reacts with one mole of cleanser), the minimum number of moles of cleanser that are required to neutralise the entire agent spill (per unit surface area of the spill) is

(3.23) \begin{align} N=\frac{v}{\chi }\quad \text{mol}\,\text{m}^{-2}. \end{align}

The amount of cleanser per unit cross-sectional area that enters the porous medium over the entire decontamination process, however, is

(3.24) \begin{align} \text{total flux in}= - \frac{D_0\phi }{\chi \beta k}\int _{t=0}^{t_D}c_y|_{y=0}\,\mathrm{d}t \quad \text{mol}\,\text{m}^{-2}\,. \end{align}

We define the decontamination efficiency to be the ratio of the minimum amount of cleanser required to the flux of cleanser into the porous medium, namely

(3.25) \begin{align} \text{efficiency}=\frac{N}{\text{total flux in}}=\frac{\gamma }{-\int _{t=0}^{t_D}c_y|_{y=0}\,\mathrm{d}t}\,, \end{align}

where we have made use of (3.1) and (3.23). The efficiency defined by (3.25) is a dimensionless value between 0 and 1, with 1 being perfectly efficient decontamination. We note that we may compute the leading-order decontamination efficiencies in the various asymptotic limits studied above and find

(3.26) \begin{align} \text{efficiency}=\begin{cases}1+O(\alpha ) &\text{ for case i, }\alpha \ll 1,\\[5pt] \dfrac{1}{1+\alpha }+O(\gamma ) & \text{ for case ii, }\alpha =O(1),\gamma \ll 1,\\[9pt] \exp \left ({-}\dfrac{\alpha \eta _*^2}{4}\right )+O\left (\gamma ^{-1}\right ) &\text{ for case iii, }\alpha =O(1),\gamma \gg 1, \end{cases} \end{align}

where $\eta _*(\alpha )$ is the solution of (3.20).

In Figure 6, we show the variation of the decontamination time (left) and decontamination efficiency (right) with the system parameters and compare numerical solutions of (3.4) to these limiting analytic expressions. In Figure 6a, we show the variation of the numerically computed decontamination times as a function of $\gamma$ . In particular, since $\gamma$ is proportional to $L_{\text{SI}}$ , we may interpret this plot as the effect of the spill depth $L_{SI}$ on the decontamination time (since for a particular agent, cleanser and porous material, $\beta k$ and $D_0$ are constant). For small $\alpha$ , we observe that $t_D$ behaves like the case-i solution, $1+\gamma/2$ , as given by (3.8). For larger $\alpha =O(1)$ the limiting behaviours (3.14) and (3.22) capture the behaviour of the numerical solution well. In particular, we see that, for both $\alpha \ll 1$ and $\alpha =O(1)$ , the dimensionless decontamination time appears to increase linearly with spill depth. Since the time scaling used in (3.3) is also linear in $L_{\text{SI}}$ , we therefore conclude that the dimensional decontamination time increases quadratically with spill depth $L_{\text{SI}}$ .

Figure 6. Sharp-interface model decontamination time $t_D$ or scaled decontamination time $t_D/\alpha \gamma$ (left) and decontamination efficiency (right) as the system parameters $\alpha$ and $\gamma$ are varied. Numerical solutions of (3.4) are shown in colour, and the analytical expressions (3.8), (3.14) and (3.22) for the decontamination time, and (3.26) for the efficiency, in the various asymptotic limits, are shown in black (dashed, dotted or solid).

We note that the time scaling used in (3.3) also depends on both the reaction rate $k$ and the applied cleanser concentration $c_0$ . In order to draw physical conclusions about the impact of these parameters on the decontamination time, in Figure 6c and e we therefore plot the scaled decontamination time

(3.27) \begin{align} \frac{t_D}{\alpha \gamma }=\frac{D_0}{\chi c_0\beta k L_{\text{SI}}}t_D, \end{align}

which has the same dependence on $c_0$ and $ k$ as the time scaling used in (3.3). Thus, the dependence of $t_D/\alpha \gamma$ on $ k$ is the same as the dependence of the dimensional decontamination time on $k$ . In Figure 6c, we see that this scaled decontamination time decreases with $\gamma$ so that (as we might expect) faster reaction rates $k$ result in faster decontamination. However, we observe limiting returns, with the scaled decontamination time approaching an ( $\alpha$ -dependent) constant for larger values of $\gamma$ . This is because as $\gamma$ increases the reaction becomes almost instantaneous and the problem is no longer reaction-limited, but instead becomes limited by the supply of cleanser. In Figure 6e, we show the variation of the scaled decontamination time $t_D/\alpha \gamma$ with $\alpha$ (which is proportional to the applied cleanser concentration $c_0$ , and so we may interpret increasing $\alpha$ with increasing the applied cleanser concentration). For larger $\alpha$ , and so stronger cleanser concentrations, we observe faster decontamination. As $\alpha$ becomes small, the scaled decontamination time grows like $\alpha ^{-1}$ : the small- $\alpha$ expression, which from (3.8) is given by $\alpha ^{-1}(\gamma ^{-1}+1/2)$ , appears to be a good approximation for $\alpha \lt 1$ .

In Figure 6b and d, we show the variation of the decontamination efficiency (3.25) with $\gamma$ and $\alpha$ , respectively. For this sharp-interface case, the efficiency can also be thought of as a measure of how much cleanser is left inside the porous medium, in $y\in [0,1]$ , at time $t_D$ . We see in Figure 6b that the efficiency increases with $\gamma$ , since we recall (from Figure 5d, say) that, for larger $\gamma$ , we expect lower values of $c$ , and thus less cleanser is contained in $y\in [0,1]$ . The variation with $\gamma$ is most noticeable for larger $\alpha$ . For $\alpha \ll 1$ , the efficiency approaches 1. This is because for small applied cleanser concentrations $c_0$ , there are fewer moles of cleanser chemical left within the porous material at the decontamination time. From Figure 6d, we see that the small- $\gamma$ and large- $\gamma$ expressions (3.26) capture the decrease of the efficiency with $\alpha$ well. Moreover, we observe that the numerically simulated efficiencies shown in Figure 6b and d appear to be bounded below and above by the small- and large- $\gamma$ asymptotic formulae, respectively.

In summary, we have seen that increasing the reaction rate $\beta k$ decreases the decontamination time (albeit with limited returns) and also improves the efficiency of the system. Increasing the concentration of the applied cleanser $c_0$ results in faster decontamination times but reduced decontamination efficiency. In practice, therefore, there is a pay-off between cleanser wastage and speed of decontamination. To illustrate this, we show the variation of the efficiency with the decontamination time (parameterised by $\alpha$ or $c_0$ ) in Figure 6f. The applied cleanser concentration should be chosen to strike a reasonable balance between fast decontamination times and limited waste of the cleanser.

4.1. Model statement

In [Reference Luckins, Breward, Griffiths and Wilmott10], we showed that in regions with $R=0$ the cleanser concentration $c$ satisfies

(4.1a) \begin{align} c_t=D_0c_{yy}, \end{align}

while if $R\gt 0$ we have the reaction–diffusion problem

(4.1b) \begin{align} c_t&=\frac{1}{\mathcal{V}(R)}\left (\mathcal{D}(R)\mathcal{V}(R)c_y\right )_y-\mathcal{F}(R)k\left (1+\chi c\right )c, \end{align}

(4.1c) \begin{align} R_t&=-\chi k c, \end{align}

where $\mathcal{D}(R)$ is the diffusivity of the cleanser through the available pore space when the agent thickness is $R$ (this is found by the numerical solution of a cell problem which is stated in Appendix A and solved in, for instance, [Reference Bruna and Chapman1]), and $D_0=\mathcal{D}|_{R=0}$ , $\mathcal{V}(R)=1-\pi (r_0+R)^2/d^2$ , is the fraction of the cell volume occupied by the cleanser solution, and $\mathcal{F}(R)=2\pi (r_0+R)/d^2\mathcal{V}(R)$ is the ratio of the length of the agent–cleanser interface on the microscale to the volume of the cell occupied by cleanser.

Initially, we suppose that $c=0$ is uniform throughout $y\gt 0$ , and that the agent is uniformly distributed with initial layer thickness $R=R_0$ , through the region $0\lt y\lt L_{\text{AW}}$ , where

(4.2) \begin{align} L_{\text{AW}}=\frac{vd^2}{\pi \left(2R_0r_0+R_0^2\right)}, \end{align}

is the spill lengthscale. We note that, since we want the cleanser-occupied region to be initially connected, we require $R_0+r_0\lt d/2$ . Since the agent layer is being removed, the cleanser region is then connected for all time. In higher dimensions or for different microscale geometries, we would require a similar constraint on the initial agent-layer thickness to ensure this connectedness. For all $t\gt 0$ , as in the sharp-interface case, we impose a constant concentration $c=c_0$ of cleanser at $y=0$ , but as $y\rightarrow \infty$ we require $c\rightarrow 0$ .

4.2. Nondimensionalisation

In addition to the spill lengthscale $L_{AW}$ , there are various additional lengthscales in this model, compared with the sharp-interface model, namely the initial agent-layer thickness $R_0$ , the radius of the circular solid inclusions, $r_0$ , and the separation distance between these circular solid inclusions, $d$ . We will assume that $\rho _0\,:\!=\,r_0/d=O(1)$ (a property of the porous material only) is not particularly small. However, we will allow that ratio

(4.3) \begin{align} \delta =\frac{R_0}{d}, \end{align}

to be $O(1)$ or smaller, so that the layer of agent may be thin relative to the pore space. We also define the ratio of the pore lengthscale to the spill lengthscale as

(4.4) \begin{align} \epsilon =\frac{d}{L_{\text{AW}}}. \end{align}

We require $\epsilon \ll 1$ in order that our homogenised model (4.1) is valid, and so we will assume this to be true throughout this paper.

We use the scalings

(4.5a) \begin{align} c=c_0\hat{c}, \qquad R=R_0\hat{R}, \qquad y=L_{\text{AW}}\hat{y}, \qquad t=\frac{R_0}{\chi k c_0}\hat{t}, \end{align}

where we have chosen the timescale that balances (4.1c). The functions $\mathcal{D},\mathcal{V}$ and $\mathcal{F}$ are simply functions of geometry, although $\mathcal{D}$ must be calculated numerically by solving the cell problem given in Appendix A (as in, for instance, [Reference Bruna and Chapman1, Reference Dalwadi, Griffiths and Bruna3]). We make the scalings

(4.5b) \begin{align} \mathcal{D}&=D_0\hat{D}\left(\rho _0+\delta \hat{R}\right), \end{align}

(4.5c) \begin{align} \mathcal{V}&=V_0\hat{V}\left(\hat{R};\delta,\rho _0\right)=V_0\left (1-\delta \frac{\pi \left(2\rho _0\hat{R}+\delta \hat{R}^2\right)}{V_0}\right ), \end{align}

(4.5d) \begin{align} \mathcal{F}&=F_0\hat{F}\left(\hat{R};\delta,\rho _0\right)=F_0\left (1+\frac{\delta \hat{R}}{\rho _0}\right )\left (\frac{V_0}{1-\pi \left(\rho _0+\delta \hat{R}\right)^2}\right ), \end{align}

where $D_0$ , $V_0=1-\pi \rho _0^2$ and $F_0=2\pi \rho _0/dV_0$ are the values of $\mathcal{D}$ , $\mathcal{V}$ and $\mathcal{F}$ , respectively, when the agent thickness $\hat{R}=0$ .

Substituting (4.5) into (4.1) and dropping the hat notation, our model becomes

(4.6a) \begin{align} A\Gamma c_t&=\frac{1}{V(R)}\left (D(R)V(R)c_y\right )_y-\Gamma H(R)F(R)\left (1+\frac{\delta A}{f_1} c\right )c, \end{align}

(4.6b) \begin{align} R_t&=-H(R)c, \end{align}

in $y\gt 0$ , where $H(R)$ is the Heaviside function. The boundary conditions are

(4.6c) \begin{align} c=1 \qquad &\text{ at } \qquad y=0, \end{align}

(4.6d) \begin{align} c\rightarrow 0 \qquad &\text{ as } \qquad y\rightarrow \infty, \end{align}

while at $t=0$ ,

(4.6e) \begin{align} c=0 \quad \text{for }y\gt 0, \qquad \text{ and } \qquad R=\begin{cases}1 & \text{if }y\in [0,1],\\ 0 & \text{otherwise}. \end{cases} \end{align}

The dimensionless parameters in the model are the lengthscale ratios $\delta$ and $\epsilon$ , given by (4.3) and (4.4), and $\Gamma$ and $A$ , given by

(4.7) \begin{align} \Gamma =\frac{\gamma f_2}{\epsilon \delta }=\frac{kL_{\text{AW}}}{\epsilon f_1D_0}, \qquad A=\frac{\alpha f_1}{\delta }=\frac{c_0\chi f_1}{\delta }, \end{align}

where the sharp-interface parameters $\gamma$ and $\alpha$ are as defined in Table 2. The model also depends on the $O(1)$ geometric factors

(4.8) \begin{align} f_1=\frac{1-\pi \rho _0^2}{2\pi \rho _0}, \qquad f_2=\frac{2 \rho _0}{\beta (2\rho _0+\delta )}. \end{align}

The dimensionless parameters in (4.6) are summarised in Table 3. We plot $D$ , $V$ and $F$ as functions of $R$ for given $\rho _0$ and $\delta$ in Figure 7, using $D(R)$ as computed in [Reference Dalwadi, Griffiths and Bruna3]. We note that $V$ and $D$ decrease with $R$ , while $F$ increases with $R$ .

Table 3. Dimensionless parameters in the agent-on-walls model (4.6)

Figure 7. The variation of (a) $V$ and $D$ and (b) $F$ with agent thickness $R$ . We take $\rho _0=0.2$ throughout, so that the model is valid for $\delta \in (0,0.3)$ . Three values $\delta =0.01$ (dash-dotted curves), $\delta =0.1$ (dashed curves) and $\delta =0.3$ (solid curves) are shown.

4.3. Comparison with the sharp-interface model

Before analysing solutions of (4.6), we first compare the model structure, scalings and dimensionless parameters with that of the sharp-interface model (3.4).

Firstly, because we are comparing the models where the initial volume of agent is the same, we note that the depth of porous medium to be decontaminated is smaller in the sharp-interface case than the agent-on-walls case, and $L_{\text{SI}}/L_{\text{AW}}=O(\delta )$ . This is because the agent occupies only a fraction of the pore space in the agent-on-walls case, rather than saturating it. The timescales of the decontamination are likewise different with $[t_{\text{SI}}]/[t_{\text{AW}}]=O(\delta/\epsilon )$ , (where $[t_{\text{SI}}]$ and $[t_{\text{AW}}]$ are the time scalings used in (3.3) and (4.5), respectively). We recall that we require $\epsilon \ll 1$ for our agent-on-pores model to be valid, while $\delta$ may be small or of order one. If $\delta/\epsilon \gg 1$ , then the decontamination timescale is much shorter for the agent-on-walls case than for the sharp-interface case. This is despite the fact that the spill is deeper in the agent-on-walls case than for the sharp-interface case. The reason the agent-on-walls decontamination is faster is because the surface area of the interfaces between the fluids is much greater in the agent-on-pores case, and so the overall decontamination can occur much more quickly. We note, however, that with our choice of timescales we have assumed that the decontamination is limited by the reaction rate, rather than by diffusion of the cleanser. In each case, this choice is valid for sufficiently shallow agent spills but for deeper spills, where the transport of cleanser to the reacting interface(s) is also important, the timescale of decontamination is longer. Indeed for deep spills we will see that the decontamination timescale may actually be longer for the agent-on-walls model than for the sharp-interface case.

The parameters $A$ and $\Gamma$ in the agent-on-walls model are closely linked to $\alpha$ and $\gamma$ (defined in Table 2), respectively, in the sharp-interface model; the differences are purely geometric scaling factors. Specifically, $A$ and $\alpha$ in the agent-on-walls and sharp-interface models, respectively, characterise the amount of cleanser in the system compared with the local amount of agent to be removed. In the sharp-interface case, this is simply $\alpha =c_0\chi$ (essentially a ratio of the chemical concentrations), but in the agent-on-walls case we have $A=\alpha f_1/\delta$ , since the volume fraction of agent within the porespace is reduced by a factor of $\delta/f_1$ . Since we expect that $c_0\lesssim \chi$ , we only consider $\alpha$ to be $O(1)$ or smaller in the sharp-interface case. However, in the agent-on-walls case, $A$ may be small, order one or large, depending on the relative sizes of $\alpha$ and $\delta$ : if $\delta \ll \alpha$ so that $A\gg 1$ , we will find in Section 4.7 that diffusion is the limiting process of the decontamination. Both $\gamma$ and $\Gamma$ are a measure of the depth of the spill or speed of the reaction: they are the ratio of the spill depth to the diffusion lengthscale of the cleanser, over the reaction timescale.

As for the sharp-interface model, we can make asymptotic progress in certain parameter regimes. These are illustrated schematically in Figure 8, and referred to as cases I–III, which, as we will show, bear some resemblance to the limits i–iii in the sharp-interface case. We note that, for all of our asymptotic analyses of the agent-on-walls model, we will use $\delta \ll 1$ .

Figure 8. Schematic of the asymptotic limits considered for the agent-on-walls model. Case I is $A\ll 1$ , with sublimits $\Gamma \ll 1$ and $\Gamma \gg 1$ referred to as cases Ia and Ib, respectively. Case II is $A\gg 1$ and $\Gamma \ll 1$ , while case III is $A\gg 1$ and $\Gamma \gg 1$ . For all these asymptotic limits, we additionally assume $\delta \ll 1$ .

4.4. Numerical method

As for the sharp-interface model (3.4), we solve the model (4.6) using the method of lines, with 100 spatial meshpoints (beyond which the results are indistinguishable), implemented in Matlab with the inbuilt Ordinary Differential Equation (ODE) solver ode15s for the timestepping. To discretise the spatial derivatives of $c$ , we use a central finite difference scheme.

We note that the lengthscale of the diffusion of cleanser in $y\gt 1$ (where $R=0$ ) may vary drastically from the lengthscale of the agent-occupied region which has been nondimensionalised to reside between $y=0$ and $y=1$ , depending on the size of the parameter group $A\Gamma$ and also on the time $t$ . Specifically, the distance that cleanser has travelled beyond $y=1$ is on the order of $\sqrt{t/A \Gamma }$ , which may be large or small depending on the parameter regime of interest. In order that our numerical scheme is accurate and efficient in both regions $y\in [0,1]$ and $y\gt 1$ , we therefore discretise (4.6) in the two regions $y\lt 1$ and $y\gt 1$ separately and impose conditions of continuity of cleanser concentration and cleanser flux between the two (at $y=1$ ). For $y\lt 1$ , we use a uniform mesh and impose the boundary condition (4.6c) straightforwardly. For $y\gt 1$ , we make the change of variables

(4.9) \begin{align} y=1+\sqrt{\frac{t}{A\Gamma }}\eta, \qquad c(y,t)=\bar{c}(\eta,t). \end{align}

The cleanser concentration $\bar{c}$ in $y\gt 1$ (or $\eta \gt 0$ ) then satisfies

(4.10) \begin{align} t\bar{c}_t-\frac{\eta }{2}\bar{c}_\eta = \bar{c}_{\eta \eta }. \end{align}

The two conditions imposed at $y=1$ ( $\eta =0$ ) to couple the two systems together are continuity of cleanser concentration and flux, namely

(4.11) \begin{align} c(1,t)=\bar{c}(0,t), \qquad D(R)V(R)c_y(1,t)= \sqrt{A\Gamma } \frac{\bar{c}_\eta (0,t)}{\sqrt{t}}. \end{align}

We solve this problem (4.10) and (4.11) on the domain $\eta =[0,5]$ , imposing the Dirichlet condition (4.6d) at $\eta =5$ (which we find to be sufficiently large for numerical accuracy). To do this, we discretise the equation (4.10) with central differences for the diffusion term, and a first-order upwind (forward) scheme for the artificial advection term due to the growing coordinate system, to ensure stability. (For the practical implementation, we rescale $\eta$ onto $[0,1]$ and so use the same mesh for $\eta$ as for $y$ .) We note that the apparent singularity at $t=0$ in (4.10) and (4.11) does not cause numerical difficulty, since $\bar{c}=\bar{c}_\eta =0$ when $t=0$ .

4.5. Preliminary remarks

4.5.1. The first time at which $R=0$ , and the emergence of a decontamination front

We first consider the early-time behaviour of (4.6) at $y=0$ , where we have $c=1$ . Initially, $R=1$ for $y\geq 0$ and, while $R\gt 0$ everywhere, from (4.6b) we have

(4.12) \begin{align} R_t=-c. \end{align}

Thus $R_t|_{y=0}=-1,$ and so $R=1-t$ while $t\lt 1$ at $y=0$ . The first occurrence of $R=0$ is therefore at time $t=1$ , regardless of the size of the dimensionless parameters. This justifies our choice of timescale used in (4.5), as this is the minimum time that the decontamination can take.

We define the position, $y=y_*(t)$ to be the largest value of $y$ for which $R=0$ , so that $R=0$ for $y\lt y_*$ and $R\gt 0$ for $y\gt y_*$ . This position is analogous to the position of $y=h(t)$ of the reacting interface in the sharp-interface model, and we view it as the decontamination front. In particular, the decontamination time, that is, the time when the decontamination process is complete, occurs at the time $t_D$ at which $y_*=1$ . We note that $y_*=0$ when $t=1$ , unlike the reaction front $y=h(t)$ for the sharp-interface case, which satisfies $h=0$ at $t=0$ .

4.5.2. Behaviour for $y\gt 1$

For $y\gt 1$ where $R=0$ for all time, we may make analytical progress. Here, $c$ satisfies the diffusion equation

(4.13) \begin{align} A\Gamma c_t=c_{yy}, \end{align}

with the requirement that $c=0$ at $t=0$ and $c\rightarrow 0$ as $y\rightarrow \infty$ . This problem has solution (using the Green’s function approach, [Reference Keener8])

(4.14) \begin{align} c(y,t)=\sqrt{\frac{A\Gamma }{4 \pi }}\int _{\tau =0}^{t}\frac{c(1,\tau )(y-1)}{(t-\tau )^{3/2}}\exp \left({-}\frac{A\Gamma (y-1)^2}{4(t-\tau )}\right)\,\mathrm{d}\tau, \end{align}

in terms of the value of $c$ at $y=1$ at all previous times. As well as continuity at $y=1$ , we require that the solution has a continuous diffusive flux of cleanser at $y=1$ . Since $R$ is discontinuous at $y=1$ (while $R\gt 0$ for $y\lt 1$ ), the diffusivity and pore-volume occupied by cleanser are discontinuous. Continuity of diffusive flux requires

(4.15) \begin{align} D(R)&V(R)c_y|_{y=1^-}=c_y|_{y=1^+}\nonumber \\ &=\sqrt{\frac{A\Gamma }{4 \pi }}\,\lim _{y\rightarrow 1}\left (\int _{\tau =0}^{t}\frac{c(1,\tau )}{(t-\tau )^{3/2}}\exp \left({-}\frac{A\Gamma (y-1)^2}{4(t-\tau )}\right)\left (1-\frac{A\Gamma (y-1)^2}{2(t-\tau )}\right )\,\mathrm{d}\tau \right ),\end{align}

This expression (4.15) will be helpful in our analytic analyses to come, as it quantifies the size of the flux of cleanser leaving the agent-occupied region through $y=1$ , namely it is on the order of $\sqrt{A\Gamma/ t}$ .

4.5.3. The limit $\delta \ll 1$

In order to make analytic progress in the following sections, we regularly consider the limit $\delta \ll 1$ , in which the initial agent-layer thickness $R_0$ is small relative to the separation between the solid circles making up the microstructure. With $\delta \ll 1$ ,

(4.16) \begin{align} D(R)=1+O(\delta ), \qquad V(R)=1+O(\delta ), \qquad F(R)=1+O(\delta ), \end{align}

are all constant at leading order in $\delta$ , and – so long as $A\delta \ll 1$ – the model equations (4.6a)–(4.6b) reduce at leading order to

(4.17a) \begin{align} A\Gamma c_t&=c_{yy}-\Gamma H(R)c, \end{align}

(4.17b) \begin{align} R_t&=-H(R)c. \end{align}

4.6. Case $A\ll 1$ , reaction-limited decontamination

In this section, we suppose $\alpha \ll \delta$ so that $A\ll 1$ . This is the dilute-cleanser limit, since $\delta$ is order 1 or smaller, and so we have $\alpha \ll 1$ . As in the sharp-interface case $\alpha \ll 1$ , we will see that for $A\ll 1$ equation (4.6a) is dominated by diffusion (regardless of the size of $\Gamma$ ), so that diffusion is a quasi-steady process.

4.6.1. Case I: analytical solutions in the limit $\delta \ll 1$

In order to make analytical progress, we further assume that $\delta \ll 1$ . In this limit, where $A\ll 1$ and $\delta \ll 1$ (case I as shown in the schematic of Figure 8), the leading-order model is (4.17) in $y\gt 0$ , with boundary and initial conditions (4.6c)–(4.6e). We note that in the analysis below we will replace the far-field boundary condition (4.6d) with the condition (4.15), and solve (4.17) on $y\in [0,1]$ only. In particular, we note that the problem (4.6c) for $c$ decouples from the problem for $R$ (so long as we know the region for which $R\gt 0$ ).

We now explore solutions of (4.17) for $A\ll 1$ , for various sizes of $\Gamma$ . We begin with case I ( $\Gamma =O(1)$ ) before then looking at the sublimit $\Gamma \ll 1$ which is case Ia. We then find additional structures in case Ib for which $\Gamma \gg 1$ .

With $\Gamma =O(1)$ , to leading order in $A\ll 1$ and $\delta \ll 1$ we have

(4.18a) \begin{align} c_{yy}&=\Gamma H(R)c,\\[-10pt]\nonumber \end{align}

(4.18b) \begin{align} R_t&=-H(R)c. \end{align}

The flux of cleanser through $y=1$ is given by (4.15) and is of order $\sqrt{A\Gamma }$ . Thus, to leading order in $\sqrt{A}\ll 1$ , we replace (4.6d) with the boundary condition

(4.18c) \begin{align} c_y=0 \qquad \text{ on }y=1. \end{align}

This, along with the boundary condition (4.6c), and the initial conditions (4.6e), closes the leading-order problem, and we need not consider $y\gt 1$ explicitly.

For $t\lt 1$ when $R\gt 0$ throughout $y\in (0,1)$ , we find that (4.18) has solution

(4.19) \begin{align} c=\frac{\cosh \left(\sqrt{\Gamma }(1-y)\right)}{\cosh \left(\sqrt{\Gamma }\right)}, \qquad R=1-\frac{\cosh \left(\sqrt{\Gamma }(1-y)\right)}{\cosh \left(\sqrt{\Gamma }\right)}t. \end{align}

We note, in particular, that $c$ is independent of $t$ , since the reaction–diffusion process is quasi-steady, and for $\delta \ll 1$ the depletion of agent has no effect on the transport of $c$ at leading order. (We also note that there is an additional early-time problem when $t=O(A)$ that describes how $c$ evolves to its quasi-steady profile, but which is not otherwise of interest.)

For $t\gt 1$ , we also have a decontaminated region $y\lt y_*(t)$ for which $R=0$ , in addition to a region $y\in [y_*(t),1]$ for which $R\gt 0$ . Solving (4.18a) for $c$ in $y\lt y_*(t)$ and $y\gt y_*(t)$ separately, and imposing the boundary conditions at $y=0$ and $y=1$ , as well as continuity of both $c$ and $c_y$ at $y=y_*(t)$ , we find the solution

(4.20) \begin{align} c=\begin{cases} 1-\dfrac{\sqrt{\Gamma }y}{\coth \left(\sqrt{\Gamma }(1-y_*)\right)+\sqrt{\Gamma }y_*} & \text{ if }y\in [0,y_*(t)],\\[18pt] \dfrac{\cosh \left(\sqrt{\Gamma }(1-y)\right)}{\cosh \left(\sqrt{\Gamma }(1-y_*)\right)+\sqrt{\Gamma }y_*\sinh \left(\sqrt{\Gamma }(1-y_*)\right)} & \text{ if }y\in [y_*(t),1], \end{cases} \end{align}

given in terms of $y_*(t)$ . The leading-order equation (4.18b) for $R$ is therefore

(4.21) \begin{align} R_t=-\frac{\cosh \left(\sqrt{\Gamma }(1-y)\right)}{\cosh \left(\sqrt{\Gamma }(1-y_*(t))\right)+\sqrt{\Gamma }y_*(t)\sinh \left(\sqrt{\Gamma }(1-y_*(t))\right)}, \end{align}

for $y\gt y_*(t)$ , which must be solved along with the initial condition $R=R_1(y)\,:\!=\,R(1,y)$ at $t=1$ from (4.19), and the fact that $R=0$ at $y=y_*(t)$ . We look for a separable solution to (4.21) of the form

(4.22) \begin{align} R(t,y)=R_1(y)-S(t)U(y). \end{align}

Substituting (4.22) into (4.21), we see that

(4.23) \begin{align} U(y)=\cosh \left(\sqrt{\Gamma }(1-y)\right), \end{align}

and $S$ satisfies

(4.24) \begin{align} S^{\prime}(t)=\frac{1}{\cosh \left(\sqrt{\Gamma }(1-y_*)\right)+\sqrt{\Gamma }y_*\sinh \left(\sqrt{\Gamma }(1-y_*)\right)} . \end{align}

Furthermore, in order that $R=0$ at $y=y_*$ , we must have that

(4.25) \begin{align} S(t)=\frac{R_1(y_*(t))}{U(y_*(t))}=\textrm{sech}\left(\sqrt{\Gamma }(1-y_*)\right)-\textrm{sech}\left(\sqrt{\Gamma }\right) . \end{align}

We note that $S=0$ when $y_*=0$ , as required for the initial condition. Substituting (4.25) into (4.24) gives an ODE for $y_*(t)$ , namely

(4.26) \begin{align} y^{\prime}_*(t)=\frac{\coth \left(\sqrt{\Gamma }(1-y_*)\right)}{\sqrt{\Gamma }\left (1+\sqrt{\Gamma }y_*\tanh \left(\sqrt{\Gamma }(1-y_*)\right)\right )} . \end{align}

This ODE can be solved numerically, with the initial condition $y_*=0$ at $t=1$ . The leading-order solutions for $c$ and $R$ are then given in terms of $y_*$ by (4.20) and

(4.27) \begin{align} R=1-\frac{\cosh \left(\sqrt{\Gamma }(1-y)\right)}{\cosh \left(\sqrt{\Gamma }(1-y_*)\right)} \qquad \text{ for }y\in [y_*(t),1] . \end{align}

We present solutions to (4.26) for various $\Gamma$ in Figures 9 and 11. Solutions with $\Gamma =1$ are shown in Figure 9a, c and e, along with numerical solutions of (4.17) with $A=10^{-3}$ and $\delta \ll 1$ , with good agreement. For $t\gt 1$ , we see in Figure 9a that the speed of the decontamination front $y^{\prime}_*$ is increasing. We note from Figure 9c that initial transients in $c$ are not captured by the quasi-steady expression (4.20), but that this is a good approximation thereafter, both for $t\lt 1$ and $t\gt 1$ . The agent-layer thickness, $R$ , shown in Figure 9e decreases from its initial value of 1 and decreases most quickly in regions with higher cleanser concentration (close to $y=0$ ). The analytic solution (4.27) for $R$ is in good agreement with the numerical solution, although the error appears to increase with time. We note that the largest source of error in the small- $A$ approximation is the $O(\sqrt{A})$ error in the flux boundary condition applied at $y=1$ (since the other error from the time-derivative term in (4.17) is of $O(A)$ ); we accordingly see much better agreement at $y=0$ than at $y=1$ , and better agreement at earlier times. We note that the solution (4.27) for $R$ bears some resemblance to a travelling wave. However, it is not a travelling wave in the traditional sense: as well as moving with the decontamination front $y_*(t)$ , the profile $R$ changes with time, decreasing in amplitude. Furthermore, the wave-speed $y^{\prime}_*(t)$ is not constant in time: the wave-speed may initially decrease for $t\gt 1$ , with $y^{\prime\prime}_*(t)\lt 0$ (this is true for the large- $\Gamma$ case considered below), but $y^{\prime}_*$ is always seen to be increasing as $y_*$ approaches the end of the contaminated region at $y=1$ .

Figure 9. Numerical and asymptotic solution of (4.6) in the case $A\ll 1$ , $\delta \ll 1$ and for $\Gamma =1$ (left) and $\Gamma =0.1$ (right). We show the decontamination front $y_*(t)$ (top), cleanser profiles $c$ (middle) and agent-layer thickness $R$ (bottom). The numerical solutions (coloured curves) are computed in the limit $\delta \ll 1$ , making the assumption that $D=V=F=1$ , and taking $A=10^{-3}$ . In (b)–(f), solutions are shown at times $0.05t_D$ , $t_D/4$ , $t_D/2$ , $3t_D/4$ and $0.95t_D$ (where $t_D$ is the decontamination time, at which $y_*=1$ ), and arrows show the direction of increasing time. Analytical solutions are shown in black. In (c) and (d), thick dotted black curves are the analytical solution (4.20) and in (e) and (f) thick dashed curves are (4.27), where $y_*(t)$ is the (numerically computed) solution of (4.26). The thin dashed lines in (f) are the approximation (4.31) to $R$ , and the thick black line $c=1$ in (d) is the approximation (4.31) to $c$ in the small- $\Gamma$ sublimit.

We may also make further analytical progress in case Ia for which $\Gamma \ll 1$ . This is a sublimit of case I ( $\Gamma =O(1)$ ) considered above. For $\Gamma \ll 1$ we see that (4.26) reduces to

(4.28) \begin{align} y^{\prime}_*(t)\sim \frac{1}{\Gamma (1-y_*)}, \end{align}

with solution

(4.29) \begin{align} y_*(t)\sim 1-\sqrt{1+\frac{2}{\Gamma }(1-t)}. \end{align}

This expression is valid only for $1\lt t\lt t_D$ where

(4.30) \begin{align} t_D=1+\frac{\Gamma }{2}, \end{align}

is the (approximate) decontamination time in this limit. We further note that, when $\Gamma \ll 1$ , the analytic solution given by (4.20) and (4.27) reduces to simply

(4.31) \begin{align} c=1, \qquad R=1-t, \end{align}

at leading order, for both $t\lt 1$ and $t\in [1,t_D]$ .

In Figure 9 (right), we compare this solution (4.31) in the $\Gamma \ll 1$ sublimit to the full asymptotic solution for $\Gamma =O(1)$ , given by (4.20) and (4.27), and to numerical solutions of (4.17) in the limit $\delta \rightarrow 0$ , for small $A=10^{-3}$ , taking the value $\Gamma =0.1$ . In Figure 9b, we see that the small- $\Gamma$ approximation (4.29) for the decontamination front $y_*(t)$ agrees well with the numerical solution and the full solution of (4.26). The small- $\Gamma$ solutions (4.31) also approximate the $c$ and $R$ profiles well for $\Gamma =0.1$ , since $c$ is close to 1 and at each time the $R$ -profile is nearly flat.

We note that the limit $A\ll 1$ is analogous to the limit $\alpha \ll 1$ in the sharp-interface model. In particular, we see that, for small $\Gamma$ (analogous to small $\gamma$ in the sharp-interface model), the form of the decontamination time $t_D$ given by (4.30) is closely related to (3.14), in the limit $\alpha \ll 1$ of the sharp-interface model, despite the very different motion of the decontamination front in each case.

Case Ib: the limit $\Gamma \gg 1$

The case $\Gamma \gg 1$ corresponds to the deep-spill limit and is referred to as case Ib in Figure 8. We might expect that the time derivative in (4.17a) would become important, but – as in Section 3.4 for the sharp-interface case – we will see that this is not the case. For $\Gamma \gg 1$ . we will show that the timescale of the overall decontamination is longer than the reaction timescale, by a factor of $\Gamma$ . We explore the decontamination process in this limit in detail by examining the behaviour separately at ‘early’ times (which is the reaction timescale previously considered), at $O(\Gamma )$ times, and near the ‘end’ time through a temporal boundary layer analysis at the end of the decontamination. A schematic of the solution structure in this case Ib is given in Figure 10.

Figure 10. An illustration of the asymptotic solution structure for case Ib. Left: the decontamination front $y_*(t)$ against time, showing the three regions “early” ( $t=O(1)$ ), “late” ( $t=O(\Gamma )$ ) and “end” ( $1-y_*(t)=O(1/\sqrt{\Gamma })$ ) times. Right: the cleanser concentration $c$ and agent-layer thickness $R$ during the late-time stage.

To understand the early-time behaviour, we first consider the sublimit of (4.17) for which $\Gamma \gg 1$ , on an $O(1)$ timescale. Initially, $R\gt 0$ throughout $y\in [0,1]$ , and so, at leading order, we see from (4.17a) and (4.6e) that we must have $c=0$ and therefore, from (4.17b), that $R=1$ . This is incompatible with the boundary condition $c=1$ at $y=0$ , and so there must be a boundary layer near $y=0$ in which $c$ and $R$ vary. Making the change of variables

(4.32) \begin{align} y=\frac{Y}{\sqrt{\Gamma }}, \end{align}

equation (4.17a) becomes (at leading order in $\delta/\Gamma$ )

(4.33) \begin{align} A c_t=c_{YY}- H(R)c, \end{align}

so that, at leading order, we have

(4.34) \begin{align} c_{YY}=H(R)c, \end{align}

along with (4.17b) describing the evolution of $R$ . As $Y\rightarrow \infty$ , we must have $c\rightarrow 0$ and $R\rightarrow 1$ to match with the far-field behaviour. For $t\lt 1$ (when $R\gt 0$ everywhere in $Y\in [0,\infty )$ ), we find the solution

(4.35) \begin{align} c=\exp ({-}Y), \qquad R=1-t\exp ({-}Y). \end{align}

For $t\gt 1$ , we also have a decontaminated region $Y\lt Y_*(t)$ for which $R=0$ . Solving for $c$ in $Y\lt Y_*(t)$ and $Y\gt Y_*(t)$ separately, and imposing continuity of both $c$ and $c_Y$ at $Y=Y_*(t)$ , we find the solution

(4.36) \begin{align} c(Y,t)&=\begin{cases} 1-\dfrac{Y}{1+Y_*(t)} & \text{if }Y\lt Y_*(t),\\[15pt] \dfrac{\exp ({-}(Y-Y_*(t)))}{1+Y_*(t)} & \text{if }Y\gt Y_*(t). \end{cases} \end{align}

In $Y\gt Y_*(t)$ , the agent thickness $R$ then satisfies

(4.37) \begin{align} R_t=-\frac{\exp ({-}(Y-Y_*(t)))}{1+Y_*}. \end{align}

As in the $\Gamma =O(1)$ case, we look for a separable solution $R=1-U(t)S(Y)$ . Following the same method, we find

(4.38) \begin{align} S(Y)=\exp ({-}Y), \qquad U(t)=\frac{1}{S(Y_*)}=\exp \left (Y_*(t)\right ), \end{align}

so that

(4.39) \begin{align} R(Y,t)=1-\exp \left ({-}(Y-Y_*(t))\right ) \qquad \text{ for }Y\in [Y_*(t),\infty ). \end{align}

with $Y_*$ the solution of

(4.40) \begin{align} Y^{\prime}_*(t)=\frac{1}{1+Y_*(t)}. \end{align}

We note that this ODE for $Y_*$ is much simpler than the equivalent equation (4.26) in the case $\Gamma =O(1)$ , where effects from the end of the agent-occupied domain impacted on the cleanser concentration profile. The explicit solution to (4.40) with initial condition $Y_*=0$ at $t=1$ is

(4.41) \begin{align} Y_*(t)=-1+\sqrt{2t-1} \qquad \text{ for }t\gt 1. \end{align}

For $t\gt 1$ , the solution is therefore

(4.42a) \begin{align} c&=\begin{cases}1-\dfrac{Y}{Y_*(t)+1} &\text{if }Y\lt Y_*(t),\\[10pt] \dfrac{1}{Y_*(t)+1}\exp \left ({-}(Y-Y_*(t))\right ) &\text{if }Y\gt Y_*(t),\end{cases}\end{align}

(4.42b) \begin{align} R&=\begin{cases} 0 &\hspace{0.6cm}\text{if }Y\lt Y_*(t),\\[5pt] 1-\exp \left ({-}(Y-Y_*(t))\right ) &\hspace{0.6cm}\text{if }Y\gt Y_*(t).\end{cases}\end{align}

Unlike the $\Gamma =O(1)$ case given by (4.20)–(4.27), the wave-like solution $R$ does not change shape as it moves into the domain, although the speed of the wave motion, $Y^{\prime}_*(t)$ , is still non-uniform.

We now consider $O(\Gamma )$ time. We note that, for large $t$ , the early-time decontamination front given by (4.41) behaves like $Y_*\sim \sqrt{2t}$ , so that the position $y_*=Y_*/\sqrt{\Gamma }$ of the decontamination front becomes $O(1)$ after an $O(\Gamma )$ time. To examine the decontamination of the entire spill, we therefore change to an $O(\Gamma )$ timescale.

Making the change of variables $t=\Gamma T$ , equations (4.17a)–(4.17b) become (at leading order in $\delta$ )

(4.43a) \begin{align} A c_T&=c_{yy}-\Gamma H(R)c,\\[-5pt]\nonumber \end{align}

(4.43b) \begin{align} R_T&=-\Gamma H(R)c. \end{align}

The position, $y=y_*(T)$ , of the decontamination front is now of $O(1)$ . For $y\lt y_*$ we see that, to leading order (in $A$ and $\delta$ ), $c$ is linear, and given by

(4.44) \begin{align} c=1-B(T)y \end{align}

for some $B(T)$ to be determined.

For $y\gt y_*(T)$ , where $R\gt 0$ , we see that (4.43a) requires $c=O(1/\Gamma )$ to be small. We therefore again look for an $O(1/\sqrt{\Gamma })$ boundary layer in which $c=O(1)$ , this time centred on the moving interface, by making the change of variables

(4.45) \begin{align} y=y_*(T)+\frac{Y}{\sqrt{\Gamma }}, \qquad \text{for }Y\gt 0. \end{align}

This boundary layer structure is illustrated in Figure 10 (right). Substituting (4.45) into (4.43), we see that, in this boundary layer at $y=y_*(T)$ ,

(4.46a) \begin{align} A \left (c_T-\sqrt{\Gamma }y^{\prime}_*(T)c_Y\right )&=\Gamma \left (c_{YY}- H(R)c\right )\\[-5pt]\nonumber \end{align}

(4.46b) \begin{align} R_T-\sqrt{\Gamma }y^{\prime}_*(T)R_Y&=-\Gamma H(R)c. \end{align}

In particular, we see from (4.46b) that $c$ must be small within this boundary layer and of $O(1/\sqrt{\Gamma })$ . Setting $c=C(Y,T)/\sqrt{\Gamma }$ in (4.46), we find that, at leading order in $O(1/\sqrt{\Gamma })$ ,

(4.47) \begin{align} C_{YY}=C \qquad \text{ for }Y\gt 0, \end{align}

with solution

(4.48) \begin{align} C=D(T)e^{-Y}, \end{align}

in order that $C\rightarrow 0$ as $Y\rightarrow \infty$ . The two integration constants $B(T)$ and $D(T)$ are fixed by imposing continuity of both the cleanser concentration and its derivative at $y=y_*$ , corresponding to $Y=0$ . Specifically, we require $c=0$ and $c_y=C_Y$ at $y=y_*$ ( $Y=0$ ), which fix

(4.49) \begin{align} B(T)=D(T)=\frac{1}{y_*(T)}. \end{align}

In summary, and returning to the original spatial variable $y$ , the leading-order solution, $c$ , is therefore

(4.50) \begin{align} c(y,T)=\begin{cases} 1-\dfrac{y}{y_*(T)} &\text{ if }y\lt y_*(T),\\[10pt] \dfrac{1}{\sqrt{\Gamma }y_*(T)}\exp ({-}\sqrt{\Gamma }(y-y_*(T))) &\text{ if }y\gt y_*(T). \end{cases} \end{align}

The leading-order equation for $R$ is therefore

(4.51) \begin{align} y^{\prime}_*(T)R_Y=\frac{1}{y_*(T)}\exp ({-}Y), \end{align}

which has the wave-like solution

(4.52) \begin{align} R=1-\exp ({-}Y)=1-\exp \left (\sqrt{\Gamma }(y-y_*(T))\right ) \qquad \text{ with } \qquad y_*(T)=\sqrt{2T}. \end{align}

We note this solution is precisely the large- $t$ limit of the early-time solution (4.42); in particular the shape of $R$ is exactly the same as (4.42b). Indeed, the early-time solution (4.42) is a higher-order approximation of the late-time behaviour than (4.50) and (4.52), since there is no new physics entering the problem for $O(\Gamma )$ times. This late-time solution is illustrated schematically on the right of Figure 10.

As an aside, we notice that the leading-order solution given by (4.50) and (4.52) satisfies a reaction-limited, deep-spill, sharp-interface problem, namely

(4.53a) \begin{align} c_{yy}&=0 & \text{ for }y\lt y_*(T), \end{align}

(4.53b) \begin{align} c&=0 &\text{ at }y=y_*(T), \end{align}

(4.53c) \begin{align} y^{\prime}_*(T)&=-c_y, &\text{ at }y=y_*(T), \end{align}

with $c=1$ at $y=0$ , which is precisely the form of the sharp-interface model in case ib. Thus, for $A\ll 1$ , $\delta \ll 1$ , and $\Gamma \gg 1$ (case Ib), the agent-on-pores model behaves precisely like a sharp-interface model (in case ib).

The leading-order decontamination time predicted by (4.52), in the large- $\Gamma$ limit, is $t_D=\Gamma/2$ . Since, as we have already noted, the early-time solution (4.42) is, in fact, valid so long as $y_*$ is an $O(1)$ distance from 1, we may use (4.42) to find the first few correction terms to the decontamination time. Setting $Y_*=\sqrt{\Gamma }$ in (4.41), we find that the decontamination time is

(4.54) \begin{align} t_D^*=\frac{\Gamma }{2}+\sqrt{\Gamma }+1. \end{align}

However, we have assumed in our analysis that the decontamination front is an $O(1)$ distance away from the end of the agent spill at $y=1$ . When $1-y_*(T)$ becomes the same size as the width of the boundary layer, $O(1/\sqrt{\Gamma })$ , we will see end effects similar to those seen in the $\Gamma =O(1)$ case studied above, and we expect that this will provide an $O(\sqrt{\Gamma})$ correction to the decontamination time predicted by (4.54). Indeed, we study the decontamination process in this end-time temporal boundary layer near to $t_D$ , in Appendix B, and find that the $O(\sqrt{\Gamma})$ correction to the decontamination time is $-\sqrt{\Gamma }$ . The improved decontamination time is therefore

(4.55) \begin{align} t_D&=t_D^*+\sqrt{\Gamma }\tau _D^*=\left (\frac{\Gamma }{2}+\sqrt{\Gamma }+1\right )-\sqrt{\Gamma }\nonumber \\ &=\frac{\Gamma }{2}+1.\end{align}

Surprisingly, this is the same as in the small- $\Gamma$ sublimit (4.30), as well as analogous to the sharp-interface model (3.8) (case i).

In summary, we have found that, in the sublimit $\Gamma \gg 1$ , the leading-order solution of (4.17) is, for $t\lt 1$ ,

(4.56a) \begin{align} c=\exp ({-}\sqrt{\Gamma }y), \qquad R=1-t\exp ({-}\sqrt{\Gamma }y), \end{align}

while for $1\lt t\lesssim \Gamma/2-\sqrt{\Gamma }$ ,

(4.56b) \begin{align} c&=\begin{cases} 1-\dfrac{\sqrt{\Gamma }y}{\sqrt{2t-1}} &\hspace{1.8cm}\text{ if }y\lt y_*=\dfrac{1}{\sqrt{\Gamma }}\left ({-}1+\sqrt{2t-1}\right ),\\[20pt] \dfrac{\exp ({-}\sqrt{\Gamma }y-1+\sqrt{2t-1})}{\sqrt{2t-1}} &\hspace{1.8cm}\text{ if }y\gt y_*=\dfrac{1}{\sqrt{\Gamma }}\left ({-}1+\sqrt{2t-1}\right ), \end{cases}\end{align}

(4.56c) \begin{align} R&=\begin{cases} 0 &\text{ if }y\lt y_*=\dfrac{1}{\sqrt{\Gamma }}\left ({-}1+\sqrt{2t-1}\right ),\\[20pt] 1-\exp ({-}\sqrt{\Gamma }y-1+\sqrt{2t-1}) &\text{ if }y\gt y_*=\dfrac{1}{\sqrt{\Gamma }}\left ({-}1+\sqrt{2t-1}\right ). \end{cases}\end{align}

For completeness, we also note that, from the analysis in Appendix B, the leading-order solution near the end-time, for $\Gamma/2-\sqrt{\Gamma }\lesssim t\leq t_D=\Gamma/2+1$ , as

(4.56d) \begin{align} c&=\begin{cases} 1-y+\dfrac{y}{\sqrt{\Gamma }}\left (\dfrac{\coth (Z_*)-Z_*}{1-Z_*/\sqrt{\Gamma }}\right ) &\text{ if }y\lt y_*=1-\dfrac{Z_*}{\sqrt{\Gamma }},\\[15pt] \dfrac{1}{\sqrt{\Gamma }}\dfrac{\left (\exp \left(\sqrt{\Gamma }(1-y)\right)+\exp ({-}\sqrt{\Gamma }(1-y))\right )}{e^{Z_*}-e^{-Z_*}} &\text{ if }y\gt y_*=1-\dfrac{Z_*}{\sqrt{\Gamma }}, \end{cases}\end{align}

(4.56e) \begin{align} R&=\begin{cases} 0 &\text{ if }y\lt y_*=1-\dfrac{Z_*}{\sqrt{\Gamma }},\\[15pt] 1-\dfrac{\exp \left(\sqrt{\Gamma }(1-y)\right)+\exp ({-}\sqrt{\Gamma }(1-y))}{e^{Z_*}+e^{-Z_*}} &\text{ if }y\gt y_*=1-\dfrac{Z_*}{\sqrt{\Gamma }}, \end{cases}\end{align}

where $Z_*(\tau )=Z_*((t-t_D^*)/\sqrt{\Gamma })$ is the solution of (B.9).

In Figure 11, we show solutions of (4.17) in the limit $\delta \rightarrow 0$ , for small $A=10^{-3}$ , and fairly large $\Gamma$ . We compare the numerical solution of (4.17) (with $\delta =0$ ) to the full asymptotic solution, (4.20) and (4.27), as well as to the reduced asymptotic solutions in the large- $\Gamma$ sublimit, given by (4.56). For brevity, here and throughout the remainder of the paper, we show $c$ and $R$ on the same axes. We note that, as in the small- $\Gamma$ case presented in Figure 9, the full asymptotic solution is in very good agreement with the numerical solution. The asymptotic solutions for the large- $\Gamma$ sublimit also capture the profiles $R$ and $c$ well but, in particular, are good approximations for the decontamination time, $t_D$ , even for $\Gamma =20$ , which is not particularly large. We note that we have chosen a very small value of $A$ ( $=10^{-3}$ ) here; this is not necessary to see good agreement between the asymptotic and numerical solutions, but by doing so we can more clearly see the errors in the asymptotic approximations that are due to the $O(\Gamma ^{-1/2})$ or $O(\Gamma ^{-1})$ correction terms.

Figure 11. Numerical and asymptotic solution of (4.6) in the case $A\ll 1$ , $\delta \ll 1$ and for $\Gamma =20$ . The numerical solutions (solid curves) are computed in the limit $\delta \ll 1$ , making the assumption that $D=V=F=1$ , and taking $A=10^{-3}$ . In (a), we show the decontamination front $y_*(t)$ , with the full, early-time and end-time asymptotic approximations of the solution of (4.26), (4.41) and (B.9), respectively. In (b), the solid lines are the numerical solution $c$ (red) and $R$ (blue) at times $0.01t_D$ , $0.1t_D$ , $t_D/4$ , $t_D/2$ , $3t_D/4$ and $0.95t_D$ , where $t_D$ is the time at which $y_*(t_D)=1$ , and the arrows indicate the direction of increasing time. Black dashed ( $R$ ) and dotted ( $c$ ) curves are the full analytical solutions (4.20) and (4.27), where $y_*(t)$ is the (numerically computed) solution of (4.26). The magenta dashed ( $R$ ) and dotted ( $c$ ) curves are the approximations (4.56) in the large- $\Gamma$ sublimit.

4.6.2. Thick agent layers: case $\delta =O(1)$

The asymptotic analysis presented in Section 4.6.1 is in the limit $\delta \ll 1$ , so that the layer of agent $R_0$ initially coating the walls of the pore is thin relative to the pore lengthscale $d$ . In that limit $\delta \ll 1$ , we saw that the equations (4.6a)–(4.6b) for $c$ and $R$ partially decouple, since variations in the thickness of the agent layer do not affect the effective cleanser diffusivity $D$ , or the geometric factors $F$ and $V$ at leading order.

Figure 12. Variation of solutions with $\delta$ in the limit $A\ll 1$ , for $\Gamma =0.01$ (top), $\Gamma =1$ (middle) and $\Gamma =100$ (bottom). On the left, we show the decontamination front $y_*(t)$ and on the right the cleanser $c$ (red) and agent-layer thickness $R$ (blue) profiles. Numerical solutions of (4.6) are shown as solid lines with $\delta =0.25$ , and in the limit $\delta \ll 1$ as dotted lines, for comparison. Throughout we take $A=0.01$ , and $\rho _0=0.2$ . Solutions in (b), (d) and (f) are shown at the times $0.1t_D$ , $0.2t_D$ , $0.4t_D$ , $0.6t_D$ , $0.8t_D$ and $t_D$ , where $t_D$ is the decontamination time for the $\delta =0.25$ solution.

We now investigate the effect of non-negligible $\delta$ in the limit of dilute cleanser $A\ll 1$ . In Figure 12, we show solutions of (4.6) for $\delta =0.25$ (solid lines) compared with solutions in the limit $\delta \ll 1$ (dotted lines), for three different values of $\Gamma$ . We note that due to the geometrical constraints we require $\delta \in (0,0.5-\rho _0)=(0,0.3)$ . We see qualitatively similar behaviours of both the decontamination front $y_*(t)$ and the profiles of $c$ and $R$ to those observed in the small- $\delta$ limit, in all of the cases $\Gamma =0.01$ , $\Gamma =1$ and $\Gamma =100$ . This is unsurprising, as, although the local values of $D$ , $F$ and $V$ vary with $\delta$ , the asymptotic balances are not expected to change until $D$ becomes very small, which occurs when $\delta \rightarrow 0.5-\rho _0$ . However, we notice that, since $D$ now decreases with $R$ , cleanser cannot diffuse so easily into regions with large amounts of agent, and so we observe lower $c$ values in regions with large $R$ for the $\delta =0.25$ solutions, compared with the $\delta \ll 1$ solutions at the same time-points. Since $c$ is reduced, the rate of reaction is slower in these regions, and so the agent-depletion rate is reduced, leading to steeper spatial gradients in the $R$ -profiles. These effects are particularly evident when $\Gamma$ is large, since $R$ exhibits greater spatial variation than for $\Gamma \ll 1$ .

We note that we have kept $\Gamma$ constant in each of the cases shown in Figure 12, in order to observe the effect of the spatial variation in $D$ , $V$ and $F$ . However, since

(4.57) \begin{align} \Gamma =\frac{vk}{\epsilon f_1D_0\pi }\cdot \frac{1}{\delta (2\rho _0+\delta )}, \end{align}

in order to hold the total initial amount of agent, $v$ , constant, we should vary $\Gamma$ simultaneously with $\delta$ . In this case, the qualitative changes as $\delta$ is varied are mainly due to the corresponding change in $\Gamma$ and are similar to those observed in the $\delta \ll 1$ limit previously studied. We therefore do not show these solutions here.

4.7. The case $A\gg 1$ : diffusion-limited decontamination

As discussed previously, we expect the dilution of the cleanser to be $\alpha = O(1)$ or smaller for physically relevant decontamination scenarios, but since $A=\alpha f_1/\delta$ we may have $A\gg 1$ if $\delta \ll \alpha$ , that is, if the layer of agent coating the walls of the porous medium is sufficiently thin. In this section, we suppose $A\gg 1$ and, since $\delta$ is therefore small, the model (4.6) reduces to

(4.58a) \begin{align} A\Gamma c_t&=c_{yy}-\Gamma H(R)\left (1+\frac{\delta A}{f_1}c\right )c,\\[-10pt]\nonumber \end{align}

(4.58b) \begin{align} R_t&=-H(R)c, \end{align}

with the initial and boundary conditions given by (4.6c)–(4.6e).

Since $A\gg 1$ , we immediately notice that the time derivative of $c$ in (4.58a) dominates over the reaction term. Thus, we expect that the system is diffusion-limited when $A\gg 1$ (compared to reaction-limited as we found for the case $A\ll 1$ ). We consider the behaviour in subcases depending on the size of $\Gamma$ , from small $\Gamma =O( A^{-1})\ll 1$ to large $\Gamma \gg 1$ .

Figure 13. Asymptotic behaviour of the decontamination front $y_*$ and decontamination time $t_D$ in case II. In (a), we show the solution $\sqrt{A\Gamma }y_*$ of (4.60c) against $t$ (black) and the large- $t$ (large- $A\Gamma$ ) limit (4.62) (red dashed). In (b), we show the solution $t_D$ of (4.61) (blue) and the large $A\Gamma$ limit (4.63) (red dotted).

4.7.1. Case II: shallow spills, $\Gamma = O(A^{-1})$

We first suppose the spill depth is shallow, taking the distinguished limit $\Gamma =O (A^{-1})\ll 1$ .

With $A\Gamma =O(1)$ , we see that, at leading order in $A^{-1}$ , (4.58a) becomes

(4.59) \begin{align} A\Gamma c_t=c_{yy} \end{align}

regardless of whether $R$ is zero or positive. The cleanser and agent equations therefore fully decouple. Solving (4.59) subject to (4.6c)–(4.6e), the leading-order cleanser concentration $c$ is given by the similarity solution

(4.60a) \begin{align} c=\text{erfc}\left (\frac{\sqrt{A\Gamma }y}{2\sqrt{t}}\right ). \end{align}

Substituting (4.60a) into (4.58b) and solving the resulting equation, we find that $R$ takes the form

(4.60b) \begin{align} R=1-\text{erfc}\left (\frac{\sqrt{A\Gamma }y}{2\sqrt{t}}\right )\left (\frac{A\Gamma y^2}{2}+t\right )+\frac{y\sqrt{A\Gamma t}}{\sqrt{\pi }}\exp \left ({-}\frac{A\Gamma y^2}{4t}\right ), \end{align}

for $y\in [y_*(t),1]$ (and $R=0$ otherwise), with $y_*(t)$ defined as the point at which $R$ is first zero, namely the solution of

(4.60c) \begin{align} 1+\frac{y_*\sqrt{A\Gamma t}}{\sqrt{\pi }}\exp \left ({-}\frac{A\Gamma y_*^2}{4t}\right )=\text{erfc}\left (\frac{\sqrt{A\Gamma }y_*}{2\sqrt{t}}\right )\left (\frac{A\Gamma y_*^2}{2}+t\right ). \end{align}

We show how $\sqrt{A\Gamma } y_*$ , given by the solution of (4.60c), varies with $t$ in Figure 13a, and note that (4.60c) has a solution $y_*\gt 0$ only for $t\gt 1$ . In particular, the decontamination time $t_D$ (at which $y_*=1$ ) is the solution of

(4.61) \begin{align} 1+\sqrt{\frac{A\Gamma t_D}{\pi }}\exp \left ({-}\frac{A\Gamma }{4t_D}\right )=\text{erfc}\left (\frac{\sqrt{A\Gamma }}{2\sqrt{t_D}}\right )\left (\frac{A\Gamma }{2}+t_D\right ). \end{align}

We show how $t_D$ , given by the solution of (4.61), varies with $A\Gamma$ in Figure 13b. This has error of order $\max (A^{-1},\delta )$ (since it is the reaction term of (4.58a), and the $O(\delta )$ corrections to $D,F,$ and $V$ that we have neglected. This solution in the distinguished limit $A\Gamma =O(1)$ is in fact valid for all $\Gamma \ll 1$ , since the loss of cleanser to the chemical reaction only affects the leading-order solution of (4.58a) when $\Gamma =O(1)$ or larger.

For very shallow spills $\Gamma \ll A^{-1}$ , the solutions $c$ and $R$ are approximately uniform in $y$ , and indeed from (4.61) we see that $t_D\rightarrow 1$ as $A\Gamma \rightarrow 0$ . For deeper spills and so larger $\Gamma$ (with $A^{-1}\ll \Gamma \ll 1$ ), we see from (4.58a) that the decontamination process is slower, since the diffusive timescale is longer. Indeed, for $A\Gamma \gg 1$ (with $A^{-1}\ll \Gamma \ll 1$ ) we find that the solution $y_*$ of (4.60c) behaves like

(4.62) \begin{align} y_*\sim \sqrt{\frac{6t}{A\Gamma }W_0\left (\frac{2t^{2/3}}{3\pi ^{1/3}}\right )}, \end{align}

where $W_0$ is the principal branch of the Lambert-W function, and from (4.61) we see that

(4.63) \begin{align} t_D\sim \frac{A\Gamma }{10 W_0\left (\frac{2}{5}\left (\frac{A\Gamma }{4\sqrt{\pi }}\right )^{2/5}\right )}. \end{align}

In Figure 14 (left), we show numerical solutions of (4.58) and the asymptotic approximations (4.60), for $\Gamma =0.1\in (A^{-1},1)$ , and we see very good agreement even for this moderate value of $\Gamma$ . Unlike in the case $A\ll 1$ studied in Section 4.6, the cleanser profile is independent of the agent thickness, with no change in behaviour between $y\lt 1$ and $y\gt 1$ .

Figure 14. Numerical and asymptotic solution of (4.6) in the case $A\gg 1$ (and $\delta =O(A^{-1})\ll 1$ ). We take $A=100$ throughout with $\Gamma =0.1$ (left) and $\Gamma =100$ (right). We show the $c$ and $R$ profiles against the spatial coordinate $y$ in (a)–(b) and against the similarity variable $\sqrt{A\Gamma } y/2\sqrt{t}$ in (c)–(d). The decontamination front $y_*(t)$ is shown in (e)–(f). The numerical solutions (solid curves) are computed in the limit $\delta \ll 1$ , making the assumption that $D=V=F=1$ . In (a)–(d), the solid curves are the numerical solution of (4.58) ( $c$ red, $R$ blue) and the black curves are the analytical solution (4.60a) for $c$ (dotted) and (4.60b) for $R$ (dashed), where $y_*(t)$ is the (numerically computed) solution of (4.60c), shown in (e)–(f), along with the numerical decontamination front.

4.8. Moderate cleanser concentration, $A=O(1)$

Having studied the model (4.6) in the limits of $A\ll 1$ and $A\gg 1$ we now investigate the case $A=O(1)$ . Since $A=\alpha f_1/\delta$ , the case $A=O(1)$ requires that $\alpha/\delta =O(1)$ , so that the local amounts of cleanser and agent are approximately balanced. Specifically, this is true for both $\alpha,\delta$ small (so that the cleanser is dilute, but also there is a small volume fraction of agent within the pore space), or equally for both $\delta$ and $\alpha$ of order 1. We show numerical solutions of (4.6) with $A=1$ in Figure 15, for various sizes of $\Gamma$ , and also for both $\delta \ll 1$ and $\delta =0.25$ . With $A=O(1)$ , we see that both the diffusion of cleanser and the reaction rate play a role, and the system is neither diffusion-limited nor reaction-limited. As in the case $A\ll 1$ , we see qualitative agreement between the cases $\delta \ll 1$ and $\delta =O(1)$ , although the non-uniform cleanser diffusivity $D(R)$ for $\delta =O(1)$ results in slower transport of $c$ through regions where $R$ is close to 1, and correspondingly steeper spatial profiles of $R$ , as well as a clear discontinuity in $c_y$ at $y=1$ .

Figure 15. Variation of solutions with $\delta$ , holding $A=1$ , for $\Gamma =0.01$ (top), $\Gamma =1$ (middle) and $\Gamma =100$ (bottom). On the left, we show the decontamination front $y_*(t)$ and on the right the cleanser $c$ (red) and agent-layer thickness $R$ (blue) profiles. Numerical solutions of (4.6) are shown as solid lines with $\delta =0.25$ , and in the limit $\delta \ll 1$ (making the assumption $D=V=F=1$ ) as dotted lines, for comparison. Throughout we take $\rho _0=0.2$ . Solutions in (b), (d) and (f) are shown at the times $0.1t_D$ , $0.2t_D$ , $0.4t_D$ , $0.6t_D$ , $0.8t_D$ and $t_D$ , where $t_D$ is the decontamination time for the $\delta =0.25$ solution. The solutions for $\delta =0.25$ show clear discontinuities in $c_y$ at $y=1$ due to the discontinuous diffusivity $D(R)$ .

4.9. Decontamination times and decontamination efficiency

As in the sharp-interface case, the decontamination time and the decontamination efficiency are of practical interest. The decontamination time $t_D$ is the time when all the agent has been removed and occurs when $y_*(t_D)=1$ . We have found the analytic expressions (4.30) and (4.54) for $t_D$ in the limits $\Gamma \ll 1$ and $\Gamma \gg 1$ , respectively, in both cases assuming that $A\ll 1$ , $\delta \ll 1$ . The efficiency is defined (analogously to the sharp-interface case) as the ratio of the minimum amount of cleanser required to react away all the agent, to the total amount of cleanser that has entered the porous medium, namely

(4.64) \begin{align} \text{efficiency}=\frac{N}{\text{total flux in}}=\left (\frac{2\rho _0+\delta }{2\rho _0}\right )\frac{\Gamma }{-\int _{t=0}^{t_D}\left (DVc_y\right )|_{y=0}\,\mathrm{d}t}\,, \end{align}

using the expression (3.23) for $N$ , and noting that

(4.65) \begin{align} \text{total flux in}= - \frac{R_0D_0V_0}{\chi k L_{\text{AW}}}\int _{t=0}^{t_D}\left (DVc_y\right )|_{y=0}\,\mathrm{d}t \quad \text{mol}\,\text{m}^{-2}\,. \end{align}

In the limits in which we have approximate analytical solutions, we compute the approximate efficiency, finding that

(4.66) \begin{align} \text{efficiency}\sim \begin{cases}1 & \text{ if }A\ll 1, \delta \ll 1, \Gamma \ll 1,\text{ case Ia},\\[4pt] 1 & \text{ if }A\ll 1, \delta \ll 1, \Gamma \gg 1,\text{ case Ib},\\[11pt] \dfrac{2\sqrt{\pi \Gamma }}{\sqrt{A}} & \text{ if }A\gg 1, \delta \ll 1, \Gamma \ll 1,\text{case II}. \end{cases} \end{align}

In Figure 16, we explore the effect of the model parameters on the decontamination time. Broadly speaking, we observe similar dependence of the decontamination time on $A$ and $\Gamma$ as we saw for $\alpha$ and $\gamma$ , respectively, in the sharp-interface case discussed in Section 3.6. In Figure 16a, we show the variation of the numerically computed decontamination times as a function of $\Gamma$ , which we may interpret as the effect of the spill depth $L_{AW}$ on the decontamination time. We recall that for $A\ll 1$ (case I), for both $\Gamma \ll 1$ (case Ia) and $\Gamma \gg 1$ (case Ib) we have found the same expression for the decontamination time, namely $t_D=1+\Gamma/2$ . We see that, for small values of $A$ , this approximation $\Gamma/2+1$ is a good estimate for the decontamination time, even when $\Gamma =O(1)$ . For moderate values of $A=O(1)$ , we see that, while the approximation remains reasonable for $\Gamma \ll 1$ , the decontamination time is greater than $\Gamma/2+1$ for intermediate and large values of $\Gamma$ . For large values of $A$ , the small- $\Gamma$ approximation (case II) to $t_D$ (the solution of (4.61)) shows good agreement while $\Gamma$ is small. For $A\gt 1$ , the growth of $t_D$ with the spill depth appears to be roughly linear in $\Gamma$ , although we note that the scaling $A\Gamma$ suggested in Section 4.7.2 for large $A$ and $\Gamma$ gives an upper bound.

Figure 16. The behaviour of the decontamination time $t_D$ (in (a)), or the scaled decontamination time $t_D/A\Gamma$ (in (b)–(c)) or (4.70) (in (d)) for the agent-on-walls model, computed from the numerical solution of (4.6). In (a)–(c), we take the small- $\delta$ limit and investigate the effect of spill depth and reaction rate on the decontamination time via $\Gamma$ , and of the applied cleanser concentration via $A$ . In (d), we show the effect of the initial agent-layer thickness via $\delta$ .

As in the sharp-interface case, we note that the time scaling used in (4.5) depends on both the reaction rate $k$ and the applied cleanser concentration $c_0$ . In order to draw physical conclusions about the impact of these parameters on the decontamination time, in Figure 16b and c we therefore plot the scaled decontamination time

(4.67) \begin{align} \frac{t_D}{A\Gamma }=\frac{R_0D_0}{\chi k c_0f_2L_{\text{AW}}^2}t_D, \end{align}

which has the same dependence on $c_0$ and $k$ as the time scaling used in (4.5). In Figure 16b, we see that the decontamination time decreases with reaction rate $k$ , so that faster reaction rates result in faster decontamination. As in the sharp-interface case, we observe limiting returns. With $A\leq O(1)$ we see that for very large $k$ (and so large $\Gamma$ ), the decontamination time approaches a constant value, $1/2A$ . With $A\gt 1$ we also see much slower variation of $t_D$ with $\Gamma$ as $\Gamma$ becomes large. In Figure 16c, we show the variation of $t_D$ with $A$ , which is proportional to the applied cleanser concentration $c_0$ , and so we may interpret increasing $A$ with increasing the applied cleanser concentration. For stronger cleanser concentrations, we observe faster decontamination. The large- $A$ asymptotic solution given by (4.61) appears to be a good approximation at large $A$ values, even for large values of $\Gamma$ .

We now consider the effect of $\delta$ on the decontamination time. We note that for a fixed total volume of agent $v$ , as $\delta$ increases the depth $L_{\text{AW}}$ of the spill is decreased, since, from the definition (4.2) of $L_{\text{AW}}$ ,

(4.68) \begin{align} L_{\text{AW}}=\frac{v}{\pi \delta (2\rho _0+\delta )}. \end{align}

Accordingly, $\Gamma$ decreases as $\delta$ increases, since

(4.69) \begin{align} \Gamma =\frac{vk}{\pi \epsilon f_1D_0}\cdot \frac{1}{\delta (2\rho _0+\delta )}. \end{align}

In order to understand how the agent thickness impacts on the physical decontamination time while $v$ is held constant, we therefore compute solutions holding $\Gamma \delta (2\rho _0+\delta )$ constant. In order to understand the effect of $\delta$ on the dimensional decontamination time, we define another scaled decontamination time, namely

(4.70) \begin{align} \frac{t_D}{A\Gamma \delta ^2(2\rho _0+\delta )^3}=\frac{D_0 R_0\pi ^2 g}{2\rho _0\chi k c_0}t_D, \end{align}

which has the same dependence on $R_0$ , and so $\delta$ , as the time scaling used in (4.5). We plot this second scaled decontamination time (4.70) against $\delta$ in Figure 16d. The dominant effect of increasing $\delta$ is to significantly reduce the spill depth $L_{\text{AW}}$ , and thus we see that for larger $\delta$ and thicker agent layers on the pore walls, the scaled decontamination time decreases. We therefore conclude that the dimensional decontamination time decreases with the agent layer thickness $\delta$ . This behaviour is, of course, limited by the requirement that the cleanser-occupied region of the porespace remains connected, which is true so long as $\delta \lt 0.5-\rho _0$ .

In Figure 17, we explore the effect of the system parameters on the efficiency (4.64). Like the decontamination time, the efficiency also behaves mostly similarly with $\Gamma$ and $A$ as the sharp-interface efficiency behaved with $\gamma$ and $\alpha$ . In Figure 17a, we show the variation of the efficiency (4.64) with $\Gamma$ , observing that, for all values of $A$ , the efficiency increases with $\Gamma$ . Thus, for faster-reacting cleansers (larger values of $k$ and so $\Gamma$ ), less cleanser is wasted. Unlike in the sharp-interface case, for small $\Gamma$ we see the efficiency decreasing to zero, rather than approaching a constant as we saw in Figure 6d. This is because, in this agent-on-walls case, cleanser can diffuse through $y=1$ and be lost from the system, unlike in the sharp-interface case where cleanser was confined to $y\in [0,h(t)]$ . From (4.15), we find that the flux of cleanser through $y=1$ scales with $\sqrt{A\Gamma/t}$ . For small $\Gamma$ (when $t_D\rightarrow 1$ ), we therefore expect the cleanser lost through $y=1$ to be of order $\sqrt{A\Gamma }$ .

Figure 17. The behaviour of the decontamination efficiency (4.64) for the agent-on-walls model, computed from the numerical solution of (4.6). In (a)–(c), we take the small- $\delta$ limit. In (a), (b) and (d), we show the effect of $\Gamma$ , $A$ and $\delta$ , respectively, on the efficiency. In (c), we show the variation of the efficiency with the decontamination time, parameterised by $A$ .

From Figure 17b, we observe that, as in the sharp-interface case, the efficiency decreases with $A$ and so, for stronger applied cleanser concentrations (larger $c_0$ and so larger $A$ ), the efficiency worsens. Thus, as in the sharp-interface case, there is a pay-off between improving the decontamination time (by increasing the applied cleanser concentration $c_0$ ) and improving the efficiency (by decreasing $c_0$ ). In Figure 17c, we illustrate this pay-off by plotting the decontamination time against the efficiency, and we observe very similar qualitative behaviour to the sharp-interface case (in Figure 6f). In Figure 17d, we explore the effect of $\delta$ on the efficiency. As for Figure 16d, we hold the total initial volume of agent fixed, by holding $\Gamma \delta (2\rho _0+\delta )$ constant. We see that, as $\delta$ increases and $\Gamma$ is reduced, the efficiency drops.