3.1. Steady states

The steady-state solutions to (1.1) satisfy

(3.1) \begin{align} \frac{d}{dx} \left (h^m \frac{d h}{d x} \right ) = h^n, \end{align}

an ODE that will subsequently arise as a quasi-steady balance in (1.1) for circumstances in which the time derivative is non-zero but negligible. On integrating with respect to $x$ , we have that

(3.2) \begin{align} \frac{1}{2} h^{2m} \left ( \frac{dh}{dx} \right )^2 = \frac{1}{m+n+1} h^{m+n+1} + C_1, \end{align}

for arbitrary constant $C_1$ . Given the scaling invariance of (3.1),

(3.3) \begin{align} x \mapsto \phi^{m-n+1} x, \quad h \mapsto \phi ^2 h, \end{align}

where $\phi$ is an arbitrary constant, we may scale $|C_1|$ in (3.2) to a convenient value. In this section, we shall set $C_1=0$ , $C_1=1/2$ and $C_1=-1/(m+n+1)$ in turn. Given (3.3), these zero, positive and negative values in effect cover the full range of possibilities, the specific values being chosen for algebraic convenience in the sense of (3.7) and (3.9) below.

Figure 3. Solutions to (3.1) for $m=3$ and $n=0$ . The black curve corresponds to the exceptional solution (1.9); the yellow curve (the leaking solution) satisfies (3.6) while the green curve (the pre-touchdown) satisfies (3.8).

The exceptional solution

For $C_1=0$ , we have the steady solution (1.9), up to translational invariance. In Sections 3.3, 5 and 6, we shall see that it is a common component of a host of relevant solutions; owing to its significance in the remainder of the analysis, we shall refer to it as the ‘exceptional solution’. Linearising (3.2) about (1.9), we obtain

(3.4) \begin{align} h \sim A_* x^{\alpha } + \nu _1 x^{\alpha -1} + \nu _2 x^{-(m+n)\alpha } \end{align}

where

(3.5) \begin{align} \alpha =\frac{2}{m+1-n}. \end{align}

The $\nu _1$ term in (3.4) reflects translational invariance in $x$ and is less singular as $x \rightarrow 0$ than the $\nu _2$ term when $m+3n+1\gt 0$ , a condition that we impose throughout.

The leaking solution

Setting $C_1=1/2$ in (3.2) and imposing the interface condition (1.5) gives

(3.6) \begin{align} \frac{1}{2} h^{2m} \left ( \frac{dh}{dx} \right )^2 = \frac{1}{m+n+1} h^{m+n+1} + \frac{1}{2}, \quad h|_{x=0}=0. \end{align}

A numerical solution is shown in yellow in Figure 3. This solution, which we shall refer to as the ‘leaking solution’, has a unit outward flux through the origin, that is, it has the property

(3.7) \begin{align} \left. h^m \frac{dh}{dx} \right |_{x=0} = 1, \end{align}

and is a candidate to form part of the asymptotic description of a solution to (1.1) that is attached to an absorbing boundary. The solution to (3.6) has asymptotic behaviour given by (3.4) as $x \to +\infty$ , with $\nu _1$ and $\nu _2$ depending on $m$ and $n$ , the $\nu _1$ term being the dominant correction term for $m+3n+1\gt 0$ . For $n=0$ , we have the explicit solution

\begin{align*} \frac{1}{m+1} h^{m+1} =\frac{1}{2} x^2 +x. \end{align*}

The pre-touchdown solution

Finally, setting $C_1=-1/(m+n+1)$ in (3.2) and imposing a boundary condition that requires the solution to have zero gradient at the origin leads to

(3.8) \begin{align} \frac{1}{2} h^{2m} \left ( \frac{dh}{dx} \right )^2 = \frac{1}{m+n+1}\left ( h^{m+n+1} -1 \right ), \quad \left. \frac{dh}{dx} \right |_{x=0} = 0. \end{align}

The solution to this problem (determined numerically) is shown in green in Figure 3. We shall refer to this solution as the ‘pre-touchdown solution’, and it has the property

(3.9) \begin{align} h|_{x=0} = 1. \end{align}

This solution, scaled by an appropriate function of time (i.e. in quasi-steady form), will play a role in the pre-touchdown dynamics. Once again, the solution (3.8) has the far-field behaviour (3.4). In the case $n=0$ , it has explicit solution

\begin{align*} \frac{1}{m+1} h^{m+1} =\frac{1}{2} x^2 +\frac{1}{m+1}. \end{align*}

Exceptionally, the far-field behaviour of (3.8) has $\nu _1=0$ in (3.4) when $n=0$ . For $n\lt 0$ , we have $\nu _1 \lt 0$ in (3.4) while for $n\gt 0$ we have $\nu _1 \gt 0$ ; this has implications for the intermediate-asymptotic behaviour. More precisely,

(3.10) \begin{align} h^{\frac{m+1-n}{2}} \sim \frac{m+1-n}{(2(m+n+1))^{1/2}} x + \frac{n}{m+n+1} B \left ( \frac{1}{2}, \frac{m+3n+1}{2(m+n+1)} \right ) \quad \text{as} \quad x \to + \infty, \end{align}

where $B$ denotes the beta function and the constraint $m+3n+1\gt 0$ again rears its head.

3.2. Travelling-wave solutions

Here we consider travelling-wave solutions to (1.1),

(3.11) \begin{align} h=h(\zeta ), \qquad \zeta =x- v t, \end{align}

where the constant $v \neq 0$ is the wave speed. Substitution of (3.11) into (1.1) yields the ODE

(3.12) \begin{align} -v\frac{ d h}{d \zeta }=\frac{ d }{d \zeta }\left ({h}^m \frac{ d h}{d \zeta } \right )-{h}^n. \end{align}

We shall require a solution that satisfies (1.3) at $\zeta =0$ , so that

(3.13) \begin{align} h|_{\zeta =0}=0, \qquad \left. h^m \frac{\partial h}{\partial \zeta } \right |_{\zeta =0}=0. \end{align}

Subsequently, for solutions with a moving interface, $v$ will be identified with $\dot{s}$ , with the travelling-wave balance (3.12), again in quasi-steady form, providing the local behaviour about the interface whenever $\dot{s}\neq 0$ .

Via the rescalings

(3.14) \begin{align} h \to |v|^{\frac{2}{m+n-1}} h, \qquad \zeta \to |v|^{\frac{m+1-n}{m+n-1}} \zeta \end{align}

we can without loss of generality set $v=\pm 1$ in (3.12). The phase plane for an advancing interface ( $v\lt 0$ ) is plotted in Figure 4 panel (a). The initial conditions (3.13) correspond to the sole trajectory that emanates from the origin, the asymptotic behaviour of which is given in terms of the unscaled variables by

(3.15) \begin{align} h \sim \left (({-}v) m \zeta \right )^{1/m} \quad \mbox{as} \quad \zeta \to 0^+ \quad \text{with} \quad v\lt 0, \end{align}

corresponding to the local behaviour (2.3). The phase plane for $v\gt 0$ (a receding interface) is plotted in Figure 4 panel (b). The conditions (3.13) are associated with the sole trajectory that emanates from the origin, the unscaled version of which satisfies

(3.16) \begin{align} h \sim \left (\frac{(1-n)}{v} \zeta \right )^{1/(1-n)} \quad \mbox{as} \quad{\zeta \rightarrow 0^+} \quad \text{with} \quad v\gt 0. \end{align}

coinciding with (2.2).

Figure 4. Phase planes for (a) the advancing ( $v\lt 0$ ) and (b) the receding ( $v\gt 0$ ) travelling waves. Here $m=3$ and $n=0$ , but qualitatively similar results apply for other exponents in the range (1.2). The red curves indicate the trajectories leaving the origin and corresponding to the solution which links the behaviours (3.15) and (3.17) in panel (a), and (3.16) and (3.17) in panel (b). Since $v \to -v$ results from $\zeta \to -\zeta$ in (3.12), (a) and (b) can be viewed as upper and lower quadrants of the same phase plane on also swapping the direction of the trajectories in one case. All the trajectories have the same far-field behaviour, (3.18).

Far-field behaviour

The phase-plane analysis demonstrates that both advancing and receding interfaces have the same far-field behaviour, with the left-hand side of (3.12) being negligible at leading order, so that

(3.17) \begin{align} h \sim A_* \zeta^{\alpha } \quad \mbox{as} \quad \zeta \rightarrow +\infty. \end{align}

We shall need a more detailed description of the behaviour; however, the relevant correction terms follow from (3.4) and from the left-hand side of (3.12) and result in the far-field behaviour

(3.18) \begin{align} h \sim A_* \zeta ^\alpha - v B \zeta^{\alpha -\beta } + \nu _1\zeta^{\alpha -1} + \nu _2 \zeta^{-(m+n)\alpha } \quad \text{as} \quad \zeta \to +\infty, \end{align}

where the constants $\beta$ and $B$ are given by

(3.19) \begin{align} \beta =\frac{m+n-1}{m+1-n}, \quad B= \frac{(2(m+n+1))^{\frac{m-1}{m+1-n}}(m+1-n)^{\frac{3-m-n}{m+1-n}}}{(1-n) (n+m+3)}. \end{align}

The $B$ term is the dominant correction term in (3.18), given the condition (1.2), while $\nu _1$ and $\nu _2$ are the requisite two degrees of freedom. The phase planes confirm that the IVP (3.12)–(3.13) has a unique solution for given $v$ , $m$ and $n$ ; the dependencies of $\nu _1$ and $\nu _2$ upon $v$ can be inferred from the above noted scaling property, their dependence on $m$ and $n$ depending on the sign of $v$ .

