Response of a gravity current in a porous medium to a perturbation of the permeability field

Abstract

One of the challenges with modelling subsurface flows is the uncertainty in measurements of geological properties, mostly due to limited resolution in observation methods. Many subsurface flows can be modelled as a gravity current, which, for uniform material properties and power-law injection rate, has a well-characterised similarity solution. The similarity solution forms a dynamical attractor that is typically approached rapidly from a host of initial conditions. Here, we consider the impact of transverse variations to the permeability field by performing a perturbation analysis of the self-similar spreading. This treats the response as perturbations to the self-similar flow. We restrict our focus to permeability fields that vary laterally, or across the flow, starting with the simple case of a sinusoidal perturbation to a uniform permeability. At early times, the height and nose position of the current are determined by the local permeability, and deviations to the height and nose grow at the same rate as the mean, and proportional to the amplitude, of the permeability variation. The transition between the early and late time regimes is governed by the wavelength of the permeability. At late times, lateral spreading between high and low permeability streaks is dominant; the height deviations decay, and the nose deviations approach a steady state. The magnitudes of both depend on the product of the wavelength and amplitude of the permeability. The single mode sets the groundwork for examining more complex, multimodal permeabilities, which are more representative of real geological structures.

1. Introduction

Gravity-driven currents in porous media are observed in many geological settings and have a wide variety of environmental and industrial applications, including carbon capture and storage (CCS), geothermal energy storage, and groundwater flow. These systems are of increasing interest as part of the energy transition for their potential to decarbonise both energy generation and many industrial processes. In most industrial CCS schemes, supercritical CO $_2$ is injected deep underground, and it is important to understand the behaviour and extent of the injected CO $_2$ to ensure that it is safely trapped. In geothermal operations, heat may be stored in the summer and extracted in the winter (Dudfield & Woods 2012). For groundwater systems, the whole flow may be of interest as well as the dispersal and remediation of pollutants in freshwater aquifers. In this paper, we primarily consider gravity currents motivated by the context of geological CO $_2$ storage, but the ideas on the impact of heterogeneities are generic to many other applications.

Subsurface flows in complex porous media are often driven by buoyancy and large-scale pressure gradients. For example, injected CO $_2$ is buoyant and will rise until it reaches an impermeable caprock, when it will spread laterally. To model the flow of fluids in the subsurface, it is important to understand the structure and variability of a range of geological properties. For example, the caprock topography, thickness of the aquifer, porosity and permeability all impact how far, how fast and in which direction the fluid will travel. However, all of these properties can be highly heterogeneous on a range of length scales, often from centimetre to metre scales within a reservoir that may extend laterally for hundreds of metres to tens of kilometres. In addition, this range of geological heterogeneity often cannot be measured at a sufficiently high resolution using remote imaging technology to capture this heterogeneity. This is a particular issue for the permeability, which can vary over several orders of magnitude across very small length scales (Eiken et al. 2011). This small-scale heterogeneity can have a large effect on CO $_2$ plume propagation (Jackson & Krevor 2020) and can potentially increase plume speed by up to 200 % (Benham et al. 2021a ). As such, we seek to understand the impact of heterogeneity in the form of spatially dependent variations to a uniform permeability field.

We use a deterministic approach for our model by imposing a permeability structure. However, there is a large body of work on the modelling of groundwater flows using stochastic permeability fields (Freeze 1975; Dagan 1984; Rubin 2003). In this approach, subsurface properties are treated as random variables, sometimes with spatially dependent properties, leading to a statistical description of the flow.

Since the flow of CO $_2$ is driven by buoyancy, and has a large aspect ratio, it is often modelled as a gravity current (Szulczewski et al. 2012; Huppert & Neufeld 2014). With uniform material properties, buoyancy driven porous gravity currents resulting from both axisymmetric and line injections admit similarity solutions (Huppert & Woods 1995; Lyle et al. 2005), giving simple power-law behaviours in time for the height and extent of the current. Previous work has modelled large-scale heterogeneity in gravity currents using simple relationships with distance. For example, Ciriello et al. (2013) analysed power-law variations with distance or height from the point of injection, and Hinton & Woods (2018) demonstrated that a linear vertical gradient in permeability can, in some cases, result in shocks and rarefaction waves.

