2.1. Characteristic scales of the system

Before progressing with our analysis, we identify the characteristic scales of the system, and construct the equivalent non-dimensional model. We first write the permeability as $k=\overline {k}\kappa$ , where $\overline {k}$ is the characteristic permeability, which here we take to be the arithmetic mean across the whole domain of width $2W$ , so that $\overline {k}=({1}/{2W})\int _{-W}^W k(y) \, \mathrm{d}y$ . For ease of notation, we introduce the mean buoyancy velocity

(2.7) \begin{equation} u_b = \frac {\Delta \rho g \overline {k}}{\mu }. \end{equation}

We use the characteristic wavelength $\lambda$ of the permeability heterogeneity to scale lengths in the $x$ - and $y$ -directions, hence a scaling analysis of (2.4) and (2.5) leads to the non-dimensional scales

(2.8) \begin{equation} x,y\sim \lambda ,\quad \boldsymbol{q} \sim q, \quad h\sim \left (\frac {q}{u_b}\right )^{1/2} \lambda ^{1/2}, \quad t\sim \phi \left (\frac {1}{q u_b}\right )^{1/2}\lambda ^{3/2}. \end{equation}

We use these to rescale the variables as

(2.9) \begin{equation} \tilde {x} = \frac {1}{\lambda } x, \quad \tilde {y} = \frac {1}{\lambda } y, \quad \tilde {\boldsymbol{q}} = \frac {1}{q}\boldsymbol{q}, \quad \tilde {h} = \left (\frac {u_b}{q \lambda }\right )^{1/2} h, \quad \tilde {t} = \frac {\left (q u_b\right )^{1/2}}{\phi \lambda ^{3/2}} t, \quad \kappa = \frac {1}{\overline {k}}k, \end{equation}

where the dimensionless quantities are denoted with a tilde. We therefore arrive at the following non-dimensional system:

(2.10) \begin{equation} \frac {\partial \tilde {h}}{\partial \tilde {t}} = {\boldsymbol{\nabla }} \boldsymbol{\cdot } \big(\kappa \tilde {h}\, {\boldsymbol{\nabla }} \tilde {h} \big), \end{equation}

subject to

(2.11a) \begin{align} \tilde {\boldsymbol{q}}\boldsymbol{\cdot }\hat {\boldsymbol{x}} = 1 \quad &\mathrm{at} \quad \tilde {x}=0, \end{align}

(2.11b) \begin{align} \tilde {h} = 0 \quad &\mathrm{at} \quad \tilde {x}=\tilde {x}_N, \end{align}

(2.11c) \begin{align} \tilde {\boldsymbol{q}}\boldsymbol{\cdot }\hat {\boldsymbol{n}} = 0 \quad &\mathrm{at} \quad \tilde {x}=\tilde {x}_N, \\[9pt] \nonumber \end{align}

which is forced by the non-dimensional permeability modulation $\kappa (y)$ . From here on, we drop the tildes for notational convenience, and all quantities can be assumed to be non-dimensional unless otherwise stated.

4.1. Homogeneous permeability

First, we characterise the mean behaviour by considering the case of a homogeneous porous medium with $\kappa =1$ everywhere. The height is therefore uniform in $y$ , reducing the system described by (2.10) and (2.11) to a one-dimensional problem in $x$ :

(4.1) \begin{equation} \frac {\partial h}{\partial t} = \frac {\partial }{\partial x}\left ( h \frac {\partial h}{\partial x} \right )\!, \end{equation}

subject to the three boundary conditions

(4.2a) \begin{align} -h\frac {\partial h}{\partial x} &= 1 \quad \mathrm{at} \quad x=0, \end{align}

(4.2b) \begin{align} h &= 0 \quad \mathrm{at} \quad x=x_N, \end{align}

(4.2c) \begin{align} -h\frac {\partial h}{\partial x} &= 0 \quad \mathrm{at} \quad x=x_N. \end{align}

The solution to this problem is well known and self-similar, as described in Huppert & Woods (Reference Huppert and Woods1995) for constant injection rate (their $\alpha =1$ ). Their solution, reviewed here, is written in terms of the self-similar variable and height scaling

(4.3a) \begin{align} \zeta &= x t^{-2/3}, \end{align}

(4.3b) \begin{align} h &= t^{1/3}\,f(\zeta ). \end{align}

Here, $f(\zeta )$ describes the profile of the current in self-similar space, and is determined by the solution of

(4.4) \begin{equation} \frac {\mathrm{d}}{\mathrm{d}\zeta }\left (f\frac {\mathrm{d}f}{\mathrm{d}\zeta }\right ) + \frac {2}{3}\zeta \frac {\mathrm{d}f}{\mathrm{d}\zeta } - \frac {1}{3}f = 0, \end{equation}

subject to

(4.5a) \begin{align} -f\frac {\mathrm{d}f}{\mathrm{d}\zeta } &= 1 \quad \mathrm{at} \quad \zeta = 0, \end{align}

(4.5b) \begin{align} f &= 0 \quad \mathrm{at} \quad \zeta = \overline {\zeta }_N, \end{align}

(4.5c) \begin{align} -f\frac {\mathrm{d}f}{\mathrm{d}\zeta } &= 0 \quad \mathrm{at} \quad \zeta = \overline {\zeta }_N, \end{align}

where $\overline {\zeta }_N$ is the nose of the current in self-similar coordinates. For constant injection rate, no closed-form solution exists, and the form of $f(\zeta )$ must be obtained numerically. Such a numerical solution provides the profile and the values for the nose position and height at injection:

(4.6a) \begin{align} \overline {\zeta }_N &\approx 1.48, \end{align}