In the light of these properties, and in the parameter range (1.2), the velocity of the travelling wave can thus be determined if the coefficient of the third term in far-field behaviour (3.18) is known, this being the dominant correction term in (3.18), since $\beta \gt 1$ .

3.3. Self-similar solutions

As we shall see below, the backward self-similar reduction for (1.1),

(3.20) \begin{align} h(x,t) = \left ({-}t \right )^{\frac{1}{1-n}} \; f(\xi ), \quad \xi = x ({-}t)^{-\frac{m+1-n}{2(1-n)}} \quad \mbox{for} \quad t \lt 0, \end{align}

coupled to the free and absorbing boundary conditions (1.3) and (1.5), provides a plausible scenario for times leading up to a singular event. As in [Reference Foster, Please, Fitt and Richardson7], substituting (3.20) into (1.1) reveals that $f$ satisfies the ODE

(3.21) \begin{align} \frac{1}{1-n} \left ({-} f + \frac{m+1-n}{2} \, \xi \frac{d f}{d \xi } \right ) = \frac{d}{d \xi } \left ( f^m \frac{d f}{d \xi } \right ) - f^n. \end{align}

The power-law solution to (3.21) is the exceptional solution (1.9), which in similarity variables reads

(3.22) \begin{align} f = A_* \xi^{\alpha }. \end{align}

As a precursor to constructing solutions to (3.21) which give rise to dynamic interfaces, it is helpful to examine the possible near-field behaviours of $f$ about the interface at $\xi = \hat{\xi }$ . We find that

(3.23) \begin{align} f \sim \left ( \frac{m (m+1-n) \hat{\xi }}{2 (1-n)} (\xi - \hat{\xi }) \right )^{1/m} \quad \textrm{as} \quad \xi \to \hat{\xi }^+ \quad \mbox{for} \quad \hat{\xi } \gt 0, \end{align}

(3.24) \begin{align} f \sim \left ( \frac{2 (1-n)^2}{(m+1-n)({-}\hat{\xi })} (\xi - \hat{\xi }) \right )^{1/(1-n)} \quad \textrm{as} \quad \xi \to \hat{\xi }^+ \quad \mbox{for} \quad \hat{\xi } \lt 0 \end{align}

are admissible. Here, $\hat{\xi } \neq 0$ is a free parameter. In the context of the original PDE, (1.1), the former behaviour corresponds to an advancing interface with $\dot{s}\lt 0$ (as, e.g., depicted by the solid curves in Figure 2(i), which show an advancing solution before reversing), while the latter corresponds to a receding interface with $\dot{s}\gt 0$ . The behaviours (3.23)–(3.24) correspond to (2.2)–(2.3). In [Reference Foster, Gysbers, King and Pelinovsky5–Reference Foster, Please, Fitt and Richardson7], the Liouville-Green method was used (which involves linearising about (3.23) in (3.21) and examining the self-consistent asymptotic behaviours) to demonstrate that there is only one degree of freedom in both of these behaviours, namely $\hat{\xi }$ .

Even solutions satisfying

(3.25) \begin{align} f \sim \hat{f} + \frac{1}{2} \hat{f}^{-m} \left (\hat{f}^n-\frac{1}{1-n} \hat{f} \right ) \xi ^2 \quad \mbox{as} \quad \xi \to 0, \end{align}

wherein $\hat{f}\gt 0$ is the sole free parameter, are also of potential relevance. This behaviour clearly does not correspond to an interface of (1.1) but, in time-dependent setting, it corresponds to a point of no flux and is a candidate to describe a solution leading up to a touchdown event.

A fourth alternative is

(3.26) \begin{align} f \sim K \xi^{1/(m+1)} \quad \mbox{as} \quad \xi \to 0^+, \end{align}

where $K\gt 0$ is the sole free parameter. This leads to a static interface of (1.1) through which mass is being lost (as, e.g., depicted by the solid ones in 2(iv)). There is a final near-field behaviour of relevance, namely the exceptional solution (1.9); in view of (3.4), this contains no degrees of freedom.

The only relevant far-field behaviour for solutions to (3.21) is

(3.27) \begin{equation} f \sim A \xi^{\alpha } \quad \mbox{as} \quad \xi \to +\infty, \end{equation}

and it has been shown in [Reference Foster, Please, Fitt and Richardson7] that the parameter $A\gt 0$ is the only degree of freedom in the large $\xi$ behaviour. Unpicking the transformation to the self-similar variables, (3.27) corresponds to

(3.28) \begin{align} h \sim A x^\alpha \end{align}

and is a candidate for the local behaviour as $x \to 0^+$ of the solution at $t=0$ , given that $\xi \to + \infty$ as $t \to 0^-$ at fixed $x$ . We note that (3.27) coincides with the exceptional solution (1.9) in the case when $A=A_*$ . We shall return in Section 6 to a detailed exploration of these self-similar forms.

4.1. Porous-medium equation balance

Neglecting the sink term in (1.1) gives the so-called porous-medium equation, that is,.

(4.1) \begin{align} \frac{\partial h}{\partial t} \sim \frac{\partial }{\partial x} \left (h^m \frac{\partial h}{\partial x} \right ) \end{align}

the similarity solutions of which (as well as its other properties) have been very widely investigated, see, for example, [Reference Lacey, Ockendon and Tayler21]. We potentially need both the backward ( $t\lt 0$ ) self-similar solutions

(4.2) \begin{align} h \sim ({-}t)^\nu f_-(\eta _-), \quad \eta _- = x ({-}t)^{-\frac{m\nu +1}{2}} \quad \text{as} \quad t \to 0^- \quad \text{with} \quad \eta _-=O(1) \end{align}

and the forward ( $t\gt 0$ ) ones

(4.3) \begin{align} h \sim t^\nu f_+(\eta _+), \quad \eta _+ = x t^{-\frac{m \nu +1}{2}} \quad \text{as} \quad t \to 0^+ \quad \text{with} \quad \eta _+=O(1). \end{align}

The case $\nu =1/m$ is of most significance to us for reasons that will become clear later: in this case

\begin{align*} f_-^m = m\hat{\eta }( \eta _- -\hat{\eta })_+, \qquad f_+^m = m\hat{\eta }( \eta _- + \hat{\eta })_+ \end{align*}

for some constant $\hat{\eta }$ ; these are each both a scaling reduction and a travelling-wave solution to (4.1) and represent an exceptional connection in the phase space of $f_-$ – see [Reference Hulshof, King and Bowen14]. Neglect of the $h^n$ terms as $|t| \to 0$ requires, for legitimacy, that

(4.4) \begin{align} \nu \lt 1/(1-n), \end{align}

which in the case $\nu =1/m$ implies that $m+n\gt 1$ , consistent with (1.2).

4.2. Absorption balance

The last of the three possible balances in which one of the terms in (1.1) is entirely neglected is

(4.5) \begin{align} \frac{\partial h}{\partial t} \sim - h^n \end{align}

so that

(4.6) \begin{align} h^{1-n} \sim h_0^{1-n}(x) + (1-n) ({-}t). \end{align}

Here, $h_0(x)$ is an arbitrary function. As we shall see, the balance (4.5) can play an asymptotic role for power law

(4.7) \begin{align} h_0(x) = M |x|^{\mu/(1-n)} \end{align}

for positive constants $M$ and $\mu$ . In such instances, it is instructive to interpret (4.6) as backward and forward similarity reductions,

(4.8) \begin{align} h = ({-}t)^{\frac{1}{1-n}} g_-(\omega _-), \quad \omega _- = x ({-}t)^{-1/\mu }, \quad g_-^{1-n} = M^{1-n} |\omega _-|^{\mu } + (1-n), \end{align}

(4.9) \begin{align} h = t^{\frac{1}{1-n}} g_+(\omega _+), \quad \omega _+ = x t^{-1/\mu }, \quad g_+^{1-n} = M^{1-n} |\omega _+|^{\mu } - (1-n), \end{align}

from which it is clear that self-consistency as $|t| \to 0$ in the neglect of the diffusion term in (1.1) requires

(4.10) \begin{align} \mu m + (\mu -2) (1-n)\gt 0. \end{align}

5.1. Linearisation about the uniform state

We linearise about the uniform state, (1.10), by setting

(5.1) \begin{align} h \sim \left ( (1-n)({-}t) \right )^{1/(1-n)} +{\mathcal H}(x,t) \end{align}

in (1.1) and neglect nonlinear terms in $\mathcal H$ to give

(5.2) \begin{align} \frac{\partial{\mathcal H}}{\partial t} = \left ( (1-n)({-}t) \right )^{m/(1-n)} \frac{\partial ^2{\mathcal H}}{\partial x^2} - \frac{n}{(1-n)({-}t)}{\mathcal H}. \end{align}

Defining

(5.3) \begin{align}{\mathcal H} = ({-}t)^{\frac{n}{1-n}}{\mathcal G}(x,\tau ), \quad -\tau = \frac{((1-n)({-}t))^{\frac{m+1-n}{1-n}}}{m+1-n} \end{align}

transforms (5.2) to the heat equation

(5.4) \begin{align} \frac{\partial{\mathcal G}}{\partial \tau } = \frac{\partial ^2{\mathcal G}}{\partial x^2}. \end{align}

The backward self-similar solutions without exponential growth as $|x| \to +\infty$ are given by the Hermite polynomials $H_k$ , namely

(5.5) \begin{align}{\mathcal G} = C_k ({-}\tau )^{k/2} H_k \left ( x ({-}\tau )^{-1/2}/ 2 \right ) \end{align}

where $k$ is a non-negative integer and $C_k$ is some constant.Footnote ²

The following properties will be needed in due course. Firstly, (5.3)–(5.5) require, for the self-consistency of (5.1) (i.e. that the second term therein be negligible in comparison to the first in the limit that $t \rightarrow 0^-$ ), that