We approach modelling heterogeneity as a perturbation to a homogeneous permeability field. This induces a corresponding perturbation to the height and nose position of the current; we split these into a base state and a perturbed component. The base state is a well-understood solution to the homogeneous gravity current equations. However, this base state is time dependent. To remove the time dependence from the base state, Mathunjwa & Hogg (2006) instead considered the perturbations in similarity space, and performed the perturbation there, demonstrating that the similarity solution is stable to perturbations in a homogeneous medium. We adopt a similar approach by exploiting this feature of similarity space.

We begin, in § 2, by re-deriving the gravity current equations where the permeability field is subject to spatial variations. We non-dimensionalise our model using the characteristic wavelength of the permeability heterogeneity, and solve this system numerically in § 3, building on code first presented in Cowton et al. (2018). We present the results of a typical simulation, in which the current demonstrates distinct early and late time behaviours for both the height and nose position of the current. To understand the behaviour further, in § 4 we briefly revisit gravity currents in homogeneous strata, and then perform a perturbation analysis in the self-similar coordinates to understand how permeability heterogeneities alter the structure of the homogeneous case. In § 5, we extend our analysis to more complex permeability fields, which we create using randomly generated, multimodal functions. We see that these perturbation solutions are additive, and observe some interesting features from modal interactions. Finally, in § 6, we close by re-dimensionalising some key parameters, and discuss the results of our analysis in the context of CCS.

2. Governing equations

We consider the injection of CO $_2$ of density $\rho$ along a line at $x=0$ into an unconfined aquifer saturated with water of density $\rho + \Delta \rho$ , where $x$ is the direction of spreading, $y$ is the transverse direction, and $z$ points vertically down, in the same direction as gravity $g$ (see the schematic in figure 1). The two-dimensional flow associated with injection along a line is considered here mainly because it simplifies much of the analysis. Nevertheless, in many settings we anticipate that locally, the edge of the CO $_2$ would only vary on a long wavelength if the properties were spatially uniform. Hence the resulting analysis may inform the sensitivity of the spreading plume to geological heterogeneity. Compared with the ambient brine, CO $_2$ is buoyant, $\Delta \rho \gt 0$ , and therefore rises through the subsurface until it reaches an impermeable barrier at $z=0$ . Due to the action of buoyancy forces, the current becomes long and thin. A reservoir can be considered unconfined if the thickness of the CO $_2$ does not exceed approximately 10 % of the reservoir thickness (Pegler, Huppert & Neufeld 2014). Under this assumption, the ambient fluid is unaffected by the injected CO $_2$ , hence its flow can be neglected.

Figure 1. Schematic showing line injection at $x=0$ of buoyant fluid into unconfined ambient.

The Darcy flux of $\textrm {CO}_{2}$ , $\boldsymbol{u}(\boldsymbol{x},t)$ , must satisfy both conservation of mass and Darcy’s law, i.e.

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

(2.1b) \begin{align} \boldsymbol{u} &= -\frac {k(\boldsymbol{x})}{\mu } ( {\boldsymbol{\nabla }} p - \rho g \hat {\boldsymbol{z}}), \end{align}