(4.6b) \begin{align} f(0) &\approx 1.30. \end{align}

Towards the nose, $f$ is approxiately linear:

(4.7) \begin{equation} f \approx \frac {2}{3} \overline {\zeta }_N(\overline {\zeta }_N - \zeta ). \end{equation}

This approximation is used subsequently in some cases as a simple approximation of the mean structure near the nose in order to understand the structure of perturbations to the base state. The structure of the whole homogeneous solution is reviewed here for reference in the perturbed solutions described later.

4.2. Response to $\kappa (y)$ in physical space

We begin by setting up the problem in physical space, then perform the perturbation analysis in similarity space, where we can exploit the constant mean behaviour. It is instructive to initially describe the influence of permeability variations on the flow physically, in the original coordinate system in which the mean flow continually evolves.

To make analytical progress, we introduce a small sinusoidal perturbation of amplitude $\varepsilon$ to the permeability field, mirroring the set-up of numerical simulations described in § 3. The permeability perturbations induce corresponding perturbations to the height $h$ and nose $x_N$ of the current, so that

(4.8a) \begin{align} \kappa &= 1 + \varepsilon \cos (2\pi y), \end{align}

(4.8b) \begin{align} h(x,y,t) &= \overline {h}(x,t) + h^{\prime}(x,t)\cos (2\pi y), \end{align}

(4.8c) \begin{align} x_N(y,t) &= \overline {x}_N(t) + x^{\prime}_N(t)\cos (2\pi y). \end{align}

(The permeability field in the simulations includes a phase shift (3.1), but since it does not impact the form of our solutions, we have chosen to exclude it from our analysis.) Here, $\overline {h}(x,t)$ and $\overline {x}_N(t)$ are mean properties, and $h^{\prime}$ and $x^{\prime}_N$ are perturbations to the height and the nose, respectively. In the following analysis, all perturbed quantities are taken to be small compared with the corresponding leading-order term, $\varepsilon \ll 1$ , $h^{\prime}\ll \overline {h}$ , $x^{\prime}_N \ll \overline {x}_N$ .

By construction, the leading-order terms solve the homogeneous permeability system

(4.9a) \begin{align} \overline {h} &= t^{1/3}\,f(x t^{-2/3}), \end{align}

(4.9b) \begin{align} \overline {x}_N &= \overline {\zeta }_N t^{2/3}, \end{align}

as described in § 4.1. There is good agreement between the self-similar mean height and calculated transverse mean height, as seen in figure 2. Substitution of (4.8) into (2.10), and keeping only the leading-order corrections to the unperturbed system, provides a description of the evolution of the height perturbations:

(4.10) \begin{equation} \frac {\partial h^{\prime}}{\partial t} = \frac {\partial ^2}{\partial x^2}\left ( \overline {h}h^{\prime} \right ) - 4\pi ^2\overline {h}h^{\prime} + \varepsilon \frac {\partial }{\partial x}\left ( \overline {h}\frac {\partial \overline {h}}{\partial x}\right ). \end{equation}

The flux condition at the origin (2.11a ) becomes

(4.11) \begin{equation} \frac {\partial }{\partial x} \left ( \overline {h} h^{\prime} \right ) = \varepsilon \quad \mathrm{at} \ x=0. \end{equation}

The location of the perturbed nose, $x_N$ , is unknown a priori, so we expand the linearised boundary condition $h(x_N)=0$ about the mean nose position $\overline {x}_N$ . The first-order terms give the condition

(4.12) \begin{equation} h^{\prime} = - x_{N}^{\prime} \frac {\partial \overline {h}}{\partial x}\quad \mathrm{at} \ \mathrm{} x=\overline {x}_N. \end{equation}

Finally, we note that global mass conservation along with a scaling of the amplitude of perturbations requires that the net flux across the mean nose position must be zero. Since we have a periodic perturbation in $y$ , we can consider this over just one wavelength:

(4.13) \begin{equation} \int _0^1 -\kappa (y) \left [ h \frac {\partial h}{\partial x}\right ]_{x=\overline {x}_N} \, \mathrm{d}y = 0. \end{equation}

We use the expressions for permeability (4.8a ) and height (4.8b ) to expand the integrand to get

(4.14) \begin{align} &\int _0^1 \left [ \overline {h}\frac {\partial \overline {h}}{\partial x} + \cos ({2\pi y }) \left ( \frac {\partial }{\partial x}\{ \overline {h} h^{\prime} \} +\varepsilon \overline {h}\frac {\partial \overline {h}}{\partial x} \right ) \right . \nonumber \\&\quad +\left .\cos ^2({2\pi y}) \left ( h^{\prime}\frac {\partial h^{\prime}}{\partial x} + \varepsilon \frac {\partial }{\partial x}\{ \overline {h} h^{\prime} \}\right ) + \text{h.o.t.} \right ]_{x=\overline {x}_N} \mathrm{d}y = 0. \end{align}

The leading-order contribution vanishes since the mean height goes to zero at the mean nose position. For the remaining terms, the only $y$ dependence is found in the $\cos (2\pi y)$ terms. The first-order contribution therefore also automatically satisfies the global mass conservation condition. The second-order term simplifies to

(4.15) \begin{equation} \frac {\partial h^{\prime}}{\partial x} = - \varepsilon \frac {\partial \overline {h}}{\partial x}\quad \mathrm{at} \ x=\overline {x}_N, \end{equation}

which closes the system.

4.3. Response to $\kappa (y)$ in similarity space

The mean behaviour is well characterised in terms of the similarity variable $\zeta =xt^{-2/3}$ , as described in (4.3a ). Hence we transform the perturbed equations to $(\zeta ,y,t)$ -space such that