(5.6) \begin{align} km + (k-2)(1-n)\gt 0; \end{align}

the case $k=0$ is therefore always excluded, unsurprisingly since it simply corresponds to a translation in $t$ . Secondly, a standard result for $H_k$ implies that (5.5) renders the matching condition

(5.7) \begin{align}{\mathcal H}(x,t) \sim C_k ({-}t)^{-\frac{n}{1-n}} x^k \quad \text{as} \quad | x ({-}\tau )^{-1/2} | \to \infty, \end{align}

so that

(5.8) \begin{align} h^{1-n} \sim (1-n)({-}t) + (1-n)^{\frac{1-2n}{1-n}} C_k x^k, \end{align}

cf. (4.6), with (5.6) then corresponding to (4.10).

5.2. Linearisation about the exceptional solution

The other linearised problem that we shall need to characterise is that about the exceptional solution, (1.9). On setting

(5.9) \begin{align} h \sim A_* x^{\alpha } + \bar{\mathcal H}(x,t) \end{align}

in (1.1), linearising in $\bar{\mathcal H}$ gives

(5.10) \begin{align} \frac{\partial \bar{\mathcal H}}{\partial t} = A_*^m \frac{\partial ^2}{\partial x^2} \left ( x^{m \alpha } \bar{\mathcal H} \right )- n A_*^{-(1-n)} x^{-(1-n) \alpha } \bar{\mathcal H}. \end{align}

Elucidating the behaviour of solutions to (5.10) will provide information crucial to what follows, specifically with regard to backward self-similarity. In order to analyse the linear problem (5.10), it is helpful to use the definitions of $A_*$ and $\alpha$ to rewrite it in the form

(5.11) \begin{align} x^{2-m \alpha } \frac{\partial \bar{\mathcal H}}{\partial t} = A_*^m \left ( x^{2-m \alpha } \frac{\partial ^2}{\partial x^2} \left ( x^{m \alpha } \bar{\mathcal H} \right )- (\alpha (m+1)-2)(\alpha (m+1)-1) \bar{\mathcal H} \right ). \end{align}

The right-hand side of (5.11) dominates as $x \rightarrow 0$ , since $2-m\alpha \gt 0$ in the range of parameters that we consider, (1.2). The possible solution behaviours are thus found by neglecting the left-hand side of (5.11) and solving the resulting Euler equation in $x$ for $\bar{\mathcal H}$ ; these are

(5.12) \begin{align} \bar{\mathcal H} \sim \nu _1(t) x^{\alpha -1}, \qquad \bar{\mathcal H} \sim \nu _2(t) x^{-(m+n)\alpha }, \end{align}

where $\nu _1(t)$ and $\nu _2(t)$ are arbitrary functions, as in (3.4). There are two distinct cases to be considered: for $m+3n +1\gt 0$ , the second of (5.12) is the more singular so it is natural to specify

(5.13) \begin{align} \bar{\mathcal H}=o\left ( x^{-(2m+1)\alpha +2} \right ) \quad \mbox{as} \quad x \to 0^+ \quad \mbox{if} \quad m+3n +1\gt 0 \end{align}

as a boundary condition on (5.10), the local behaviour then being given by the first of (5.12), corresponding to a (small) translation of $x$ in (5.9). By contrast, if $m+3n +1\lt 0$ the first of (5.12) is the more singular so it is natural instead to specify

(5.14) \begin{align} \bar{\mathcal H}=o\left ( x^{\alpha -1}\right ) \quad \mbox{as} \quad x \to 0^+ \quad \mbox{if} \quad m+3n +1\lt 0. \end{align}

As already noted, henceforth we only consider the case $m+3n +1\gt 0$ , so that (5.13) applies (though in (II) below we record an instructive result for $m+3n+1\lt 0$ ); the analysis of the converse case proceeds on similar lines, but will be omitted here (in part because it appears of less physical relevance).

Two transformations, (I) and (II) below, reduce (5.10) to standard forms. In both cases, we define

(5.15) \begin{align} r = x^{\frac{2-m \alpha }{2}}, \end{align}

the definitions of $G_i$ being guided by (5.12). We have

(5.16) \begin{align} \text{(I): } \quad \bar{\mathcal H}(x,t) = x^{\alpha -1} G_1(r,t), \quad m+3n+1\gt 0, \end{align}

(5.17) \begin{align} \frac{\partial G_1}{\partial t} = \frac{1}{4} (2-m \alpha )^2 A_*^m \left ( \frac{\partial ^2 G_1}{\partial r^2} + \frac{N_1 - 1}{r} \frac{\partial G_1}{\partial r} \right ), \quad N_1 = 2\frac{(m+2)\alpha - 1}{2-m \alpha }. \end{align}

(5.18) \begin{align} \text{(II): } \quad \bar{\mathcal H}(x,t) = x^{-(m+n)\alpha } G_2(r,t), \quad m+3n+1\lt 0, \end{align}

(5.19) \begin{align} \frac{\partial G_2}{\partial t} = \frac{1}{4} (2-m\alpha )^2 A_*^m \left ( \frac{\partial G_2}{\partial r^2} + \frac{N_2 - 1}{r} \frac{\partial G_2}{\partial r} \right ), \quad N_2= 2\frac{5-(3m+2)\alpha }{2-m \alpha }. \end{align}

These representations allow application of standard theory for the radially symmetric heat equation.Footnote ³ As was pursued in [Reference Foster, Gysbers, King and Pelinovsky5], the standard forms for the backward similarity solutions to (5.17) and (5.19) are

(5.20) \begin{align} G = ({-}t)^M g(y), \quad y = \frac{2 r}{ A_*^{m/2} (2-m \alpha ) ({-}t)^{1/2}} \end{align}

(where $M$ is an arbitrary constant) and satisfy

(5.21) \begin{align} z \frac{d^2 g}{d z^2} + \left ( \frac{N}{2} - z \right ) \frac{dg}{dz} + Mg=0, \quad z = y^2, \end{align}

so that the solutions having the necessary regularity as $z \to 0$ (namely $g$ finite as $z \to 0$ ) and growing algebraically, rather than exponentially (of the form $\log g \sim z$ as $z \to +\infty$ ), as $z \to +\infty$ are given by the generalised Laguerre polynomials,

(5.22) \begin{align} g \propto L_P^{(N-2)/2} (z), \end{align}

where $P = M + (N-2)/2$ ; since $(N-2)/2$ need not be an integer, this is not the standard case, but $M$ is required to be a positive integer, leading to $g$ being an $M$ -th order polynomial in $z$ (i.e. in (5.20)–(5.21), $M$ should be viewed as an eigenvalue, leading to the requirement that it be a non-negative integer). The first few eigenmodes (radially symmetric heat polynomials, representing a complete set of eigenmodes) are thus multiples of the following:

(5.23) \begin{align} M = 0, \quad g_0 = 1, \end{align}

(5.24) \begin{align} M = 1, \quad g_1 = z - N/2, \end{align}

(5.25) \begin{align} M = 2, \quad g_2 = z^2-(N+2)z+N(N+2)/4, \end{align}

(5.26) \begin{align} M=3, \quad g_3 = z^3 - 3(N+4)z^2/2 + 3 (N+2)(N+4)z/4 - N(N+2)(N+4)/8. \end{align}

We shall subsequently need to view the leading term in (5.9) as a similarity solution of the form (4.9), while in case (I) above (5.20) implies that $\bar{\mathcal H}$ is of the self-similar form

(5.27) \begin{align} \bar{\mathcal H} = ({-}t)^L F(\xi ), \quad \text{where} \quad L = M+\frac{1+n-m}{2(1-n)}, \end{align}

so self-consistency of (5.9) requires that

(5.28) \begin{align} M \gt \frac{m+1-n}{2(1-n)}. \end{align}

Hence, $M=0$ and $M=1$ are excluded (in the latter case (5.28) would require $m+n\lt 1$ for $n\lt 1$ , again identifying familiar borderline cases). This is to be expected in view of the invariance of (1.1) under translations of $x$ and $t$ , illustrating a more general principle (we touch on some of the broader implications in Section 9). Equation (5.10) inherits $t$ , but not $x$ , translation invariance, and we pursue the reasoning in two stages. It follows from $t$ translation invariance that the $M$ th mode is proportional to the time derivative of the $(M+1)$ th mode, so that

(5.29) \begin{align} g_M \propto z \frac{d g_{M+1}}{dz} - (M+1) g_{M+1}. \end{align}

Hence, we have

(5.30) \begin{align} \bar{\mathcal H}(x,t+t_0) = \sum _{k=0}^M a_{k,M} t_0^k \bar{\mathcal H}_{M-k}(x,t) \end{align}

(where $\bar{\mathcal H}_p$ is the solution to (5.10) associated with $g_p$ ) with

(5.31) \begin{align} a_{0,M} = 1, \quad a_{k,M} = \prod _{L=0}^{k-1} \alpha _{M-L}/ k!, \quad k \geq 1, \end{align}

where

(5.32) \begin{align} \bar{\mathcal H}_{M-1} (x,t) = \alpha _M \frac{\partial }{\partial t} \bar{\mathcal H}_M (x,t). \end{align}

Writing the solution to (5.10) in the form

(5.33) \begin{align} \bar{\mathcal H} \sim \sum _{p=0}^\infty K_p \bar{\mathcal H}_p(x,t), \end{align}

if we now choose $t_0$ in (5.30) such that

(5.34) \begin{align} K_1 = \sum _{M=1}^\infty a_{M-1,M} t_0^{M-1} \end{align}