where $k$ is a scalar permeability, $\mu$ is the dynamic viscosity of CO $_2$ , and $p$ is the pressure. Consistent with our gravity current analysis, we consider a depth-integrated permeability, and to be amenable to the lubrication approximation, consider only slow horizontal variations. In general, the permeability and porosity will both be functions of space, $k(\boldsymbol{x})$ and $\phi (\boldsymbol{x})$ , and importantly will include vertical structure. This could be due to layered sedimentary sequences or graded bedding, for example. The main effect on the flow will be to introduce vertical gradients in the horizontal velocity, resulting in somewhat different depth-integrated models of flow (Hinton & Woods 2018; Benham et al. 2021b ). Additionally, such vertical permeability structures and the associated vertical variations in the horizontal velocity are likely to enhance or promote mixing (in the case of miscible flows, see Nijjer, Hewitt & Neufeld 2019), or lead to complex saturation patterns (see Krevor et al. 2011). For the sake of simplicity, here we neglect such vertical variations of the permeability field, and consider only lateral variations in the permeability field, $k(\boldsymbol{x})=k(y)$ . Although porosity typically correlates with permeability, here we assume that it is constant, again for simplicity.

Since the flow within the long, thin current is predominantly horizontal, the pressures within both the CO $_2$ and ambient fluid are approximately hydrostatic. Within the ambient, the pressure is therefore as $p_a(z) = p_0 + (\rho +\Delta \rho )gz$ , where $p_0$ is a constant reference pressure, when $h$ is much less than the depth of the reservoir (the unconfined limit). Similarly, in the CO $_2$ , $p(\boldsymbol{x},t)=p_1(\boldsymbol{x},t) + \rho g z$ , where $p_1$ is the pressure at $z=0$ . These pressures must be equal at the interface $z=h(x,y,t)$ , giving an expression for $p_1$ in terms of the height and reference pressure, and therefore

(2.2) \begin{equation} p=p_0 + (\rho + \Delta \rho ) g h + \rho g (z-h), \quad 0\lt z\lt h. \end{equation}

Hence the pressure gradients driving the flow can be expressed in terms of gradients in the thickness of the current. Because of the large aspect ratio of the flow, the vertical component of velocity within a gravity current is small compared with the horizontal, so it is reasonable to vertically average the Darcy velocity and the mass conservation equation. The former gives us, by definition, the depth-integrated fluid flux:

(2.3) \begin{equation} \boldsymbol{q}(\boldsymbol{x},t) = - \frac {\Delta \rho \textit{gk}(\boldsymbol{x})}{\mu } h {\boldsymbol{\nabla }} h. \end{equation}

The latter gives a relationship between the evolution of height with time and the divergence of fluid flux. Combining (2.3) with the vertically averaged mass conservation provides a single partial differential equation describing the evolution of the thickness of the current:

(2.4) \begin{equation} \phi \frac {\partial h}{\partial t} = \frac {\Delta \rho g}{\mu }\, {\boldsymbol{\nabla }} \boldsymbol{\cdot } \left [ k(\boldsymbol{x}) h {\boldsymbol{\nabla }} h \right ]. \end{equation}

To complete the description of the system, we impose a constant injection flux $q$ across the domain $-W\lt y\lt W$ at the origin, $x=0$ , vanishing height and zero flux at the nose:

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

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

(2.5c) \begin{align} \boldsymbol {q}\boldsymbol{\cdot }\hat {\boldsymbol {n}} &= 0 \quad \mathrm{at} \quad x=x_N, \end{align}

where $\hat {\boldsymbol{x}}$ is the unit vector in the $x$ -direction, $x_N=x_N(y,t)$ is the nose of the current, and the normal $\hat {\boldsymbol {n}}$ to the nose position is given by

(2.6) \begin{equation} \hat {\boldsymbol {n}} = \frac {(1,x_{N,y})}{\sqrt {(1+x_{N,y}^2)}}, \end{equation}

where $x_{N,y}=\partial x_N/ \partial y$ .

In many geological settings, the depositional environment is either known or can be inferred, which gives a dominant length scale over which the permeability varies. Here, we will treat this length scale as a wavelength $\lambda$ , and consider permeability fields that are a combination of a mean and a small sinusoidal perturbation. We will use this dominant perturbation wavelength to non-dimensionalise our governing system of equations, and note that for more complex permeability fields, a pragmatic choice of non-dimensionalisation may be used.

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.