(4.16a) \begin{align} \kappa &= 1 + \varepsilon \cos (2\pi y), \end{align}

(4.16b) \begin{align} h(\zeta ,y,t) &= t^{1/3}\,f(\zeta ) + h^{\prime}(\zeta ,t)\cos (2\pi y), \end{align}

(4.16c) \begin{align} \zeta _N(y,t) &= \overline {\zeta }_N + \zeta '_N(t) \cos (2\pi y). \\[9pt] \nonumber \end{align}

Thus the mean height and nose position are both independent of time, simplifying the analysis of first-order effects.

The perturbed system therefore becomes

(4.17) \begin{equation} t\frac {\partial h^{\prime}}{\partial t} = \frac {2}{3}\zeta \frac {\partial h^{\prime}}{\partial \zeta } + \frac {\partial ^2}{\partial \zeta ^2}(f h^{\prime}) - 4 \pi ^2 t^{4/3} f h^{\prime} + \varepsilon t^{1/3} \frac {\mathrm{d}}{\mathrm{d}\zeta }\left ( f \frac {\mathrm{d} f}{\mathrm{d} \zeta }\right )\!, \end{equation}

which is the perturbed version of (2.10). This solution is subject to the boundary conditions

(4.18a) \begin{align} \frac {\partial }{\partial \zeta }\left (f h^{\prime}\right ) = \varepsilon t^{1/3} \quad &\mathrm{at} \quad \zeta = 0, \end{align}

(4.18b) \begin{align} h^{\prime} = - \zeta '_N t^{1/3} \frac {\mathrm{d}f}{\mathrm{d}\zeta } \quad &\mathrm{at} \quad \zeta = \overline {\zeta }_N , \end{align}

(4.18c) \begin{align} \frac {\partial h^{\prime}}{\partial \zeta } = - \varepsilon t^{1/3} \frac {\mathrm{d}f}{\mathrm{d}\zeta } \quad &\mathrm{at} \quad \zeta = \overline {\zeta }_N. \end{align}

4.4. Early time regime, $t \ll 1$

At early times, $t\ll 1$ , the spreading in the $y$ -direction is negligible compared to the flow along permeability streaks in the $x$ -direction (figure 2), and we observe deviations from the mean height throughout the current. We will show that these perturbations to the height and nose grow at the same rate as the respective mean quantities. We do so using two different approaches. First, we use the perturbed evolution (4.17) derived above, neglecting the spreading in the $y$ -direction a priori. Second, since $\kappa =\kappa (y)$ and there is negligible spreading in $y$ , we use the solution to the homogeneous permeability problem (4.3) scaled by the local permeability.

When $t\ll 1$ , the lateral spreading is insignificant, therefore we can neglect the term $t^{4/3}fh^{\prime}$ in (4.17). To balance the remaining terms, we take $h^{\prime}\sim t^{1/3}$ so that

(4.19) \begin{equation} h^{\prime} = t^{1/3}\, H(\zeta ). \end{equation}

The boundary condition at $\overline {\zeta }_N$ , (4.18b ), indicates that in similarity space, the nose perturbation is independent of time, so that $\zeta '_N$ is a constant to be determined. Substituting (4.19) into (4.17), we now have a modified equation to solve:

(4.20) \begin{equation} \frac {1}{3}H = \frac {2}{3}\zeta \frac {\mathrm{d} H}{\mathrm{d} \zeta } + \frac {\mathrm{d}^2}{\mathrm{d} \zeta ^2}(f H) + \varepsilon \frac {\mathrm{d}}{\mathrm{d}\zeta }\left ( f \frac {\mathrm{d} f}{\mathrm{d} \zeta }\right ). \end{equation}

This evolution equation is subject to

(4.21a) \begin{align} \frac {\mathrm{d}}{\mathrm{d}\zeta }\left (f H\right ) &= \varepsilon \quad \mathrm{at} \quad \zeta = 0, \end{align}

(4.21b) \begin{align} H & = - \zeta '_N \frac {\mathrm{d}f}{\mathrm{d}\zeta } \quad \mathrm{at} \quad \zeta = \overline {\zeta }_N, \end{align}

(4.21c) \begin{align} \frac {\mathrm{d} H}{\mathrm{d} \zeta } &= -\varepsilon \frac {\mathrm{d}f}{\mathrm{d}\zeta } \quad \mathrm{at} \quad \zeta = \overline {\zeta }_N. \end{align}

Equation (4.20) subject to (4.21a ) and (4.21c ) provides the solution for $H$ , hence we can find the deviation to the nose position $\zeta '_N$ using (4.21b ). We solve this numerically for a range of $\varepsilon$ values, and find good agreement with simulations, as shown in figure 3 for $\varepsilon =0.1$ . Using (4.21b ) gives value $\zeta '_N=0.49\varepsilon$ . Translating back to physical space, this means that the perturbations to the mean nose position grow according to

(4.22) \begin{equation} x^{\prime}_N = 0.49\varepsilon t^{2/3}, \end{equation}

at early times. That is, the amplitude of the nose perturbations is half the amplitude of the permeability perturbations.

Figure 3.

Height and height perturbations scaled by $t^{-1/3}$ in similarity space. Simulations shown with symbols and the theoretical result from solving (4.20) are shown with the black solid line. Heights in the simulation are evaluated at midpoints of grid cells, and the grid spacing at injection increases with time, which creates the small gap in simulation data towards $\zeta =0$ for $t=0.35$ .

The power law behaviours of both the height and nose perturbations are the same as the mean behaviours. At early times, the buoyant fluid advances at a speed set by the local permeability. Together, the observation of the power-law growth and sensitivity set by the local permeability motivate a different approach, which we now outline. We can write the full height profile as a scaled version of the uniform permeability solution (4.3):

