The effect of pore-scale contaminant distribution on the reactive decontamination of porous media

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)