3. Simulations

We solve the dimensionless gravity current (2.10), subject to a constant injection along $x=0$ and no flux at all other edges of the computational domain (the boundary conditions (2.11) will be used in our subsequent analysis). We use a numerical scheme similar to that presented in Cowton et al. (2018) and Gilmore (2022). This scheme is written in a fixed volume, flux conservative manner, using a Crank–Nicolson discretisation for the flux. A predictor–corrector step is used to deal with the nonlinear diffusivity. An alternating direction implicit scheme is used to advance in $x$ , then $y$ (Darms et al. 2002). We use adaptive time stepping. At each time step, we compare the result of performing one time step with two half time steps, and compute the L2 norm, scaled by the maximum height and number of grid cells. If this measure of the two heights is below our error threshold $10^{-4}$ , then the time step is either multiplied by 1.5 or set by the estimated speed of the nose, whichever is smaller. If the error is more than twice the threshold, then it is divided by 1.5. We save the full output at approximately logarithmically spaced intervals to capture the power-law behaviour of the current.

The grid spacing in $x$ varies with distance and time in a prescribed manner so as to maintain high resolution at the nose of the current in order to fully resolve the perturbations. The regridding process is triggered when the leading edge of the current occupies over $90\,\%$ of the grid cells. The grid cells from injection to the trailing edge of the nose, minus a buffer number of cells, are then averaged, and the corresponding grid spacing is doubled. A convergence test showed that for a uniform permeability field, the simulations converged for ${\rm d}x_{min}\leqslant 2^{-8}$ and ${\rm d}y\leqslant 0.5$ . The simulations presented here have ${\rm d}y\leqslant 0.125$ and either ${\rm d}x_{min}=2^{-10}$ or ${\rm d}x_{min}=2^{-11}$ .

We solve (2.10) for exemplar permeability fields of the form

(3.1) \begin{equation} \kappa = 1 + \varepsilon \cos (2\pi y + \pi ), \end{equation}

where $\varepsilon$ is an imposed amplitude. The phase shift ensures that the permeability is at its minimum at the edges of the computational domain, but it is not important to the form of the solution or the subsequent analysis, so we will exclude it from our theoretical analysis.

The result of a typical simulation is presented in figure 2. The thickness of the current decreases with distance from the origin and forms a nose whose shape echoes the sinusoidal shape of the permeability field. Comparing the height with its mean in the $y$ -direction, we see that at early times, the height deviates from the mean throughout the current. However, at late times, these deviations are strongest at the origin and towards the nose.

Figure 2. Simulation results for $\varepsilon =0.1$ . The two columns show the same simulations at two different times: (a,c,e) $t=0.215$ and (b,d,f) $t=62.4$ . Graphs show (a,b) the height, $h$ , (c,d) the deviations from the transverse mean height, $h-\overline {h}$ , and (e,f) three slices of the full height: the minimum (orange circles), transverse mean (magenta crosses) and maximum (green plusses). The minimum and maximum are taken from the corresponding orange and green lines in (a,b). An inset shows the deviation of the minimum and maximum from the mean. The black solid line is the self-similar solution (4.3b ), which the simulation mean matches.

In general, we find that the behaviour of the current depends on the relative magnitude of the lateral and longitudinal spreading. At early times, lateral spreading is small, and the current flows predominantly in the $x$ -direction, and each line of the current propagates independently, governed by its local permeability. The injection is uniform across the whole interval, so initially the current contains approximately the same volume of fluid per unit width in the $y$ -direction. The current moves faster along lines of higher permeability: at a given time, it will have travelled further and hence will have a lower height at the injection point. The opposite is true along lines of lower permeability: the current will travel more slowly, so go less far and have a higher height at the injection point. This means that the height varies with respect to the mean profile throughout the current. The form of the permeability perturbation controls the local perturbation to the height and nose position, therefore perturbations of the current thickness will have the same functional form as the permeability field.