then the time-translated version of $\mathcal H$ has $K_1=0$ . This translation will modify the value of $K_0$ , but here $x$ translation invariance comes into play: since $\bar{\mathcal H}_0 \propto d x^\alpha/ dx$ , translating suitably in (5.9) allows $K_0$ also to be set to zero.

The condition (5.28) will play a crucial role in what follows, the borderline cases

(5.35) \begin{align} m=(2M-1)(1-n), \qquad M=2,3, \cdots \end{align}

identifying the bifurcation points in Figure 7 below.

Figure 5. Results of the shooting scheme to detect viable solutions (3.21) for $m=2.5$ and $n=0$ . The upper panel, (a), shows the variation of $f$ , $f^m df/d\xi$ and $\xi _{\text{term}}$ at the termination point ( $\xi =\xi _{\text{term}}$ ) with the shooting parameter $A$ . Candidate receding, exceptional and attached trajectories are indicated by the blue, black and yellow arrows, respectively. The lower panels, (b) and (c), show the candidate receding and attached trajectories (solid curves), respectively, along with a few trajectories ‘near’ to these candidates (dashed curves) – “near” in the sense that they have values of $A$ close to the candidates.

6.1. Backward self-similarity

We now analyse in detail some of the phase space of the ODE (3.21) associated with the full backward similarity reduction (1.7) of the PDE (1.1). In order to construct solutions to the ODE (3.21) that give rise to viable dynamic solutions to (1.1), it remains to determine whether it is possible to connect one of the near-field behaviours (3.23)–(3.26) with the far-field behaviour (3.27). We shall address this question by setting up an IVP (i.e. a shooting problem) from $\xi = \infty$ . It was proved in Lemma 4.4 of [Reference Foster and Pelinovsky6] that trajectories emanating from (3.27) with $A\gt 0$ will, when integrated in the direction of decreasing $\xi$ , either reach $f=0$ or $f^m df/d\xi =0$ at some finite $\xi$ . In lieu of an analytical means to solve the connection problem, the result from [Reference Foster and Pelinovsky6] motivates its numerical study by means of a shooting scheme in which $A\gt 0$ is varied as the shooting parameter. Once a value for $A$ has been adopted, the behaviour (3.27) can be used to set initial conditions to integrate (3.21) backwards in $\xi$ , and this process can be terminated at the value $\xi =\xi _{\text{term}}$ at which $f=0$ or $f^m df/d\xi =0$ (whichever occurs first). Once this has been done, it is then possible to assess whether one of the behaviours (3.23)–(3.26) has been reached by reading off the values of $f|_{\xi =\xi _{\text{term}}}$ , $f^m df/d\xi |_{\xi =\xi _{\text{term}}}$ and $\xi _{\text{term}}$ . If the integration stopped at a point where

1. (a) $f|_{\xi =\xi _{\text{term}}}=0$ , $\displaystyle f^m \left. \frac{df}{d\xi } \right |_{\xi =\xi _{\text{term}}}=0$ with $\xi _{\text{term}}\gt 0$ , the behaviour (3.23) has been attained,

2. (b) $f|_{\xi =\xi _{\text{term}}}=0$ , $\displaystyle f^m \left. \frac{df}{d\xi } \right |_{\xi =\xi _{\text{term}}}=0$ with $\xi _{\text{term}}\lt 0$ , the behaviour (3.24) has been attained,

3. (c) $f|_{\xi =\xi _{\text{term}}}\gt 0$ , $\displaystyle f^m \left. \frac{df}{d\xi } \right |_{\xi =\xi _{\text{term}}}=0$ with $\xi _{\text{term}}=0$ , the behaviour (3.25) has been attained,

4. (d) $f|_{\xi =\xi _{\text{term}}}=0$ , $\displaystyle f^m \left. \frac{df}{d\xi } \right |_{\xi =\xi _{\text{term}}} \gt 0$ with $\xi _{\text{term}}=0$ , the behaviour (3.26) has been attained.

The additional constraints in (a)–(d), namely a pair of boundary conditions in (a) and (b) and the requirement that $\xi _{\text{term}}$ be zero in (c) and (d), lead to the selection of a discrete set of admissible values of $A$ for given $m$ and $n$ , as exemplified in Figure 7. The results of this shooting scheme are summarised in Figures 5 and 6 for $m=2.5$ , $n=0$ and $m=4.5$ , $n=0$ , respectively. Similar bifurcation diagrams have been presented in the authors’ previous work, [Reference Foster, Gysbers, King and Pelinovsky5, Reference Foster and Pelinovsky6], but this is the first time that solutions reaching states (c) or (d) have been found. In the former case, one viable solution with a receding interface (blue), satisfying (3.24), was identified, as well as one with an attached interface (yellow), satisfying (3.26). In the latter case, two viable solutions were found, one with the local behaviour (3.26), which is a potential profile prior to touchdown (green), and another with an advancing interface (red) and the near-field behaviour (3.23). Similar experiments were carried out for many values of $m \in (1,8)$ while retaining $n=0$ and the locations of the viable solutions are summarised in the top right panel of Figure 7. Further computations were carried out for other values of $n$ , namely $-1,-1/2$ and $1/2$ , and similar results were found, see the other panels in Figure 7. In summary, we observe that for a particular pair of values of $m$ and $n$ it is possible to find connections between the far-field behaviour (3.27) and the near-field behaviours (3.23)–(3.26) for isolated values of $A$ .

Figure 6. Results of the shooting scheme to detect viable solutions (3.21) for $m=4.5$ and $n=0$ . The upper panel, (a), shows the variation of $f$ , $f^m df/d\xi$ and $\xi _{\text{term}}$ at the termination point ( $\xi =\xi _{\text{term}}$ ) with the shooting parameter $A_*$ . Candidate exceptional, touchdown and advancing trajectories are indicated by the black, green and red arrows, respectively. The lower panels, (b) and (c), show the candidate touchdown and attached trajectories (solid curves), respectively, along with a few trajectories ‘near’ to these candidates (dashed curves).

Figure 7. Locations of viable solutions to the self-similar ODE for (3.21). The exceptional trajectory is indicated by the presence of a black dot. Advancing and receding interface solutions satisfying (3.23) and (3.24) are indicated by red and blue dots, respectively. A touchdown candidate with the near-field behaviour (3.25) is indicated by a green dot, and a leaking candidate, satisfying (3.26), is indicated by a yellow dot.

6.2. Forward self-similarity

When a backward similarity solution exists, as described in Section 6.1, it is a candidate for describing the local behaviour as the relevant type of singularity is approached. The immediate post-singularity behaviour is then expected to be a forward similarity solution of the form (1.8) with

(6.1) \begin{align} f(\xi ) \sim A \xi^{\alpha } \qquad \mbox{as} \quad \xi \rightarrow +\infty, \end{align}

where the admissible values of $A$ are determined by the pre-singularity backward solution; (6.1) is supplemented by the relevant boundary conditions at the interface. Importantly, the forward problem is not of eigenvalue type: in contrast to (1.7), an exponentially decaying correction to (6.1) is available and can be constructed by the Liouville-Green method, so the boundary value problem is correctly specified for given $A$ , rather than the possible $A$ being determined as part of the solution (the corresponding quantity then being exponentially growing rather than decaying). It was proved in Corollary 4.2 of [Reference Foster and Pelinovsky6] that reversing and anti-reversing solutions emanating from suitable behaviour at finite $\xi$ all reach (6.1), for some value of the parameter $A$ , when integrated in the direction of increasing $\xi$ . Associated numerical results are shown in Figure 8, for $n=0$ , $m=2,3,4$ , the ODE being solved in the direction of increasing $\xi$ , the approach confirming that arbitrary $A\gt 0$ can be attained in this fashion.

Figure 8. The connection problem for the self-similar ODE for $t\gt 0$ ; illustrating that connections can be made to arbitrary positive values of $A$ (viable reversers have $A \lt A_*$ and viable anti-reversers have $A \gt A_*$ ).

7.1. Preliminaries

We shall exploit the $x$ and $t$ translation invariance of (1.1) to locate the singularity at $x=0$ and $t=0$ , without loss of generality (in particular, where a fixed boundary is involved, this will be taken to be at $x=0$ ). These symmetries will play further important roles in what follows.

Section 5 gives a detailed description of the two relevant linearisations in terms of special functions. Here we summarise the results in more informal terms and in a self-contained fashion to highlight the information we shall need in the rest of the section. Firstly, the linearisation (5.1) about the spatially uniform solution

(7.1) \begin{align} h=((1-n)({-}t))^{\frac{1}{1-n}}, \end{align}

leads to the linear PDE (5.2). The relevant backward similarity solutions to this PDE take the form

(7.2) \begin{align}{\mathcal H}_M= ({-}t)^{\frac{n}{1-n}}({-}t)^{\frac{(m+1-n)M}{2(1-n)}} \Psi _M(\eta ), \quad \eta =\sqrt{\frac{m+1-n}{2}} (1-n)^{-\frac{(m+1-n)}{2(1-n)}} \frac{x }{({-}t)^{\frac{(m+1-n)}{2(1-n)}}}, \end{align}

for non-negative integer $M$ , $\Psi _M$ being an $M$ th order polynomial in $\eta$ that satisfies

(7.3) \begin{align} -M \Psi _M + \eta \frac{d \Psi _M}{d \eta }=\frac{d^2 \Psi _M}{d \eta ^2}. \end{align}

We normalise $\Psi _M(\eta )$ by requiring that

(7.4) \begin{align} \Psi _M \sim \eta ^M \qquad \mbox{as} \quad |\eta | \rightarrow \infty, \end{align}

so that

(7.5) \begin{align} \Psi _0=1, \quad \Psi _1=\eta, \quad \Psi _2=\eta ^2-1, \quad \Psi _3=\eta ^3-3 \eta, \quad \Psi _4=\eta ^4-6 \eta ^2 +3. \end{align}