(4.23a) \begin{align} h &= \kappa ^{-1/3} t^{1/3} f\left ( \kappa ^{1/3}\zeta \right )\!, \end{align}

(4.23b) \begin{align} \zeta _N &= \kappa ^{1/3} \overline {\zeta }_N, \end{align}

where the powers for the permeability are determined by the scaling (2.8). Fluid flowing along a high permeability region travels further and has a lower height at injection compared with a low permeability region at the same time, as seen in figure 2.

Comparing this with our expression for the nose perturbation (4.8c ) in physical space, this can be expressed as

(4.24a) \begin{align} x^{\prime}_N \cos y &= \overline {\zeta }_N t^{2/3} (1-\kappa ^{1/3}) \nonumber \\ &= \overline {\zeta }_N t^{2/3}\left ( 1 - \left [1 + \dfrac {1}{3}\varepsilon \cos (2\pi y) + \text{h.o.t.} \right ] \right )\!, \end{align}

(4.24b) \begin{align} \Rightarrow\quad x^{\prime}_N &= \dfrac {1}{3}\varepsilon \overline {\zeta }_N t^{2/3} \approx 0.49 \varepsilon t^{2/3}, \end{align}

in agreement with the result of the perturbation analysis (4.22). This result for $x^{\prime}_N$ matches very well with simulations shown in figure 4 for $\varepsilon \leqslant 0.75$ .

Figure 4.

Evolution of nose perturbation for a range of $\varepsilon$ values. Symbols show simulations, calculated as $( \max {x_N} - \min {x_N} )/2$ . Solid lines show theory of early behaviour (4.22) and late behaviour (4.43).

As we observed in figure 2, this early solution starts to break down as lateral spreading begins to dominate. figure 3 shows the start of this breakdown, where the later height perturbation profiles do not collapse onto the solution of (4.20). The transition between the two regimes is governed by the distance between a permeability minimum and maximum, which here is $L=1/2$ . To quantify when this transition occurs, we balance lateral spreading with time derivative in (2.10) to find that $T\sim L^{3/2}$ , hence the transition time is

(4.25) \begin{equation} \tau =\left (\frac {1}{2} \right )^{3/2}\approx 0.35. \end{equation}

4.5. Late time regime, $t \gg 1$

As the current increases in size, the velocity of the front slows down. Once the current slows, transverse spreading becomes more significant, and the current transitions to a regime where spreading in the $y$ -direction smooths out the thickness of the current in the bulk, and the permeability heterogeneity is felt most strongly at the nose. The main body of the current is close to the mean behaviour. Towards the nose, the deviations from the mean are more significant, as seen in figure 2. Here, we use the late time asymptotics to examine where the dynamics of the current is sensitive to permeability variations and with what magnitude. A balance is achieved between lateral spreading and forcing of the flow by gradients in the permeability, which gives a solution for the interior of the current. To describe the behaviour towards the mean nose position, we introduce a boundary layer description. As part of this, we show that the nose perturbation tends to a constant value, which we calculate.

First, we examine the interior region away from the nose. To distinguish between the two regions, we express the height perturbations in the interior as $h_{i} = h^{\prime}$ . At late times, the lateral spreading and permeability forcing are the dominant terms in (4.17); balancing the two gives a scaling $h_{i} \sim t^{-1}$ . This motivates the transformation

(4.26) \begin{equation} h_{i} (\zeta ,t) = H(\zeta )\, t^{-1}, \end{equation}

which removes all the derivatives to the height perturbation $H$ in (4.17), and the relationship may be simply written as

(4.27) \begin{equation} H = \frac {\varepsilon }{4\pi ^2 f} \frac {\mathrm{d}}{\mathrm{d}\zeta }\left ( f \frac {\mathrm{d} f}{\mathrm{d}\zeta }\right )\!. \end{equation}

When we compare this with the simulations, we find that (4.27) is an underestimate, and the relative error approaches a constant. To capture the behaviour more accurately, and to compute the approach to this late time behaviour, we propose the ansatz

(4.28) \begin{equation} h_i = \alpha (t)\, \frac {1}{f} \frac {\mathrm{d}}{\mathrm{d}\zeta }\left ( f \frac {\mathrm{d} f}{\mathrm{d}\zeta }\right )\!, \end{equation}

and fit for the prefactor $\alpha (t)$ . The result is shown in figure 5(a), where we see that there is very good agreement for the functional form. The ratio of the expected prefactor ${\varepsilon }/({4\pi ^2 t})$ to the fitted one $\alpha (t)$ goes to a constant value as time increases, as shown in figure 5(b). This shows that the $1/t$ behaviour from (4.26) is correct, but with an additional constant factor $0.571$ . The difference may be explained by including higher-order terms, but the precise nature is left for future work. As the solution for the interior approaches the mean nose position, $f$ is approximately linear, as given in (4.7), so the height decays with distance from the mean nose position, $H\sim (\overline {\zeta }_N - \zeta )^{-1}$ .

Figure 5.

(a) Height deviations back from the mean nose position. Simulations are shown with symbols, theory for the nose (4.41) and interior (4.28) with solid lines, and the composite solution (4.42) with dashed lines. (b) Ratio of the expected prefactor from (4.27) to the one fit to simulations, $\alpha (t)$ .

Therefore, this solution breaks down towards the mean nose position, where the perturbation to the permeability is comparable with the size of the current, and we must introduce a boundary layer at $\overline {\zeta }_N$ . To understand what happens at the nose, we change coordinates to a frame of reference travelling with the mean nose position, $\eta = \overline {\zeta }_N- \zeta$ . The height may now be written as