When lateral spreading is significant at late times, the effect of the variations in permeability is reduced in the bulk of the current. Instead, the variations in permeability are felt most strongly at the origin and nose of the current. At the origin of the current, freshly injected fluid has not had enough time to spread in the $y$ -direction, and so is significantly impacted by permeability variations. At the nose, the current is thinner, and its height is comparable with the wavelength of the permeability perturbations. We found that at late times in our simulations, the height deviates from the mean most strongly at the origin and nose of the current.

4. Influence of permeability variations

We now explore analytically the implications of perturbing the permeability field. Motivated by the distinct early and late time behaviours that we observe in the numerical simulations, we decompose the permeability, height and nose position into mean and perturbed components. We start in § 4.1 by considering a uniform permeability field to give us the form of the mean solution about which we will consider variations. This mean solution is self-similar, a feature that we take advantage of when considering the influence of perturbations on the underlying permeability field. Then a perturbation to the permeability field is introduced in § 4.2, and we derive a system of equations to describe the resulting height and nose perturbations. In § 4.3, we transform the perturbation equations to self-similar coordinates, so that we can proceed with our analysis using a time-independent base state. In § 4.4, we look at the early time behaviour, when flow is predominantly in the $x$ -direction. In § 4.5, we examine the late time behaviour, when spreading in the $y$ -direction dominates. Finally, we briefly estimate the second-order effects and highlight the nonlinear coupling between modes, seen particularly with higher amplitude permeability forcing.

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 (1995) 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. Generalised permeability fields

The analysis in § 4 can readily be extended to more general and more realistic permeability fields, such as those formed in complex depositional environments. Our analysis is still focused on permeability fields with transverse permeability streaks, which we now represent by a Fourier series

(5.1) \begin{equation} \kappa (y) = 1 + \sum _{m=1}^{\infty } \varepsilon _m \cos \left (\frac {m \pi }{W} y\right )\!, \end{equation}

with amplitudes $\varepsilon _m$ , wavelengths $\lambda _m = 2W/m$ , and wavenumbers $n_m = m/2W$ . For simplicity, we initially consider the case where all $\varepsilon _m\ll 1$ are of the same order. The single mode analysis in § 4 forms a basis for our analysis here. We use the same methods, and show that each mode is independent, to leading order, and the solution for each is very similar to the single mode case but with the additional dependence on wavelength. Simulations of two randomly generated multimodal permeability fields are presented in § 5.2, with amplitudes of the permeability modes spanning two orders of magnitude. In these, we see that our analysis matches well at early times, and for the largest mode at late times too. The higher modes exhibit interesting behaviour at late times, which we discuss.

Considering only the leading-order perturbation for each mode is still able to capture the leading-order impact of the permeability variations. To capture the full behaviour, the addition of these modes is not sufficient, and higher-order nonlinear terms created by the interaction of different modes ought to be considered. Accounting for weakly nonlinear mode coupling has been studied in the context of Hele-Shaw cell models, e.g. by Miranda & Widom (1998). However, given the large order of magnitude variations seen in geological settings, weakly nonlinear models are unlikely to be able to capture the leading-order effects. Instead, we anticipate that the identification of dominant modes, as outlined in the following analysis, is likely to be the most fruitful application of this analysis.

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.

6. Discussion and conclusion

In this analysis, we have sought to mathematically illuminate the impact of permeability variations on the observed variance in the spreading of buoyancy-driving flows in porous media. For a gravity current flowing through a permeability field with a sinusoidal perturbation, we find that the height and nose perturbations initially grow at the same rate as the mean, leading to a locally self-similar profile (here in $y$ ). At late times, the height perturbations are localised around the mean nose position, and the nose perturbations approach a steady state. The transition time between these two states is governed by the wavelength of the permeability variations, with smaller wavelengths transitioning more rapidly. For modest permeability forcing, these solutions are additive, provided that the amplitudes of the permeability modes are comparable. For these multimodal permeabilities, the height and nose perturbations exhibit behaviour similar to that observed for a single mode. The separability of these solutions breaks down with larger permeability amplitude variations. When this happens, the late time behaviour is more complex. The largest product of wavelength and amplitude is the dominant feature, and the smaller modes are suppressed.