The coefficients in $\Psi _M$ alternate in sign, so that the constant term in $\Psi _{2N}$ has the same sign as $({-}1)^N$ , as does the linear term in $\Psi _{2N+1}$ , an observation that will be relevant in due course. In view of the power of $({-}t)$ in the prefactor in (7.2), self-consistency, in the sense that the second term in (5.1) is much smaller than the first in the limit $t \rightarrow 0^-$ , requires

(7.6) \begin{align} m M +(1-n)(M-2)\gt 0, \end{align}

which is automatically satisfied for $M \geq 1$ in the range of interest here but fails for $M=0$ ; the full complement of $M$ is required to furnish a complete set of eigenfunctions to (5.2), an issue to which we shall also implicitly need to return.

Similar issues arise for the linearisation about the steady state

(7.7) \begin{align} h=A_* x^{\frac{2}{m+1-n}}, \end{align}

for which $\bar{\mathcal H}$ , the perturbation about (7.7), satisfies the linear PDE (5.10), which is more conveniently written in the form

(7.8) \begin{align} \frac{\partial \bar{\mathcal H}}{\partial t} = A_* ^m \left (\frac{\partial ^2}{\partial x^2} \left ( x^{\frac{2m}{m+1-n}} \bar{\mathcal H} \right ) - \left (\frac{2 n (m+n+1)}{(m+1-n)^2}\right ) x^{-\frac{2(1-n)}{m+1-n}} \bar{\mathcal H} \right ). \end{align}

Since we limit ourselves to the range $m+3n+1\gt 0$ , the required backward self-similar solutions to (5.10) (i.e. those with maximal regularity as $\eta \rightarrow 0^+$ , this being a requirement for the subsequent matching) take the form

(7.9) \begin{align} \bar{\mathcal H}_M= ({-}t)^{\frac{n+1-m}{2(1-n)}} ({-}t)^M \eta^{\frac{n+1-m}{m+1-n}} \Theta _M(\eta ), \qquad \eta = A_*^{-\frac{m(m+1-n)}{2(1-n)}} x ({-}t)^{-\frac{m+1-n}{2(1-n)}}, \end{align}

where $\Theta _M(\eta )$ is a polynomial in $\eta^{\frac{2(1-n)}{m+1-n}}$ of $M$ th order, with

(7.10) \begin{align} -M \Theta _M + \frac{m+1-n}{2(1-n)} \eta \frac{d \Theta _M}{d \eta }= \eta^{-\frac{n+1-m}{m+1-n}}\frac{d^2 }{d \eta ^2} \left ( \eta^{\frac{m+1+n}{m+1-n}} \Theta _M \right )- \frac{2 n (m+1+n)}{(m+1-n)^2} \eta^{-\frac{2(1-n)}{m+1-n}} \Theta _M. \end{align}

We note that $\eta$ in (7.9) has the same time dependence as in (7.2), being that of the full similarity solution of Section 6.1, though the constant factor for convenience differs. We normalise $\Theta _M(\eta )$ according to

(7.11) \begin{align} \Theta _M \sim \eta^{\frac{2(1-n) M}{m+1-n}} \qquad \mbox{as} \quad \eta \rightarrow + \infty ; \end{align}

the coefficients again alternate sign, so that $\Theta _M(0)$ has the same sign as $({-}1)^M$ ; in contrast to $\Psi _M$ , for which the polynomials are alternately odd and even, $\Theta _M$ contains all the relevant powers. It follows that

(7.12) \begin{align} \Theta _0=1, \qquad \qquad \qquad \Theta _1=\eta^{\frac{2(1-n)}{m+1-n}}-2 \frac{(1-n)(n+m+3)}{(m+1-n)^2}, \end{align}

(7.13) \begin{align} \Theta _2=\eta^{\frac{4(1-n)}{m+1-n}}-4 \frac{(1-n)(m+5-n)}{(m+1-n)^2} \eta^{\frac{2(1-n)}{m+1-n}} +4 \frac{(1-n)^2(m+5-n)(n+m+3)}{(m+1-n)^4 }. \end{align}

Given that (7.7) is proportional to $({-}t)^{1/(1-n)}$ for $\eta =O(1)$ , self-consistency of (7.9) requires that

(7.14) \begin{align} m\lt (2M-1)(1-n) \end{align}

7.2. Extinction behaviour

We start with this case because it involves only one class of self-similar solutions (these being those of the limit problem in Section 4.2, rather than those of the full PDE) and does not involve any post-singular behaviour, $h$ being identically zero following extinction; moreover, it illustrates some principles that will be of more general relevance. We refer to [Reference Grundy13] for a very closely related approach to this case.

The extinction behaviour is ‘flat’ in the sense that it is selected by the linearisation (5.1), with $M$ even in (7.2); in other words, we have ‘inner’ expansion

(7.15) \begin{align} h \sim \left ((1-n) ({-}t)\right )^{\frac{1}{1-n}} - a_N ({-}t)^{\frac{n+(m+1-n)N}{1-n}} \Psi _{2N}(\eta ),, \end{align}

with $N\geq 1$ and where $a_N$ is an arbitrary constant that is required to be positive in view of the matching to $b_N$ noted below. In view of (7.4), the expression (7.15) disorders for $|x| \geq O(({-}t)^{\frac{1}{2N}})$ (i.e. the final ‘linearised’ term in (7.15) becomes comparable with the leading $\left ((1-n) ({-}t)\right )^{\frac{1}{1-n}}$ term). On this outer scale, the dominant balance in (1.1) is given by

\begin{align*} \frac{\partial h}{\partial t} \sim - h^n \end{align*}

so that (cf. Section 4.2)

(7.16) \begin{align} h \sim ({-}t)^{\frac{1}{1-n}} \left ((1-n) \left (1-b_N \frac{x^{2N}}{({-}t)} \right )\right )^{\frac{1}{1-n}} \end{align}

with $b_N \propto a_N$ ; the matching condition (‘inner of outer’) that follows from (7.16) can readily be inferred. From (7.16) it follows that the interfaces satisfy

\begin{align*} x \sim \pm \left ( \frac{({-}t)}{b_N} \right )^{\frac{1}{2N}} \qquad \mbox{as} \quad t \rightarrow 0^-, \end{align*}

wherein $b_N$ is an arbitrary positive constant.

We now discuss the status of $N$ in the analysis above in terms of the stability (reflecting the completeness of the eigenfunctions) and of the necessity for an entire family of possible asymptotic behaviours ( $N$ being an arbitrary positive integer). The case $N=1$ is expected to provide the generic (i.e. stable) extinction behaviour; while the solution to (5.2) will in general also contain ${\mathcal H}_0$ and ${\mathcal H}_1$ contributions, each of these being eliminated by identifying $x=0$ , $t=0$ with the place and time of extinction (it is no coincidence that there are two translation symmetries and that two eigenfunctions need to be suppressed); indeed, ${\mathcal H}_0$ corresponds to a small shift in $t$ in (1.6), while ${\mathcal H}_M \propto \partial{\mathcal H}_{M+1}/\partial x$ , so an $x$ translation in ${\mathcal H}_2$ can suppress ${\mathcal H}_1$ . Similar arguments apply in the singular scenarios outlined below, where we typically omit such details (see the end of Section 5.2, however). The necessity for the non-generic cases $N \geq 2$ can be argued as follows: for initial data containing a minimum between the two maxima, two generic behaviours can occur prior to extinction:

1. 1. diffusion may eliminate the minimum to leave a single maximum or

2. 2. the sink term may lead to the minimum touching down (as described in the next subsection), with the support of $h$ breaking into two disconnected sets.

The case $N=2$ (which has a minimum at $\eta =0$ ) describes intermediate-asymptotic behaviour on the borderline between these two, with yet higher $N$ corresponding to borderlines between scenarios that are themselves non-generic (the inevitability of the latter being implicit in the above identification of the necessity of a scenario with $N=2$ , enabling a bootstrapping argument to higher and higher $N$ given the presence of the comparable families associated with other types of singular behaviour).

In summary, as well as characterising a particularly significant type of intermediate-asymptotic behaviour, the above illustrates a central element of the subsequent subsections, namely the existence of a (discrete) set of asymptotic possibilities, the leading one of which is generic, this being particularly transparent in the current case given that the relevant solution to the heat equation can be written as a sum over the Hermite polynomials (5.5) with $k \geq 2$ . The number of modes to which a given asymptotic scenario is unstable is readily identified from such a representation in terms of the omitted contributions and the non-generic behaviours play an identifiable role in the dynamics: for example, and as noted above, the case $C_2 = C_3 = 0$ can be interpreted as providing the borderline between solutions whose support breaks up prior to extinction and those for which it does not.

7.3. Touchdown

It is noteworthy that there are three distinct intermediate-asymptotic scenarios in this case. Firstly, the results of the previous subsection can be immediately revisited with the only changes being $a_N \rightarrow -a_{N}$ and $b_N \rightarrow -b_{N}$ , with $a_N$ and $b_N$ then still positive. Hence, in particular, (7.16) becomes

(7.17) \begin{align} h \sim ({-}t)^{\frac{1}{1-n}} \left ( (1-n) \left ( 1 + b_N \frac{x^{2N}}{({-}t)} \right ) \right )^{\frac{1}{1-n}} \end{align}

so that the local profile at $t=0$ is given by

\begin{align*} h(x,0) \sim \left ( (1-n) b_N x^{2N} \right )^{\frac{1}{1-n}} \qquad \mbox{as} \quad |x| \rightarrow 0. \end{align*}

For $t \rightarrow 0^+$ (i.e. after the touchdown and film rupture at $t=0$ ), one then has the asymptotic behaviour