(4.29) \begin{equation} h^{\prime} = h^{\prime}(\eta ,t). \end{equation}

This transforms (4.17) to

(4.30) \begin{equation} t\frac {\partial h^{\prime}}{\partial t} = \frac {2}{3}(\eta - \overline {\zeta }_N) \frac {\partial h^{\prime}}{\partial \eta } + \frac {\partial ^2}{\partial \eta ^2}(f h^{\prime}) - 4\pi ^2 t^{4/3}fh^{\prime} + \varepsilon t^{1/3} \frac {\mathrm{d}}{\mathrm{d}\eta }\left ( f \frac {\mathrm{d} f}{\mathrm{d} \eta }\right )\!. \end{equation}

In the interior, we find that $h^{\prime}\sim t^{-1}\eta ^{-1}$ , meaning that the dominant terms grow according to $t^{1/3}$ . We also know that towards the mean nose position, the mean height is linear, $f\sim \eta$ . As $\eta \rightarrow 0$ , the advective term grows fastest, $h^{\prime}_\eta \sim t^{-1}\eta ^{-2}$ . This is comparable to the interior balance when $t^{-1}\eta ^{-2}\sim t^{1/3}$ , or

(4.31) \begin{equation} \eta \sim t^{-2/3}. \end{equation}

Motivated by this scaling of the boundary layer at the nose, we take $\delta =t^{-2/3}$ as the size of the boundary layer, and close to the nose define a local coordinate $\xi =\eta /\delta (t)$ and let $h_{N}=h_{N}(\xi ,t)$ be our solution for $h^{\prime}$ near the nose. Since we are close to the mean nose position, it is reasonable to take the linear approximation for the mean height $f=({2}/{3\delta })\overline {\zeta }_N\xi$ as in (4.7). These changes of variable transform (4.30) into

(4.32) \begin{align} t\left (\frac {\partial h_N}{\partial t} +\frac {2\xi }{3t}\frac {\partial h_N}{\partial \xi }\right ) &= \frac {2}{3}\frac {1}{\delta }\left (\delta \xi -\overline {\zeta }_N\right ) \frac {\partial h_N}{\partial \xi } + \frac {1}{\delta }\frac {2}{3}\overline {\zeta }_N\frac {\partial ^2}{\partial \xi ^2}(\xi h_N) - 4\pi ^2 t^{4/3}h_N \nonumber \\&\quad + \frac {4}{9}\varepsilon \overline {\zeta }^2_N t^{1/3}, \end{align}

subject to

(4.33a) \begin{align} h_N & =\frac {2}{3}\overline {\zeta }_N \zeta '_N t^{1/3} \quad \kern7pt \mathrm{at} \quad \xi =0, \end{align}

(4.33b) \begin{align} \frac {\partial h_N}{\partial \xi } & = - \frac {2}{3}\varepsilon \overline {\zeta }_N t^{-1/3} \quad \mathrm{at} \quad \xi =0, \end{align}

at the mean nose position. Since $t\gg 1$ , the latter can be approximated as

(4.34) \begin{equation} \frac {\partial h_N}{\partial \xi } \approx 0 \quad \mathrm{at} \quad \xi =0. \end{equation}

We replace the boundary condition at the origin with a matching condition. The interior and nose heights must have the same power-law dependence on time when $\eta \sim t^{-2/3}$ , so we impose

(4.35) \begin{equation} h_N \sim t^{-1/3}. \end{equation}

Finally, the nose perturbation must be continuous, so the late time nose perturbation must be equal to the early time behaviour (4.22) at the transition time,

(4.36) \begin{equation} \zeta '_N=\frac {1}{3}\varepsilon \overline {\zeta }_N \quad \mathrm{at} \quad t=\tau =2^{-3/2}, \end{equation}

giving us a patching condition in time. While this last condition is not strictly rigorous, we will show that it gives us a sensible solution for the height towards the nose.

The leading-order balance to $\textit {O}(\delta ^{-1})$ in (4.32) is

(4.37) \begin{equation} \frac {2}{3}\overline {\zeta }_N \frac {\partial h_N}{\partial \xi } = \frac {2}{3}\overline {\zeta }_N \frac {\partial ^2}{\partial \xi ^2}(\xi h_N), \end{equation}

which integrates to give

(4.38) \begin{equation} h_N = A(t) \log \xi + B(t), \end{equation}

where $A$ and $B$ are to be determined. The height perturbation must be finite at the mean nose position, where $\xi \rightarrow 0$ , so we set $A=0$ and find that $h_N$ is independent of the distance to the nose, consistent with (4.34). Therefore, the boundary condition (4.33a ) gives the form of the whole inner solution,

(4.39) \begin{equation} h_N=\frac {2}{3}\overline {\zeta }_N \zeta '_N t^{1/3}, \end{equation}

and we need only find the form of the nose perturbation.

To satisfy the matching power law (4.35), we must have $\zeta '_N \sim t^{-2/3}$ . Together with continuity of the nose perturbation (4.36), these give a relationship for the perturbed nose position,

(4.40) \begin{equation} \zeta '_N =\frac {1}{6}\varepsilon \overline {\zeta }_N t^{-2/3}, \end{equation}

which in turn provides a closed-form expression for the height,

(4.41) \begin{equation} h_N=\frac {1}{9}\varepsilon \overline {\zeta }^2_N t^{-1/3}. \end{equation}

The interior and nose solutions can be combined to give a composite solution for $h^{\prime}$ of the form

(4.42) \begin{equation} h^{\prime} = \big ( h_i^{-1} + h_N^{-1} \big )^{-1} = \left ( \frac { f}{\alpha (t)\,(ff_{\zeta })_{\zeta }} + \frac {9 t^{1/3}}{\varepsilon \overline {\zeta }^2_N} \right )^{-1}. \end{equation}