To consider these results in the context of geological carbon storage, we use the parameter value ranges for different storage sites collated by Benham, Neufeld & Woods (2022) and shown in table 2. For geological features with length scale 10–100 m, the transition time $\tau$ is between 4.8 hours and 280 days, which is short compared with the predicted lifetime of an industrial CCS project, where the injection phase is expected to last approximately 20 years (Heddle, Herzog & Klett 2003).

Table 2. Typical reservoir values from Benham et al. (2022).

Carbon storage sites are monitored using time lapse seismic imaging, which has typical resolution approximately 10 m, corresponding to the quarter wavelength at common storage depths. The final nose perturbation, dimensionally, is $x^{\prime}_N = 0.25\varepsilon L$ , i.e. at most a quarter of the length scale of variations. Therefore, the smallest permeability variations that can be observed are approximately 40 m. Any smaller features would likely be below the resolution of the seismic imaging. Nevertheless, even when it cannot be resolved seismically, the interaction of the buoyant CO $_2$ current with geological heterogeneities across many wavelengths will lead to a highly convoluted front. This increases the contact area between CO $_2$ and brine, which also potentially enhances the rate at which CO $_2$ dissolves, and is thereby stably trapped. A consequence of this analysis is that the magnitude of this trapping by dissolution is to some extent mitigated by buoyant spreading, particularly at late times when the buoyancy of the current limits the extent to which heterogeneities drive fingering.

Prior to injection, data on reservoir heterogeneity are similarly horizontal. While detailed information is often available from the injection well itself, in practice, seismic data are used to interpolate further into the reservoir. The present analysis suggests that uncertainty in geological heterogeneity may be linked to uncertainty in plume location, but that this uncertainty first increases and then saturates (early and late time, respectively, in our analysis). Importantly, this also suggests that if the transition in plume morphology can be observed, e.g. seismically, then this may be used to further infer the wavelength and magnitude of the dominant underlying heterogeneity.

Throughout the analysis, we have assumed that the permeability variations are close to the mean. In reality, we often see very large changes of permeability. We observe pinching out of the trailing edge of the nose when $\varepsilon \geqslant 0.5$ . Even for the larger amplitudes of permeability perturbation, the nose reaches a steady shape. This kind of periodic permeability does not lead to continual growth. Most importantly, our multimodal analysis suggests that it is the modes that have the greatest product of permeability perturbation and wavelength that will dominate the variance in the flow. Hence very localised high-permeability streaks are expected to have less effect than broad, larger-wavelength variations.

In stochastic models, permeability (or more generally, hydraulic conductivity) is typically modelled as a log-normal distribution (Rubin 2003) with variance $\sigma _Y^2$ . Under this assumption, Dagan (1984) showed, with a uniform flow in the $x$ -direction, that the longitudinal dispersion coefficient and lateral displacement coefficient tend towards constants after the flow has travelled tens of the characteristic length scale. While not directly analogous to the work presented here, both demonstrate a maximum effect that permeability variations can have. An important extension of this model, motivated in part by stochastic models, but also from field observations of permeability heterogeneity, is to explore the role of longitudinal and lateral variations in the permeability field.

There are many complicating factors that we have neglected for simplicity. Periodicity of the permeability used allows us to consider a global uncertainty to a permeability, but it does not let us say much about a single high or low permeability streak, for example. Topography, both variability in caprock and overall slope, will also have a large impact on plume speed and location. We leave these intriguing complexities for future study.