(7.18) \begin{align} h \sim \left ( (1-n) \left (b_N x^{2N} -t \right ) \right )^{\frac{1}{1-n}} \end{align}

with retracting interfaces given by

\begin{align*} x\sim \pm \left ( \frac{t}{b_N} \right )^{\frac{1}{2N}}, \end{align*}

and with $h \equiv 0$ for smaller $|x|$ ; since (7.18) is consistent with (2.2), the former needs no supplementation by an inner scaling.

Secondly, the green curves in Figure 7 are of the required type for the pre-touchdown local behaviour, while the associated post-touchdown local behaviour is then expected to be a forward similarity solution of the form (1.7).

Thirdly, and clarifying the bifurcation structure in Figure 7, the black branch therein, which gives the exceptional solution, can be followed in the sense of the linearisation (5.10), so that

(7.19) \begin{align} h \sim \left \{ \begin{array}{l@{\quad}l@{\quad}l} \displaystyle A_* x^{\frac{2}{m+1-n}} + c_M^+ \bar{\mathcal H}_M(x,t) & \mbox{for} & \eta \gt 0, \\[9pt] \displaystyle A_* ({-}x)^{\frac{2}{m+1-n}} + c_M^- \bar{\mathcal H}_M({-}x,t) & \mbox{for} & \eta \lt 0, \end{array} \right. \end{align}

provides the outer (i.e. $\eta =O(1)$ ) expansion, with arbitrary constants $c_M^{\pm }$ and $\bar{\mathcal H}_M$ and $\eta$ as in (7.9). Here, the pre-touchdown solution of Section 3.1, see (3.8), comes into play via its quasi-steady generalisation: the resulting inner scalings are

(7.20) \begin{align} x=({-}t)^M \xi, \qquad h=({-}t)^{\frac{2M}{m+1-n}} g(\xi ) \end{align}

implying the quasi-steady balance (as in Section 3.1)

(7.21) \begin{align} \frac{d }{d \xi } \left ( g^m \frac{d g }{d \xi } \right ) \sim g^n \end{align}

at leading order. The matching condition

(7.22) \begin{align} g \sim A_* |\xi |^{\frac{2}{m+1-n}} \qquad \mbox{as} \quad |\xi | \rightarrow \infty \end{align}

specifies the solution to (7.21) up to arbitrary constants $\xi _0$ and $\alpha$ in the form

(7.23) \begin{align} g=\alpha ^2 G\left ( \frac{\xi +\xi _0}{\alpha^{m+1-n}} \right ) \end{align}

where $G(\sigma )$ is an even function of $\sigma$ that can be specified uniquely via

(7.24) \begin{align} & G(\sigma ) \sim A_* \sigma^{\frac{2}{m+1-n}} +\sigma^{\frac{n+1-m}{m+1-n}} \qquad \mbox{as} \quad \sigma \rightarrow +\infty \qquad \mbox{for} \quad n \gt 0, \nonumber \\[5pt] & G(\sigma ) \sim A_* \sigma^{\frac{2}{m+1-n}} -\sigma^{\frac{n+1-m}{m+1-n}} \qquad \mbox{as} \quad \sigma \rightarrow +\infty \qquad \mbox{for} \quad n \lt 0,\end{align}

The switch in sign follows (as in (3.10) above) from an analysis of the integral that arises on solving the separable equation

\begin{align*} \frac{1}{2} G^{2m} \left ( \frac{d G }{d \sigma } \right )^2 -\frac{1}{m+n+1} \left ( G^{m+n+1} -G_0^{m+n+1} \right )=0 \end{align*}

that results as the first integral of (7.21), where

\begin{align*} G_0=G(0). \end{align*}

In the special case $n=0$ already noted above, we have

\begin{align*} G=\left ( G_0 + \frac{m+1}{2} \sigma ^2 \right )^{\frac{1}{m+1}} \end{align*}

and the correction terms in (7.24) are instead absent; we omit further analysis of this exceptional case here.

Given that

(7.25) \begin{align} \bar{\mathcal H}_M \sim K_M ({-}1)^M ({-}t)^M x^{\frac{n+1-m}{m+1-n}} \qquad \mbox{as} \quad \eta \rightarrow 0. \end{align}

for some positive constant $K_M$ that is to be regarded as known, given the analysis of Section 5.2, matching to (7.20) and (7.24) implies that $\alpha$ and $\xi _0$ are given by

(7.26) \begin{align} & \frac{2}{m+1-n} \xi _0 A_* \pm \alpha^{m+1-n} = ({-}1)^M c_M^+ K_M \nonumber\\[5pt]& -\frac{2}{m+1-n} \xi _0 A_* \pm \alpha^{m+1-n} = ({-}1)^M c_M^- K_M \end{align}

where the signs on the left-hand sides correspond to those in (7.24). Self-consistency of such a scenario requires that (7.14) hold. Before turning to questions of stability, it is helpful to make connections between (7.26) and Figure 7; the green branches in the latter bifurcate off when equality holds in (7.14) and are symmetric. For $m$ slightly larger than $(2M-1)(1-n)$ , the expressions (7.26) remain relevant with $\xi _0=0$ , $c_M^+=c_M^-=c_M$ ; since $\alpha$ is necessarily positive (a constraint with broader implications in (7.26)), this implies $c_M\gt 0$ and hence $A\gt A_*$ for $n\lt 0$ , $M$ odd and $n\gt 0$ , $M$ even, with $c_M\lt 0$ and $A\lt A_*$ in the converse cases. This is consistent both with the alternating directions of the green branches in Figure 7 and with their swapping sides as $n$ passes through zero. Turning now to stability, $N=1$ in (7.17) leads to a generic cases, as in the previous subsection; $N \geq 2$ scenarios are again necessary but unstable. We signpost that the cases akin to (7.17), but with $M$ odd, will also be relevant in due course. In terms of (7.19), the cases $\eta \gt 0$ and $\eta \lt 0$ in effect operate independently, so suppressing the $M=0$ mode in (5.10) on both sides uses up both translational invariants. Since (7.14) is violated for $M=1$ , $m+n\gt 1$ , the first legitimate instance of (7.19) is doubly unstable (though in the symmetric case, or with zero-Neumann boundary condition on $x=0$ (for which $x$ translational invariance cannot be exploited), it is only singly so). However, for related reasons the left-most green branch in panel of Figure 7 with $n \leq 0$ is expected to be stableFootnote ⁴ above the fold, but unstable below it: in these regimes there are thus two stable scenarios, that in Figure 7 for a narrow range of $m\gt 1-n$ , while (7.17) with $N=1$ is applicable throughout. In Figure 7, successive branches (of any colour) are expected to gain an additional unstable mode on each occasion a bifurcation point is passed through from left to right.

At $t=0$ , the spatial profiles in the third scenario take the form

(7.27) \begin{align} h \sim A_* |x|^{\frac{2}{m+1-n} } + c_M^{\pm } |x|^{\frac{n+1-m+2M(1-n)}{m+1-n} } \qquad \mbox{as} \quad |x| \rightarrow 0. \end{align}

There seems subsequently to be the possibility of non-uniqueness, but we shall not discuss the associated post-touchdown behaviour given that these are in any case unstable.

7.4. Attachment, detachment and coalescence

Attachment, either to a fixed boundary or through two interfaces colliding (either through both advancing or by one advancing faster than the other is retracting) – termed coalescence above – raises no pre-singularity questions since the solution has no prior knowledge that it is about to attach (unlike the other cases analysed here).Footnote ⁵ Post-attachment raises no significant challenges either: one has a porous-medium equation similarity solution of the form

\begin{align*} h \sim t^{\frac{1}{m}} f \left ( \frac{x}{t} \right ) \qquad \mbox{as} \quad \mbox{with} \quad x=O(t), \end{align*}

with additional structure when an advancing interface overuns a receding one of a form akin to that analysed in the generic anti-reverser case addressed below. One related observation is in order; however, the pre-attaching behaviour of an advancing interface takes the form

\begin{align*} h \sim ({-}t)^{\frac{1}{m}} f \left ( \frac{x}{({-}t)} \right ), \qquad f^m(\eta )=m\eta _0(\eta _0-\eta ) \end{align*}

for some constant $\eta _0$ . This simply corresponds to a Taylor expansion in $t$ of (2.3), reflecting the lack of prior knowledge referred to above, but is also an exceptional connection for the porous-medium equation in a sense elaborated upon in [Reference Hulshof, King and Bowen14], the significance of which here is that this relates to the probably unexpected branches present in the bottom left-hand corners of Figure 7, though we shall not go into details here.

There are of course non-generic attaching behaviours, such as an interface reversing at the same instant as attaching, but the no-prior-knowledge argument applies to these also. For the hole filling case in higher dimensions, the situation is more complicated, however, since the interface is aware of its increasing curvature; see, for example, [Reference Angenent, Aronson, Betelu and Lowengrub1] for the porous-medium equation – including a sink term, as in (1.1), leads to open questions in this regard.

The detachment case is also simple in the sense that we conjecture that the left-most yellow branch in each panel of Figure 7 provides the generic behaviour, with post-detachment being described by an associated forward similarity solution. For $n\gt 0$ , we expect this branch to start from $A\lt A_*$ as $m+n \rightarrow 1^+$ (with a fold in $m+n\lt 1$ connecting the branch to $A=A_*$ at $m+n=1$ ), whereas for $n\leq 0$ we have $A \rightarrow A_*$ in that limit. Non-generic cases arise both as the other yellow branches and from tracking along the black branch prior to the secondary and subsequent bifurcations, as in the first of (7.19), being subject to (7.14). Now a quasi-steady version of the leaking solution of Section 3.1, satisfying (3.6), describes the inner behaviour – the inner scalings are again given by (7.20), but now (7.21) is subject to

\begin{align*} g=0 \qquad \mbox{at} \quad \xi =0 \end{align*}

so that

\begin{align*} g=\alpha ^2 G \left ( \frac{\xi }{\alpha^{m+n-1}} \right ) \end{align*}

with

\begin{align*} \frac{1}{2} G^{2m} \left ( \frac{d G}{d \sigma } \right )^2 -\frac{1}{m+1} G^{m+n+1} =\frac{1}{2} J_0^2 \end{align*}

where

\begin{align*} J_0=G^m \left. \frac{d G}{d \sigma } \right |_{\sigma =0} \gt 0 \end{align*}

determines the outward flux. For $n=0$ , we have

\begin{align*} G=\left ( (m+1) \left (J_0 \sigma + \frac{1}{2} \sigma ^2 \right )\right )^{\frac{1}{m+1}} \end{align*}

and for any $n$ the sign in the first of (7.24) applies here – correspondingly, the yellow branches of Figure 7 are on the opposite sides from the green ones for $n\lt 0$ and on the same side for $N\gt 0$ . The local behaviour at detachment is then again as in (7.27). Instability follows even for $M=2$ because the $x$ -translation mode is not available given the boundary condition (i.e. $h=0$ at $x=0$ ), and the $M=0$ and $M=1$ modes both need to be absent; accordingly, we omit any discussion of the post-detachment behaviour in these non-generic cases.

7.5. Reversers and anti-reversers

We start with scenarios in which the two have closely related behaviour, before discussing a (generic) anti-reversing structure distinct from the other categories arising here. We take the support to lie in $x\gt 0$ at $t=0$ and start from the first of (7.19) (writing $c_M=c_M^+$ ) with, as always, (7.14) being required. Since here we can apply both $t$ and $x$ translational invariance, the $M=0$ and $M=1$ modes can both be suppressed (in contrast to the situation in the previous section), leaving $M=2$ as a potential stable scenario for $m\lt 3(1-n)$ . The matching condition

(7.28) \begin{align} h \sim A_* x^{\frac{2}{m+1-n}} + ({-}1)^M K_M c_M ({-}t)^M x^{\frac{n+1-m}{m+1-n}} \end{align}

then applies on a travelling-wave inner region, for which a translation

(7.29) \begin{align} x=s(t)+z, \qquad s(t) \sim \frac{({-}1)^{M+1} K_M c_M}{A_*} ({-}t)^M \end{align}

is required, with the dominant balance being (as in Section 3.2)

(7.30) \begin{align} -\dot{s} \frac{\partial h}{\partial z} \sim \frac{\partial }{\partial z} \left ( h^m \frac{\partial h}{\partial z} \right ) -h^n, \end{align}

which implies that the inner scalings are

(7.31) \begin{align} z=O \left ( ({-}t)^{\frac{(m+1-n)(M+1)}{m+n-1}} \right ), \qquad h=O \left ( ({-}t)^{\frac{2(M-1)}{m+n-1}} \right ). \end{align}

Given that the outer scaling associated with (7.19) is the full self-similar one, $x=O \left ( ({-}t)^{\frac{m+1-n}{2(1-n)}} \right )$ , consistency conditions on the above read

(7.32) \begin{align} \frac{(m+1-n)(M-1)}{m+n-1} \gt M \gt \frac{m+1-n}{2(1-n)}, \end{align}

both inequalities simply reproducing (7.14). Whether (7.28) is associated with a reverser or an anti-reverser depends upon the sign of $c_M$ ; if $c_M ({-}1)^{M+1}$ is positive the above asymptotic behaviour corresponds to a reverser (advancing interface for $t\lt 0$ ), while if it is negative it is an anti-reverser (receding interface for $t\lt 0$ ),Footnote ⁶ subject to the following qualification – for the stable case $M=2$ and other even $M$ the interface indeed changes direction (see (7.29)), while for $M$ odd it recommences its original direction after an instantaneous pause (the singly unstable case $M=3$ corresponding to the borderline between (i) reversers and anti-reversers occurring in rapid succession (in either order) and (ii) no change of direction).

A first integral of (7.30) is available in the case $n=0$ , namely

\begin{align*} -\dot{s} h = h^m \frac{\partial h}{\partial z} -z. \end{align*}

With regard to the transcritical red/blue bifurcations in Figure 7, for $n\lt 0$ the numerics provide evidence that the primary such bifurcation has red unstable and blue stable close to the bifurcation point, though with folds rapidly swapping the stability of each, whereas for $n\gt 0$ the red branch is expected to be stable throughout and the blue unstable (except perhaps in the lower left-hand corner: we leave the stability of the branches there as an open question). Similar comments apply with respect to the number of unstable model arising from the subsequent bifurcations.

Collecting together the implications of the above for the generic cases, we summarise our conjectures as follows.

(I) Reversers

For $m\lt 3(1-n)$ , (7.28)–(7.30) with $M=2$ provide a stable solution, changing direction smoothly (i.e. $s \propto ({-}t)^2$ in (7.29)). At $t=0$

(7.33) \begin{align} h \sim A_* x^{\frac{2}{m+1-n}} + \kappa x^{\frac{5-3n-m}{m+1-n}} \end{align}

holds for some constant $\kappa$ and the $t \rightarrow 0^+$ behaviour follows directly from the $t \rightarrow 0^-$ behaviour (the equivalent polynomial solutions provide the correction term, though their coefficients then all take the same sign). For $n\gt 0$ , the generic behaviour in $m\gt 3(1-n)$ is a full similarity solution, corresponding to the first red branch in Figure 7. For $n\lt 0$ , the situation is slightly more complicated, with a bistable regime being present over a small range of $m$ in $m\lt 3(1-n)$ , the red branch being unstable below the fold (being the borderline between the two stable scenarios) but stable above it and continuing into the range $m\gt 3(1-n)$ . When a full similarity solution applies, $s$ is in general not smooth:

(7.34) \begin{align} s \sim S_- ({-}t)^{\frac{m+1-n}{2(1-n)}} \quad \mbox{as} \quad t \rightarrow 0^-, \qquad s \sim S_+ t^{\frac{m+1-n}{2(1-n)}} \quad \mbox{as} \quad t \rightarrow 0^+, \end{align}

with a forward similarity solution holding as $t \rightarrow 0^+$ and with the constants $S_-$ and $S_+$ presumably being different and positive on the primary branch.

(II) Anti-reversers

The description in (I) above again applies for $m\lt 3(1-n)$ , $M=2$ (whether a reverser or anti-reverser pertains depends on the sign of $\kappa$ ). For $n\lt 0$ there is expected to be a narrow range of bistability in $m\gt 3(1-n)$ also (one of the stable scenarios being the full similarity solution, as in Figure 7, the other being described next), but the first blue branch is otherwise unstable, necessitating the identification of a different class of intermediate-asymptotic behaviour. The range $m\lt 3(1-n)$ is thus bistable, with the blue branch expected to provide the borderline between the two generic cases.

An important clue into the other generic scenario is provided by a comparison between (2.2) and (2.3): since

\begin{align*} 1-n\lt m \quad \mbox{for} \quad m+n\gt 1, \end{align*}

retreating interfaces are associated with local behaviour smaller than that of the advancing ones, so we may anticipate that the latter can overrun the former with little modification to itself. We now formalise that intuition; related analyses of the porous-medium equation, [Reference King and Please19, Reference Lacey20], represent precedents. Two outer regions are present on the scaling $x=O({-}t)$ , namely

(7.35) \begin{align} \begin{array}{lll} \displaystyle h^m \sim m q (x-q({-}t)) &\mbox{for} & \displaystyle \frac{x}{({-}t)} \gt q,\\ \displaystyle h^{1-n} \sim (1-n) \left ({-}t+\frac{x}{Q} \right ) &\mbox{for} & \displaystyle q\gt \frac{x}{({-}t)} \gt -Q, \end{array} \end{align}

wherein $q$ and $Q$ are each positive constants, the former expression being a similarity solution of the porous-medium equation (as in Section 4.1) and the latter being one of the pure absorption limit (as in Section 4.2). The two expressions in (7.35) are comparable for

\begin{align*} x=q({-}t)+ O\left ( ({-}t)^{\frac{m}{1-n}} \right ), \end{align*}

self-consistency here again requiring that $m+n\gt 1$ , which implies the inner scalings

\begin{align*} x=q({-}t) + ({-}t)^{\frac{m}{1-n}} \xi, \qquad h=({-}t)^{\frac{1}{1-n}} g, \end{align*}

leading to

\begin{align*} q \frac{\partial g}{\partial \xi } \sim \frac{\partial }{\partial \xi } \left ( g^m \frac{\partial g}{\partial \xi } \right ) \end{align*}

and hence to

\begin{align*} q \left ( g- \left ( (1-n) \left ( 1+ \frac{q}{Q} \right )\right )^{\frac{1}{1-n}}\right ) \sim g^m \frac{\partial g}{\partial \xi }, \end{align*}

thereby matching with (7.35) in both limits $|\xi | \rightarrow \infty$ . In consequence, we simply have

\begin{align*} h(x,0) \sim ( m q x)^{\frac{1}{m}} \qquad \mbox{as} \quad x \rightarrow 0^+ \end{align*}

so the interface obeys

(7.36) \begin{align} x \sim Q t \qquad \mbox{as} \quad t \rightarrow 0^-, \qquad x\sim -q t \qquad \mbox{as} \quad t \rightarrow 0^+, \end{align}

its velocity being discontinuous (indeed swapping sign) at $t=0$ .

There is one final, but unstable, case to which we have already alluded and for which relevant results can be read off directly from (7.17) on replacing the even exponent by an odd one, so that

\begin{align*} h \sim ({-}t)^{\frac{1}{1-n}} \left ( (1-n) \left ( 1+b_N \frac{x^{2N+1}}{({-}t)} \right ) \right )^{\frac{1}{1-n}} \end{align*}

with $N \geq 1$ and with $b_N\gt 0$ without loss of generality. This has a retreating interface for $t\lt 0$ at

\begin{align*} x\sim -\left ( \frac{({-}t)}{b_N} \right )^{\frac{1}{2N+1}}, \end{align*}

with $h \equiv 0$ for smaller $x$ , and

\begin{align*} h(x,0) \sim \left ( (1-n) b_N x^{2N+1} \right )^{\frac{1}{1-n}} \qquad \mbox{as} \quad x \rightarrow 0^+. \end{align*}

Since the interface has

\begin{align*} x \sim \left ( \frac{t}{b_N} \right )^{\frac{1}{2N+1}} \qquad \mbox{as} \quad t \rightarrow 0^+ \end{align*}

in this case, no change of direction occurs, with $N =1$ being singly unstable and providing the borderline between solutions that lose a local maximum/minimum pair while retreating and those which touchdown. These cases are striking in having unbounded, rather than zero, interface velocity at $t=0$ and do not fall naturally under any of the headings above.

8.1. Validation of the asymptotic predictions

Here we provide evidence to substantiate the observability of some of the local descriptions identified as generic above. First, we summarise the findings of the numerical experiments, and then, we give the details on how solutions to (1.1) were found by direct numerical simulation. All results in this section have $n=0$ but similarly good agreement is expected for other values of $n$ in the range (1.2).

(i) Reversing events

Panel (ia) in Figure 9 indicates that for $m\lt 3(1-n)$ we have $s \sim \text{const.} ({-}t)^2$ as $t \to 0^-$ with the constant dependent upon the initial conditions. This is in agreement with (7.33). Panel (ib) indicates that for $m\gt 3(1-n)$ we have $s \sim \hat{\xi } ({-}t)^{(m+1-n)/2(1-n)}$ where the constant $\hat{\xi }$ is that obtained from the ODE solutions described in Section 6 on the stable red branch: we find

(8.1) \begin{align} \hat{\xi } \approx 0.3864 \quad \mbox{for} \quad m=4, \quad n=0, \end{align}

(8.2) \begin{align} \hat{\xi } \approx 0.5014 \quad \mbox{for} \quad m=5, \quad n=0, \end{align}

the PDE simulations being consistent with these values up to three significant figures.

Figure 9. Comparison between the asymptotic results, shown in grey, and direct numerical simulation of (1.1), shown in black, for $n=0$ . Panels (ia) and (ib) show the interface position leading up to a reversing event for different choices of initial condition. In (ia) dashed curves are for $m=2$ , whereas solid ones are $m=3$ . In (ib), dotted curves are for $m=4$ and solid curves for $m=5$ . Panels (ii) and (iii) show the interface position leading up to anti-reversing and attachment events, respectively, for different choices of initial condition. In both panels (ii) and (iii), dashed curves indicate $m=2$ , solid curves $m=3$ and dotted curves $m=4$ .

(ii) Anti-reversing events

Panel (ii) in Figure 9 indicates that the interface obeys $s \sim \text{const.} ({-}t)$ with the constant selected by the initial data, in agreement with (7.36). We omit here the other regime of stable anti-reversing identified in Section 7.5, namely pure self-similarity.

(iii) Attaching events

Panel (iii) in Figure 9 indicates that the interface obeys $s \sim \text{const.} ({-}t)$ with the constant selected by the initial data in agreement with the discussion in Section 7.4; in short, this is expected because the interface has no way to anticipate its collision with the absorbing boundary.

(iv) Detaching events

Panel (iv) in Figure 10 indicates that the outward flux at the interface leading up to a detachment event obeys

(8.3) \begin{align} \left. h^m \frac{\partial h}{\partial x} \right |_{x=0} = ({-}t)^{(m+1)/2} \left. f^m \frac{\partial f}{\partial \xi } \right |_{\xi =0} \end{align}

where the quantity $\left. f^m{\partial f}/{\partial \xi } \right |_{\xi =0}$ can be determined via the shooting method described in Section 6. Shooting on the self-similar ODE yields

(8.4) \begin{align} \left. f^m \frac{\partial f}{\partial \xi } \right |_{\xi =0} = 0.4912 \quad \mbox{for} \quad m=2, \quad n=0, \end{align}

(8.5) \begin{align} \left. f^m \frac{\partial f}{\partial \xi } \right |_{\xi =0} = 0.8786 \quad \mbox{for} \quad m=3, \quad n=0, \end{align}

(8.6) \begin{align} \left. f^m \frac{\partial f}{\partial \xi } \right |_{\xi =0} = 1.3076 \quad \mbox{for} \quad m=4, \quad n=0, \end{align}

with which the PDE solutions are again consistent up to three significant figures. This verifies the claim in Section 7.4 that the left-most yellow branch in Figure 7 provides the generic description of detachment.

Figure 10. Comparison between the asymptotic results, shown in grey, and direct numerical simulation of (1.1), shown in black, for $n=0$ ; results for $m=2$ are indicated by dashed curves, for $m=3$ with solid curves and for $m=4$ with dotted curves. Panel (iv) shows the flux at the interface leading up to a detachment event for different choices of initial condition. Panels (va) and (vb) show $h$ and $h_{xx}$ at the minimum leading up to a touchdown event.

(v) Touchdown events

Panels (va) and (vb) in Figure 10 indicate that the description leading up to a touchdown event are in agreement with (7.18) with $N=1$ , verifying that this is the generic behaviour; we note that the curvature does not tend to a constant at $t \to 0^-$ in cases where $N\gt 1$ .

(vi) Coalescence events

Panel (vi) in Figure 10 indicates that the interfaces obey $s \sim \text{const.} ({-}t)$ with the constant selected by the initial data, in agreement with the discussion in Section 7.4 and again following from the no-prior-knowledge argument.

8.2. Numerical methods

We furnish numerical solutions (1.1) on a finite domain by supplementing it with an interface condition at the left-hand end of its compact support, that is, one of the conditions (1.3)–(1.5), as well as a symmetry condition at some location $x=x_0\gt 0$ , that is,

(8.7) \begin{align} \left. \frac{\partial h}{\partial x} \right |_{x=x_0} =0. \end{align}

The value of $x_0$ is chosen arbitrarily; the spatial coordinates are translated post hoc such that the singular event occurs at $x=0$ .

Numerical solution is made awkward by two features, namely (a) that the domain is dynamic and (b) that the solution often exhibits an infinite spatial gradient at the interface. We alleviate these difficulties by transforming the spatial coordinates using a similar strategy to that in [Reference Foster and Pelinovsky6]. As we have seen, the solution local to an interface takes the form

(8.8) \begin{align} h \sim \lambda (t) (x-s(t))^{\omega } \quad \mbox{as} \quad x \to s(t)^+ \end{align}

for some $\lambda$ , where $\omega =1/m,1/(1-n),1/(m+1)$ or $2/(m+1-n)$ depending on whether the interface is advancing, receding, in contact with an absorbing boundary or stationary and mass preserving. These local behaviours motivate the coordinate transformations

(8.9) \begin{align} y = \left ( 1-\frac{x_0 - x}{x_0 - s(t)} \right )^\omega \end{align}

such that the behaviour (8.8) is replaced with $h \sim \hat{\lambda }(t) y$ as $y \to 0^+$ . The PDE (1.1) becomes

(8.10) \begin{align} \frac{\partial h}{\partial t} = \frac{wy^{(w-1)/w}}{x_0-s} \frac{\partial }{\partial y} \left ( \frac{wy^{(w-1)/w}}{x_0-s} h^m \frac{\partial h}{\partial y} \right ) + \frac{w \left (1-y^{1/w}\right ) y^{(w-1)/w}}{x_0 -s} \frac{ds}{dt} \frac{\partial h}{\partial y} - h^n. \end{align}

This is closed by the boundary conditions

(8.11) \begin{align} h |_{y=0} =0, \quad \left. \frac{\partial h}{\partial y}\right |_{y=1} =0, \end{align}

where the latter results from transforming (8.7) and the former results from (1.3a) for a mobile interface or (1.5) for an absorbing boundary. When an absorbing boundary is being considered we have $\dot{s}=0$ , but for a dynamic interface the motion is solved for as part of the solution to the problem using the transformed counterpart of (1.4), which reads

(8.12) \begin{align}{\dot{s} = \lim _{y \searrow 0} \left \{ \begin{array}{ll} \displaystyle -\frac{wy^{(w-1)/w}}{x_0-s} h^{m-1} \frac{\partial h}{\partial y} & \textrm{if} \;\; \dot{s}(t) \leq 0, \\[9pt] \displaystyle \frac{h^n (x_0-s)}{w y^{(w-1)/2}} \left ( \frac{\partial h}{\partial y} \right )^{-1} & \textrm{if} \;\; \dot{s}(t) \geq 0. \end{array} \right .} \end{align}

The spatial derivatives in the PDE (8.10) and the boundary conditions (8.11) (along with (8.12) in scenarios where an interface is moving) are discretised using finite differences. Second-order central differences are used on internal points, while the ghost-point method is used to incorporate (8.11). To retain second-order accuracy, the spatial derivative in (8.12) is computed using a forward-finite-difference approximation involving values of $h$ at the node in question and two neighbouring nodes. Once this has been done, the resulting large system of coupled ODEs is integrated forward in time using an adaptive Runge–Kutta–Fehlberg (RKF45) scheme implemented in MATLAB.