This is compared with simulations in figure 5. Towards the nose ( $\eta =0$ or $\zeta =\overline {\zeta }_N$ ), the theory and simulations are in very good agreement. They match less well towards the interior, but (4.42) captures the overall shape well.

Transforming back to physical space, the nose perturbation is given by

(4.43) \begin{equation} x^{\prime}_N = \frac {1}{6}\varepsilon \overline {\zeta }_N \approx 0.25 \varepsilon , \end{equation}

in good agreement with simulations (figure 4). The amplitude of the nose perturbations approaches a value that is approximately a quarter of the permeability perturbation.

4.6. Second-order effects and nonlinear coupling

So far, we have considered only the leading-order terms in our perturbation analysis, and have demonstrated good agreement with simulations of the full equations, especially for small $\varepsilon$ and at early times. The deviations of the height from the mean observed in simulations match the shape of our perturbation analysis. For the nose, the deviations from the mean agree very closely with the linearised values, even up to permeability fields varying by up to $\pm 50 \,\%$ .

The imposed single mode permeability field generates higher-order terms in the nose profile. For larger permeability variations in particular, higher-order terms can be included to more accurately capture the response of the flow to these permeability variations. Here, we briefly estimate the error between the simulations and the leading-order perturbation analysis above, showing that this error scales with $\varepsilon ^2$ . We use two typical nose profiles to show that these effects increase with increasing $\varepsilon$ , and each order has an associated wavenumber.

Figure 6.

Evolution of (a,c) $x^{\prime}_N$ and (b,d) $x^{\prime\prime}_N$ with time for permeability perturbation amplitudes (a,b) $\varepsilon =0.1$ and (c,d) $\varepsilon =0.5$ .

To calculate the size of the second-order contribution, we take the nose profiles from the simulations, subtract their mean, and subtract the linearised theory for the perturbations $x^{\prime}_N$ given by (4.22) and (4.43):

(4.44) \begin{equation} x^{\prime\prime}_N = x_N - \overline {x}_N - x^{\prime}_N \cos (2\pi y). \end{equation}

Figure 7.

Amplitude of higher-order modes forced by single mode permeability field. Solid line shows linear theory. The amplitude of the permeability perturbations is (a) $\varepsilon =0.1$ and (b) $\varepsilon =0.5$ .

The impact of the amplitude of the permeability perturbation on the shape of the profiles of $x^{\prime}_N$ and $x^{\prime\prime}_N$ is shown in figure 6. For both $\varepsilon =0.1$ and $\varepsilon =0.5$ , the amplitude of $x^{\prime\prime}_N$ is approximately a multiple of the corresponding $\varepsilon$ smaller than $x^{\prime}_N$ ; the error estimate is order $\varepsilon ^2$ . For the smaller value, $\varepsilon =0.1$ , the shape of $x^{\prime}_N$ from the simulations is close to a $\cos (2\pi y)$ curve; the size of $x^{\prime\prime}_N$ is an order of magnitude smaller than for $\varepsilon =0.5$ . For the larger value, $\varepsilon =0.5$ , $x^{\prime}_N$ pinches out in the trailing edge, a similar shape to a $\cos ^2(2\pi y)$ curve. This suggests that higher-order modes are being created in the height and nose. This effect can be readily seen by introducing a second-order correction to height $h^{\prime\prime}$ that can be included in the expression for the current height:

(4.45) \begin{equation} h(x,y,t) = \overline {h}(x,t) + h^{\prime}(x,t) \cos ({2\pi y}) + h^{\prime\prime}(x,y,t). \end{equation}

Substituting this into (2.10), the mean and first order behaviour are as above, and we find that the second-order evolution equation is

(4.46) \begin{align} \frac {\partial h^{\prime\prime}}{\partial t} - \frac {\partial ^2}{\partial x^2}\left ( \overline {h} h^{\prime\prime} \right ) - \frac {\partial ^2}{\partial y^2}\left ( \overline {h} h^{\prime\prime} \right ) ={}& \frac {1}{2} \left ( 1 + \cos ({4 \pi y})\right )\left [ \frac {\partial }{\partial x} \left ( h^{\prime} \frac {\partial h^{\prime}}{\partial x}\right ) + \varepsilon \frac {\partial ^2}{\partial x^2}\left ( \overline {h} h^{\prime} \right ) \right ] \nonumber \\ &{}- 4 \pi ^2 \cos ({4 \pi y}) \left [ \left ( \frac {\partial h^{\prime}}{\partial x} \right )^2 + \varepsilon \overline {h}h^{\prime} \right ]. \end{align}

In order to balance the $y$ dependence on the right-hand side, $h^{\prime\prime}$ must depend on $\cos ({4 \pi y})$ ; the wavenumber $n=2$ is created at second order.

Both $x^{\prime\prime}_N$ profiles are periodic in $y$ with wavelength $\lambda =1/2$ . This fits with the $\cos ^2(2\pi y )$ shape observed for $\varepsilon =0.5$ . Performing a fast Fourier transform (FFT) of the nose profiles over time shows clear peaks at the integer wavenumbers. The amplitudes of the first four modes are shown in figure 7. For $\varepsilon =0.1$ , the oscillations in figure 6 are dominated by wavenumber 2, which echoes what we see with the FFT; the amplitude of this contribution is $\sim 10^{-3}$ . The contributions from $n=3,4$ are negligible and noisy. For $\varepsilon =0.5$ , there is a clear $n=2$ repeat, but with contributions from higher-order modes to give the asymmetric shape. The FFT echoes with contributions from wavenumbers 2, 3 and 4. At each order in $\varepsilon$ , a new higher-order wavenumber is induced; order $\varepsilon ^n$ will create wavenumbers $n$ . For $\varepsilon =0.5$ , the third-order terms $\varepsilon ^3=0.125$ are comparable to those induced at leading order for $\varepsilon =0.1$ .

5.1. Response to multimodal $\kappa (y)$

We extend our previous analysis to a general Fourier series representation of the perturbed quantities. The same features as observed for a single mode still apply. At early times, the current is governed by the local permeability, and flows predominantly in the $x$ -direction. At late times, the permeability variations are felt most strongly at the nose, and lateral spreading has smoothed out their effect in the interior of the current. To leading order, our system (4.17) is linear in the perturbed quantities. Therefore, we are able to separate the modes and arrive at an analogous system for a general permeability field. We then propose the same power-law behaviour for early and late times as in §§ 4.4 and 4.5, respectively, and demonstrate that the early height and nose perturbations depend only on the corresponding permeability amplitude $\varepsilon _m$ . The transition time of each mode is governed solely by its wavelength. At late times, the height and nose variations depend on both the amplitude and wavelength of the permeability perturbation.

Echoing the variations in the permeability, we also describe the height and nose perturbations, $h^{\prime}$ and $x^{\prime}_N$ , using Fourier series in physical space

(5.2a) \begin{align} h &= \overline {h}(x,t) + \sum _{m=1}^{\infty } H_m(x,t) \cos \left (\frac {m \pi }{W} y\right )\!, \end{align}

(5.2b) \begin{align} x_N &= \overline {x}_N(t) + \sum _{m=1}^{\infty } X_m(t) \cos \left (\frac {m \pi }{W} y\right )\!, \end{align}

which have the following forms in similarity space,

(5.3a) \begin{align} h &= t^{1/3}\,f(\zeta ) + \sum _{m=1}^{\infty } H_m(\zeta ,t) \cos \left (\frac {m \pi }{W} y\right )\!, \end{align}

(5.3b) \begin{align} \zeta _N &= \overline {\zeta }_N + \sum _{m=1}^{\infty } Z_m(t) \cos \left (\frac {m \pi }{W} y\right )\!, \end{align}

where $X_m = Z_m t^{2/3}$ .

When the amplitude of the variations is small, the same process for a single mode can be repeated here to arrive at a linearised system similar to (4.17) written now in terms of the Fourier coefficients $H_m$ and $\varepsilon _m$ . To do so, we substitute the new expressions for permeability, height and nose into (2.10). The leading-order terms are single sums, and integration over the $y$ domain isolates the Fourier coefficients to produce a system to solve for a general mode $m$ :

(5.4) \begin{equation} t\frac {\partial H_m}{\partial t} = \frac {2}{3}\zeta \frac {\partial H_m}{\partial \zeta } + \frac {\partial ^2}{\partial \zeta ^2}\left (f H_m\right ) - \frac {\pi ^2 m^2}{W^2}t^{4/3} f H_m + \varepsilon _m t^{1/3}\frac {\mathrm{d}}{\mathrm{d} \zeta }\left (f\frac {\mathrm{d} f}{\mathrm{d} \zeta }\right )\!, \end{equation}

which is subject to

(5.5a) \begin{align} \frac {\partial }{\partial \zeta }\left (f H_m\right ) &= \varepsilon _m t^{1/3} \quad \mathrm{at} \quad \zeta =0, \end{align}

(5.5b) \begin{align} H_m &= -Z_m t^{1/3}\frac {\mathrm{d}f}{\mathrm{d}\zeta } \quad \mathrm{at} \quad \zeta =\overline {\zeta }_N, \end{align}

(5.5c) \begin{align} \frac {\partial H_m}{\partial \zeta } &= - \varepsilon t^{1/3}\frac {\mathrm{d}f}{\mathrm{d}\zeta } \quad \mathrm{at} \quad \zeta =\overline {\zeta }_N. \end{align}

The transition time $\tau _m$ is when the spreading in the $y$ -direction becomes comparable with the flow in the $x$ -direction. This is governed by the distance between a permeability minimum and maximum for the corresponding mode:

(5.6) \begin{equation} \tau _m = \left (\frac {\lambda _m}{2}\right )^{3/2} = \left (\frac {W}{m}\right )^{3/2}. \end{equation}

For times $t\ll \tau _m$ and for a given mode, the height scales according to

(5.7) \begin{equation} H_m = \hat {H}_m(\zeta )\, t^{1/3}, \end{equation}

where $\hat {H}_m$ solves

(5.8) \begin{equation} \frac {1}{3}\hat {H}_m = \frac {2}{3}\zeta \frac {\mathrm{d}\hat {H}_m}{\mathrm{d}\zeta } + \frac {\mathrm{d}}{\mathrm{d} \zeta }\left (f\hat {H}_m\right ) + \varepsilon _m \frac {\mathrm{d}}{\mathrm{d} \zeta }\left (f\frac {\mathrm{d} f}{\mathrm{d} \zeta }\right )\!, \end{equation}

and the nose perturbation grows in time as

(5.9) \begin{equation} X_m = \frac {1}{3}\overline {\zeta }_N \varepsilon _m t^{2/3} \approx 0.49 \varepsilon _m t^{2/3}. \end{equation}

Thus we see that for each small-amplitude mode, the early height and nose perturbations are independent of the wavelength of the corresponding permeability perturbation, and depend only on the amplitude of that mode. Hence the mode with the largest amplitude is the most significant feature initially.

When $t\gg \tau _m$ and lateral diffusion dominates, we consider the regions away from and close to the mean nose position separately to arrive at the following composite expression for the height perturbation profile:

(5.10) \begin{equation} H_m = \varepsilon _m \left ( \frac {\pi ^2 m^2 t f}{W^2 (ff_{\zeta })_{\zeta }} + \frac {9 m t^{1/3}}{2 W \overline {\zeta }^2_N} \right )^{-1}. \end{equation}

In physical space, the nose perturbation goes to a constant:

(5.11) \begin{equation} X_m = \frac {W \varepsilon _m \overline {\zeta }_N}{3m} \approx 0.25 \lambda _m \varepsilon _m. \end{equation}

The magnitude of the resultant nose perturbation increases with both increasing wavelength and amplitude of the corresponding permeability mode. Hence the dominant behaviour is not necessarily from the largest wavelength or amplitude, but is from the largest product, $\lambda _m \varepsilon _m$ .

5.2. Results

To test the validity of these linearised results, we simulate the evolution of the full system (2.10) with randomly generated multimodal permeability fields on a domain of width 4 ( $W=2$ ). For illustration, we choose to have six modes, with the largest wavelength equal to 1 and with amplitude 0.5, and the remaining five modes generated randomly. The Fourier modes $m$ are random integers taken from a uniform distribution, $m\in U(2W+1,30W)$ , to give wavenumbers $1+1/2W \leqslant n \leqslant 15$ . The amplitudes are taken alternately from two log-normal distributions, the first with mean $\log 0.1$ , the second with mean $\log 0.01$ , and both with standard deviation 1 (see figure 8). Two example simulations, runs A and B, have been chosen to compare with the linearised theory. The parameter values for these runs are shown in table 1.

Figure 8.

Probability distributions from which the $\varepsilon _m$ values are taken alternately for the simulations.

Table 1.

Permeability parameters for runs A and B, with $W=2$ .

From the linearised theory, we expect to see the nose perturbations initially grow like $t^{2/3}$ . The time at which each mode transitions to the late regime will depend on its wavelength. When $t\gt \max {\tau _m}$ , we expect $x^{\prime}_N$ to approach a constant shape. At late times, successive profiles should plot on top of each other. To compare with simulations, we examine the behaviour of the nose, both the perturbation profile over time, and its decomposition into the separate modes using an FFT. The latter we compare with the expected nose perturbation from (5.9) and (5.11). The results of this comparison can been seen in figures 9 and 10 for runs A and B, respectively.

Figure 9.

Run A. (a) The permeability profile used in this simulation. (b) The evolution of the deviations from the mean nose position from early to late times (only $0\lt y\lt 2$ is shown since permeability, and hence the nose, is symmetric about $y=0$ ). (c) The FFT of the nose profile, showing amplitude of contribution from each mode over time, with results from simulations shown with symbols, and theory with solid lines. Only modes imposed in the permeability variations are shown.

Figure 10.

Run B. (a) The permeability profile used in this simulation. (b) The evolution of the deviations from the mean nose position from early to late times (only $0\lt y\lt 2$ is shown since permeability, and hence the nose, is symmetric about $y=0$ ). (c) The FFT of the nose profile, showing amplitude of contribution from each mode over time, with results from simulations shown with symbols, and theory with solid lines. Only modes imposed in the permeability variations are shown.

For runs A and B, the dominant shape of both is from the first mode, with $\lambda _4=1$ and $\varepsilon _4=0.5$ . The shapes of the nose profiles are very similar to the single mode run with $\varepsilon =0.5$ , seen in figure 6. The amplitude of this mode agrees well with our linearised theory; the final value in the simulation is 95 % of the expected value.

The higher wavenumbers exhibit more interesting behaviour. In both runs, almost all of the individual mode amplitudes $X_m$ match well with our theory at early times, and agree with the expected transition times. At late times, the modes first approach the value that we predict, then decrease, generally to a lower steady value. This final value in the simulations is $33{-}76\,\%$ lower than predicted in (5.11) (see figure 11), with a weak dependence on wavenumber. The lower-wavenumber modes are generally closer to the expected value than higher-wavenumber modes. However, for the smallest $\varepsilon _m$ values in run B, we observe that the magnitude of these modes keeps decreasing and does not reach a steady value within the simulated time.

Figure 11.

Ratio of observed $X_{m,{sim}}$ to expected final nose perturbation value $X_{m,{theory}}$ given by (5.11). Simulation value taken as average of final five points.

The notable exception to this is the mode with wavenumber $n=3$ ( $m=12$ ) in run B, where the simulation overpredicts our theory. This overprediction is slight at early times, and large at late times; the final value is 7.6 times higher than the linearised theory. Here, the forcing from the permeability mode $n=3$ ( $m=12$ ) is small compared to the third-order term generated by the largest mode, $n=1$ . The amplitude of this first mode is large enough that even at third order, it is dominant over the forcing from the $n=3$ mode: $\varepsilon _4^3 = 0.125 \gt \varepsilon _{12}=0.017$ .

Our perturbation analysis assumed that all the $\varepsilon _m$ were of similar magnitude to arrive at a linear solution. At higher orders, or if the amplitudes of the modes $\varepsilon _m$ are orders of magnitude different, higher-order terms would need to be fully resolved in order to model the modal interactions. The amplitudes that we have used in these simulations range across two orders of magnitude, so deviations from the linear theory are expected. Interestingly, this large range of $\varepsilon _m$ values has little impact on the early time flow, where the current evolves according to its local permeability, but does matter at late times, where lateral spreading is dominant. The precise mechanics of these higher-order effects at late times is beyond the scope of this work. However, the additional variability in the permeability may allow for more lateral movement before reaching an equilibrium, late time configuration.

Throughout, we still see that the dominant mode in the permeability creates the dominant feature in the nose structure.