Buoyancy segregation suppresses viscous fingering in horizontal displacements in a porous layer

We consider the axisymmetric displacement of an ambient fluid by a second input fluid of lower density and lower viscosity in a horizontal porous layer. If the two fluids have been vertically segregated by buoyancy, the flow becomes self-similar with the input fluid preferentially flowing near the upper boundary. We show that this axisymmetric self-similar flow is stable to angular-dependent perturbations for any viscosity ratio. The Saffman-Taylor instability is suppressed due to the buoyancy segregation of the fluids. The radial extent of the segregated current is inversely proportional to the viscosity ratio. This horizontal extension of the intrusion eliminates the discontinuity in the pressure gradient between the fluids associated with the viscosity contrast. Hence, at late times viscous fingering is shut down even for arbitrarily small density differences. The stability is confirmed through numerical integration of a coupled problem for the interface shape and the pressure gradient, and through complementary asymptotic analysis, which predicts the decay rate for each mode. The results are extended to anisotropic and vertically heterogeneous layers. The interface may have steep shock-like regions but the flow is always stable when the fluids have been segregated by buoyancy, as in a uniform layer.


Introduction
Displacement flows in porous media occur in many important environmental settings including geological CO 2 sequestration, underground hydrogen storage, geothermal energy generation and the infiltration of sea water into fresh water aquifers. If the displacing fluid is of lesser viscosity than the ambient or displaced fluid, then the fluid-fluid interface can become unstable. This phenomenon is known as the 'Saffman-Taylor instability' or as 'viscous fingering' owing to the fingers of low viscosity input fluid that penetrate the ambient fluid (Saffman & Taylor 1958;Paterson 1981;Homsy 1987). In many applications, understanding this instability is a key concern. For example, sequestered carbon dioxide is much less viscous than the host brine in the aquifer and viscous fingering may significantly reduce the fraction of the aquifer that can be accessed by the CO 2 , which in turn reduces the storage efficiency (Bachu 2015). Multiple strategies have been developed to suppress viscous fingering including controlling the flow geometry and varying the volume flux of the input fluid (Al-Housseiny et al. 2012;Zheng et al. 2015a). Significant miscibility between the fluids has also been shown to shut down fingering at later times (Nĳjer et al. 2018;Sharma et al. 2020).
In the case that the two fluids have different densities, as is relevant to CO 2 sequestration, variations in the hydrostatic pressure provide an extra horizontal force that drives (or opposes) the fluid displacement; such flows are known as 'gravity currents' (Huppert & Woods 1995). There has been much research into gravity currents in porous media, particularly in confined layers in which the ambient fluid must be displaced (Hesse et al. 2007;Juanes et al. 2010;Zheng & Stone 2022). Viscous fingering can occur at early times if the input fluid is of lower viscosity. However, it is generally assumed that at later times a stable interface between the input and ambient fluids develops, which seems to be corroborated by laboratory experiments (Nordbotten & Celia 2006;Pegler et al. 2014;Guo et al. 2016). It has also been shown that the combination of buoyancy and mixing can suppress the Saffman-Taylor instability at long times (Tchelepi 1994;Riaz & Meiburg 2003;Nĳjer et al. 2022). However, for sharp-interface gravity currents in confined porous media, stability of the gravity current solution has not yet been confirmed and the flow physics that suppresses viscous fingering is not fully understood.
In the present paper, we consider injection of a relatively buoyant fluid into an aquifer confined above and below by two impermeable layers. The injection source is on the upper boundary and the aquifer is initially filled with a second fluid of greater viscosity (see figure 1). We show that if the fluids have segregated owing to buoyancy, the gravity current solution is stable to both axisymmetric and angular-dependent perturbations for any values of the viscosity ratio and the input flux relative to the buoyancy velocity.
The classical Saffman-Taylor instability is associated with the discontinuity in the pressure gradient across the fluid-fluid interface (Paterson 1981). For example, in the case of flow in a porous medium with no vertical dimension, the base radial displacement has a circular interface. The total radial flux in each fluid in a circular cross-section is given by where is the permeability, is the viscosity, is the radial coordinate and is the pressure. This flux is continuous across the interface so the pressure gradient either side of the interface is proportional to the viscosity of the fluid there. If the pressure gradient is larger in the ambient fluid then perturbations to the interface into the ambient fluid experience a larger pressure gradient and tend to grow. Hence viscous fingers can arise if the ambient fluid is of greater viscosity. In a three-dimensional porous layer with fluids of different densities that have been vertically segregated by buoyancy, the input fluid will form an intrusion with large radial extent at later times. The radial flux at the interface no longer needs to be continuous because the input intrusion advances into the ambient fluid. Indeed, we will show that the vertical structure of the flow ensures that there is no destabilising pressure gradient across the interface. This stabilisation occurs for any density difference between the two fluids that eventually leads to buoyancy segregation.
The combination of buoyancy segregation and an intrusion of large radial extent is a fundamentally different stability mechanism to the related problem of two-layer viscous gravity currents where two viscous fluids co-flow driven by gravity alone. In this case, buoyancy segregation does not guarantee stability. Instead, stability is controlled by competing contributions to the pressure gradient from buoyancy and the viscosity contrast (Kowal & Worster 2019). The Saffman-Taylor instability is suppressed provided that the density difference is sufficiently large. A similar phenomena occurs in unconfined flows with a longitudinal viscosity contrast (Kowal 2021).
The stability of single-layer unconfined gravity currents (in which the aquifer is bounded only on one side by a single impermeable layer), has been investigated extensively (Pattle 1959;Grundy & McLaughlin 1982;Newman 1984;Bernoff & Witelski 2002;Chertock 2002;Mathunjwa & Hogg 2006). In this case, the motion of the ambient fluid is unimportant and the gravity current is not resisted by its displacement; the solutions are always stable. The stability of confined gravity currents is different because the displacement of the more viscous ambient fluid influences the physics. For these confined flows, stability has been established for perturbations that do not depend on the transverse or azimuthal coordinate (Nordbotten & Celia 2006). The stability analysis of both unconfined gravity currents and axisymmetric variations to confined gravity currents is significantly simplified because the driving pressure gradient in the input fluid can be eliminated from the problem. The present study of three-dimensional perturbations to confined currents requires the solution to a coupled system for both the pressure gradient and the interface shape. The paper is structured as follows. The model for the buoyancy-segregated displacement flow is formulated in §2. The coupled system that governs the background pressure and the interface shape is integrated numerically in §3 with the results demonstrating that the axisymmetric solutions are stable to angular-dependent perturbations for all values of the viscosity ratio and the density difference. In §4, the physical mechanism that suppresses the Saffman-Taylor instability is analysed. A linear stability analysis is carried out in §5 for the case of buoyancy-segregated fluids with an arbitrarily small density difference. In §6, we show that the stabilising effect of buoyancy segregation generalises to layers with vertically varying permeability for which the interface may contain steep shock-like regions. Concluding remarks are made in §7.

Model
The setup is shown in figure 1. A similar geometry was studied by Nordbotten & Celia (2006) and Guo et al. (2016). Fluid of density and viscosity is injected into a porous medium of thickness 0 and porosity . The medium is initially filled with a second fluid of greater density + Δ and different viscosity . The top and bottom boundaries of the layer are impermeable. The input fluid is injected with volume flux from a point source at the upper boundary, which we take to be the origin of our coordinate system. The coordinate is measured downwards from the upper boundary at which = 0. We use cylindrical polar coordinates with denoting the distance from the axis and denoting the angular coordinate. Time is denoted by . We neglect any intermingling of the fluids and assume that there is a sharp interface at = ( , , ) between the fluids. We also assume that the two fluids are segregated by buoyancy so that is a single-valued function.
When the input fluid has spread far from the source, the flow is predominantly in the radial direction. We may then apply the shallow flow approximation and the pressure is approximately hydrostatic, in the input and ambient fluids, respectively, where¯ ( , , ) is the pressure at the upper boundary, = 0, which is to be determined. The Darcy velocities are given by (2.4) in the input and ambient fluids, respectively. The governing equations of the flow are + ∇ · U = 0 (2.5) and ∇ · U + ( 0 − )U = 0, (2.6) which represent mass conservation integrated over the thickness of the input fluid, and mass conservation integrated over the thickness of the layer, respectively. At = 0, the porous layer is entirely filled with the ambient fluid and constant-flux injection of the input fluid begins. Global mass conservation of the input fluid takes the form where ( , ) denotes the frontal contact line where = 0 (see figure 1b). Eventually, the fluid-fluid interface will touch the lower boundary (Nordbotten & Celia 2006). It is well-known that the predicted interface shape of an unconfined axisymmetric gravity current has an unphysical singularity at the origin, whilst in confined layers, the flow will always encompass the entire layer near the source, regardless of the thickness of the layer (Lyle et al. 2005;Guo et al. 2016). Recently, Benham et al. (2022) carried out a detailed analysis of the flow near the source, incorporating the pressure-driven vertical flow there. Their results corroborate the conclusion that a confined layer becomes filled with input fluid at sufficiently late times (see also Huppert & Pegler 2022). After the fluid-fluid interface has touched the lower boundary, there is a second contact line at = ( , ) where = 0 . The flow then has three regions: there is the 'fully-flooded' zone in which the input fluid fills the layer ( = 0 ; 0 < ( , )); next is the partially filled zone in which the fluid-fluid interface lies between the top and bottom boundaries (0 < < 0 ; ( , ) < < ( , )); finally, beyond this there is an uninvaded zone in which the ambient fluid fills the layer ( = 0; > ( , )); see figure 1b. When the input fluid has filled the layer near the source, the boundary conditions are as follows. At the source, where the latter arises from volume conservation. In the far field, the boundary condition is (2.9) The dependent variables,¯ ( , , ) and ( , , ), are defined on ∈ [0, 2 ) with periodic boundary conditions:¯ ( , 0, ) =¯ ( , 2 , ) and ( , 0, ) = ( , 2 , ). We also require that 0 0 . Equations (2.5) and (2.6) combined with global mass conservation (2.7) and the boundary conditions (2.8), (2.9) define a complete coupled system for determining the two dependent variables¯ ( , , ) and ( , , ).
Although the system derived above is sufficient to solve the problem numerically (see §3), the stability investigation requires detailed treatment of the behaviour in the vicinity of the two contact lines where the interface touches the layer boundaries. Across these lines, the gradient of the interface and the gradient of the pressure at the top of the layer can be discontinuous. Indeed, it is well-known that the stability of the fluid-fluid interface for equally dense fluids is controlled by the jump in the pressure gradient at the interface (see §1). Therefore, we derive boundary conditions for the interface shape and the pressure gradient across each contact line.
At the trailing contact line, = ( , ), continuity of the interface, continuity of the pressure and continuity of the flux of the input fluid take the following forms respectively, where [·] + − denotes the discontinuity in the quantity in square brackets at the contact line = ( , ). In addition, we have a kinematic boundary condition; at the trailing contact line, the radial velocity of the ambient fluid is equal to the radial velocity of the interface, which takes the following form (2.11) At the leading contact line, = ( , ), the corresponding boundary conditions are (2.14) where [·] + − denotes the discontinuity in the quantity in square brackets at the contact line = ( , ).

Transformed coordinate system
Mass conservation of the input fluid suggests that the radial extent of the currents grows with 2 ∼ /( 0 ). This motivates introducing the following transformed dimensionless coordinates where = ref > 0 is a reference time, which corresponds to = 0. We also scale the dependent variables to obtain the following dimensionless quantities The two mass conservation equations (2.5, 2.6) become is the gradient operator with respect to the transformed coordinates and we have introduced the following two dimensionless groups which are the viscosity ratio and gravity number, respectively. The latter represents the magnitude of pressure gradients associated with buoyancy relative to pressure gradients associated with injection of fluid. In the case of equally dense fluids, the gravity number is = 0. The boundary conditions at the source and in the far field are and The trailing contact line in the transformed coordinates is denoted by = ( , ) and the boundary conditions there (2.10, 2.11) are recast as whilst at the leading contact line, = ( , ), Mass conservation (2.7) takes the form 1 2 (2.28)

Self-similar axisymmetric solution
The so that the contact lines are circles, fixed in similarity space. We are interested in analysing the stability of and convergence to these solutions in § §3, 4 and 5. In the present section, we derive a single ordinary differential equation which governs the axisymmetric self-similar interface shape and we recall an important analytical solution to this equation from Nordbotten & Celia (2006) and Guo et al. (2016). For axisymmetric flow, the pressure at the top boundary, P 0 ( ), can be eliminated from the problem as follows. We integrate equation (2.18) with respect to and apply the source flux boundary condition (2.21) to obtain Upon substituting (2.31) into (2.17), we obtain the following ordinary differential equation for (2.32) This axisymmetric system can be solved numerically for H 0 ( ) as described in Guo et al. (2016) for a wide range of and . For ≫ 1 (buoyancy-dominated flow), the input fluid forms a thin current near the top of the layer and the displacement of the ambient fluid is unimportant. The flow is effectively 'unconfined' and single phase. Hence, we do not expect Saffman-Taylor instabilities in this regime. Instead, we focus on the case of relatively larger input flux ( O(1)) and a less viscous input fluid ( < 1) for which stability is not well-understood. For ≪ 1 and < 1, (2.32) has an approximate analytic solution. We neglect the right-hand side of (2.32). Although this removes the highest derivative of H 0 , the boundary conditions (2.27) and (2.24) become independent of H 0 for ≪ 1 and simply give the locations of the contact points so the problem remains complete. Integrating the left-hand side of (2.32) furnishes the following interface shape (Guo et al. 2016) where the contact lines are given by Note that these provide inner bounds on the locations of the contact lines for general > 0 provided that < 1: 37) with equality as → 0 + . These inequalities arise because buoyancy acts to extend the interface (its effect is proportional to ).
The corresponding pressure at the top of the layer takes the form where P 0 ( ) is defined up to addition of an arbitrary constant.
The numerical solution to (2.32) for H 0 is plotted as a blue line in figure 2a for = 0.1 and = 0.2. The small-analytical solution (2.34) is shown as a black dashed line and there is good agreement with the numerical result. Figure 2b shows the analytical solutions for various values of < 1. When the input fluid is of relatively lower viscosity (smaller ), it forms a finger that intrudes further into the ambient fluid at the top of the layer.
The axisymmetric analytic solution (2.34) is relevant to the regime in which pressure gradients associated with buoyancy are weak compared to those associated with injection ( ≪ 1). It is worth considering the behaviour in the limit → 0, corresponding to equally dense fluids and no buoyancy force. As → 0, the axisymmetric analytic solution (2.34) becomes a more accurate approximation of the numerical solution, with the input fluid preferentially flowing near the upper boundary (figure 2). However, = 0 corresponds to no buoyancy segregation and so the input fluid should have no preference for the top of the layer. Indeed, in the case of equally dense fluids or an arbitrarily thin layer ( = 0), and < 1, one would observe classical Saffman-Taylor fingering with angular dependence and no vertical preference (Paterson 1981). The derivation of the axisymmetric analytical solution (2.34) implicitly assumes that the input and ambient fluid have been segregated by buoyancy, which requires > 0 and takes a time proportional to −1 . Hence, the case of = 0 has qualitatively different late-time behaviour to the case of small but non-zero and the limit → 0 is singular. In this paper, we analyse the stability when is small but positive; we exclude the case = 0 for which buoyancy has no effect.
For small but non-zero, the interfaces in figure 2b occur provided that (i) buoyancy segregation has occurred and (ii) the radial extent of the flow is much greater than the layer thickness so that the shallow approximation applies. For < 1, in dimensional terms, the former requires Prior to buoyancy segregation, Saffman-Taylor fingering can occur. The present paper is concerned with the post-segregation stability of the axisymmetric solutions to (2.32) to angular-dependent perturbations with Δ > 0. For ≪ 1, the first condition (2.41) implies the second (2.42). Notice that for small density differences, Δ , the time for buoyancy segregation is very large. It is also important to note that in an anistropic porous medium, the permeability that appears in segregation timescale (2.41) will depend only on the vertical permeability, whereas the gravity number (and hence the self-similar solutions) depends only on horizontal permeability (cf. Benham et al. 2022). The results we derive in this paper concerning flow stability extend to porous layers that have different vertical and horizontal permeabilities (as is common in many aquifers); see also §6. The case of more complicated anisotropy such as different horizontal permeabilities in different directions is beyond the scope of the present work.
For simplicity, we have modelled the case where the source of input fluid is located at a point on the upper boundary. However, the analysis and results of this paper apply for any location of the injection source (for example, at the lower boundary, or over a vertical line within the layer). This is because the input fluid will always eventually fill the thickness of the layer near the source. Subsequently, the flow becomes predominantly radial. There will, of course, be different early-time transient behaviour for different source locations but once (2.41) and (2.42) apply, the exact source location becomes unimportant.

Numerical results
In order to investigate the response of the flow to perturbations with angular dependence, we develop a numerical method to solve the coupled system for the evolution of the fluid-fluid interface, H ( , , ), and the pressure at the upper boundary, P ( , , ). This consists of equations (2.17) and (2.18) with boundary conditions (2.21) and (2.22), and an appropriate initial condition.  We solve this system numerically using a finite-difference scheme, details of which are given in Appendix A.
The numerical method is applied to a wide range of axisymmetric and non-axisymmetric initial conditions over many values of the parameters, and . In all cases, the solution converges to the axisymmetric similarity solution. This demonstrates both the veracity of the numerical method and that the similarity solution provides the intermediate asymptotics for many initial conditions (cf. Mathunjwa & Hogg 2006). Examples are shown in figures 3 and 4 illustrating how mode six and mode three fingers are suppressed even though the input fluid is of lesser viscosity than the ambient. We find that any perturbation is suppressed, including radial perturbations, provided that the fluids remain buoyancy segregated (if the interface becomes nonmonotonic as a function of vertical coordinate then the model breaks down).
The numerical results also demonstrate that there is a jump in the radial pressure gradient at the upper boundary across the leading contact line (see figure 5). This pressure gradient jump is associated with the density difference between the two fluids and the discontinuity vanishes as → 0. The discontinuity has a stabilising effect because the magnitude of the pressure gradient is smaller in the ambient fluid (see §4). However, this discontinuity in the hydrostatic pressure is not required for the interface to be stable. Indeed, in §5 we show that the interface is stable even in the limit → 0.

Stability mechanism
In the present section, we describe the mechanism by which buoyancy segregation suppresses viscous fingering. In the classical case of radial flow with no vertical dimension (discussed in §1), the pressure gradient either side of the fluid-fluid interface is proportional to the fluid viscosity owing to continuity of the radial flux, which drives the well-known instability.
We now analyse the more complicated three-dimensional case in which the fluids have different densities and are segregated by buoyancy. Between the contact lines, the pressure gradient in the input fluid is dP 0 /d and the pressure gradient in the ambient fluid is d(P 0 − H 0 )/d , where the − H 0 contribution arises from the density difference between the fluids. Since both dP 0 /d and dH 0 /d are negative, the magnitude of the pressure gradient is larger in the input fluid (in contrast to the classical case with no vertical dimension), which implies that between the contact lines the interface is stable for > 0. The stability at the contact lines must be treated separately because there is a discontinuity in dP 0 /d and dH 0 /d between the input and ambient fluids as shown in figure 5. We consider the driving pressure gradients in the input and ambient fluids either side of the contact lines (these locations are labelled (i)-(iv) in figure 2a).
At location (i), just ahead of the leading contact line, at = + 0 , the interface gradient and pressure gradient at the top boundary are The latter is calculated from (2.31). The pressure gradient driving the ambient fluid at (i) is At location (ii), just behind the leading contact line, = − 0 , the interface and pressure gradients The former is calculated by taking H → 0 in (2.32). There is a discontinuity in the pressure gradient at the top boundary, dP 0 /d across = 0 , which can be observed in the numerical results in figure 5. This discontinuity vanishes as → 0 even for fluids with different viscosities; the jump in the pressure gradient associated with the viscosity contrast is eliminated by the radial extension of the interface. The pressure gradient driving the input fluid at (ii) is given by (4.3)b, and its magnitude is greater than (or equal to) the magnitude of the pressure gradient in the ambient fluid (4.2) because 0 2/ (see 2.37). Hence small perturbations to the contact line decay because they experience unfavourable pressure gradients, which suggests interfacial stability provided that > 0. The driving pressure gradients either side of the contact line are equal in the limit → 0; this case is discussed in more detail in §5.
Similar analysis applies at the trailing contact line. At location (iii), just downstream of the trailing contact line, = + 0 , the interface gradient, the pressure gradient at the upper boundary, and the pressure gradient driving the ambient fluid are given by respectively. At location (iv), just behind the trailing contact line, = − 0 , the interface gradient and the pressure gradient at the upper boundary (which drives the input fluid) are given by The magnitude of the pressure gradient driving the input fluid, 1/ 0 , is greater than (or equal to) the magnitude of the pressure gradient driving the ambient fluid, 0 /(2 ) since 0 √ 2 with equality when = 0. This suggests that the trailing contact line is stable to small perturbations for > 0.
It is worth noting that for small , at the location where H 0 = 0, the pressure gradient in the ambient fluid (4.1)b is proportional to the relative viscosity of the ambient fluid divided by the radial distance and similarly at H 0 = 1, the pressure gradient in the input fluid (4.5)b is proportional to the relative viscosity of the input fluid divided by the radial distance. This dependence on the relative viscosity and radial location is as in the classical case (see 1.1). However, the key difference is that, in contrast to the classical instability, the contact lines are separated so that the radial coordinate is different at H 0 = 0 and H 0 = 1. The change in the pressure gradient between location (iv) and location (i) (figure 2a) is controlled by the viscosity contrast (relatively larger ambient viscosity leads to a larger increase in the magnitude of the pressure gradient) and how far apart the contact lines are, which reduces the magnitude of the pressure gradient downstream. These two effects act in opposition. If the intrusion of input fluid has a relatively short radial extent then the flow is unstable and the intrusion grows radially until the distance between the contact lines acts to stabilise the interface. Hence the self-similar axisymmetric solutions have larger radial extent at smaller viscosity ratios. In summary, the radially extensive intrusion of the interface along the upper boundary associated with buoyancy segregation counteracts the destabilising discontinuity in the pressure gradient associated with the viscosity contrast. The vertically-segregated radial intrusion could be thought of as a single axisymmetric viscous finger, which is stable to angular-dependent fingers.
We have shown that it is not gravity that acts against the pressure gradient to stabilise the contact lines. Indeed, the radial extent of the interface is insensitive to the size of for small . The role of gravity is simply to segregate the fluids and stability occurs after the segregation has occurred. In section 5, we show that provided that buoyancy has segregated the fluids, the small-axisymmetric self-similar solutions are linearly stable. This formally confirms the results in the present section.

Linear stability for small density difference ( ≪ 1)
In the present section, we demonstrate the linear stability of the ≪ 1 axisymmetric selfsimilar solutions (given by 2.34). Although buoyancy does not appear in the form of the solution, some buoyancy is required to segregate the fluids. Confirming the stability of the ≪ 1 solution will demonstrate stability for > 0 as discussed in §4. The linear stability analysis also gives the rate of decay of each mode and its dependence on the viscosity ratio.
We consider -dependent perturbations to the axisymmetric self-similar solutions of the form where ≪ 1 and is an integer. We seek to determine the stability of these perturbations. Note that the perturbation corresponds to ∼ in terms of real time (see 2.15). Often linear stability analyses of viscous gravity currents requires rescaling the spatial domain owing to the singular behaviour of the interface at the contact lines (e.g. Mathunjwa & Hogg 2006;Kowal & Worster 2019). For the present porous gravity current, the interface is linear at the contact lines so such a transformation is not required. Instead, we linearise the boundary conditions about the leading order locations of the contact lines, = 0 and = 0 . We consider perturbations with -dependence ( 1) first as the stability in this case has not yet been investigated. Note that global mass conservation of the input fluid (2.28) is identically satisfied by the form of the perturbations for 1. The case of axisymmetric perturbations ( = 0) is included at the end of this section for completeness.
In the single-phase regions (upstream of the trailing contact line and downstream of the leading contact line), the pressure satisfies Laplace's equation and the hence pressure perturbation is given by the solution to Given that the pressure perturbation remains finite as → ∞ and as → 0, the solution in the single-phase regions is P 1 = for < 0 (5.6) where and are constants. After some algebra, the linearised governing equations (2.17), (2.18) in the interface region ( 0 < < 0 ) become We can then eliminate H 1 and obtain a single equation for P 1 in the interface region between the -1 -0.9 -0.8 -0.7 -0.6 -0.5 -1 -0.5 0 0.5 1 F 6. The dispersion function, F ( , , ) (5.21) plotted versus the growth rate for = 0.5. It has been normalised by its maximum value over the range of , F max . The zeros of the curves correspond to admissible decay rates (5.20) with the exception of the largest zero (shown as stars). The next largest zero (shown as circles) will be the decay rate observed (further details in the text).
contact lines: To determine the pressure perturbation in the interface region, we solve (5.10) with three boundary conditions at each contact line, which arise from the linearised version of continuity of the pressure, continuity of the flux, and the kinematic boundary condition (see 2.23, 2.24, 2.25, 2.27). At = + 0 , these boundary conditions take the form where we have used the upstream behaviour (5.6). At = − 0 , the analogous boundary conditions take the form where we have used the downstream behaviour (5.7). The system for P 1 between the contact lines is linear and has six boundary conditions with six unknown constants: , , 1 , 1 and two constants arising from solving the ordinary differential equation (5.10). The system is an eigenvalue problem for the growth rate . The general solution to equation (5.10) is characterised by the value of relative to the critical value = −2 2 − 1 4 2 + 1 . This critical value of corresponds to a repeated root in the auxiliary equation for the ordinary differential equation (5.10). There are no non-trivial solutions that satisfy all the boundary conditions for . For −1 < < , the solution takes the form and 1 and 2 are constants. We apply the six boundary conditions to obtain the following dispersion relation for The dispersion equation (5.20) can be solved to obtain the growth rate, for any 1 and 0 < < 1. The corresponding eigenfunctions are given by (5.18). Since < < 0, the growth rate is always negative and hence the axisymmetric similarity solutions are linearly stable.
The dispersion equation (5.20) can have multiple solutions. The function F ( , , ) is plotted against in figure 6 for = 0.5 and four values of . For each , the largest zero for is = (these solutions are indicated by stars in figure 6). However, the case = corresponds to a trivial solution for P 1 , which is not admissible. The next largest solution of (5.20) corresponds to the slowest decaying non-trivial solution and these values are indicated by circles in figure 6. We expect this to be the decay rate that is observed and we take this solution as the correct value of .
The eigenfunctions (5.18) that correspond to this correct solution are defined up to a multiplicative constant. The form of the interfacial perturbation, H 1 for each eigenfunction (5.18) can be calculated from (5.9). These interfacial perturbation eigenfunctions are shown in figure 7 for different values of the viscosity ratio, and the mode, . The shapes are more oscillatory at higher values of .
The growth rate, , calculated for different values of and using the dispersion relation (5.20) is plotted using continuous lines in figure 8. The rate of decay is slower for larger values of and smaller viscosity ratios, . The crosses in figure 8 denote the decay rate calculated from the numerical method for = 0.5 and = 0.8 for three values of with = 0.01; details of this calculation are given in Appendix A.1. There is excellent agreement between the numerically-derived prediction for and the values obtained from the linear stability analysis.
Finally, we consider the case = 0. In < 0 , the pressure perturbation is constant, P 1 = 0 , whilst in > 0 , P 1 = 0 . In the interface region, the pressure perturbation is linear; P 1 = + . The kinematic boundary conditions (5.13, 5.16) furnish = −1. The thickness perturbation takes the form (see 5.9) The numerical results for axisymmetric perturbations to the self-similar solutions corroborated that = −1 (shown as a cross at = 0 in figure 6). Axisymmetric perturbations decay faster than -dependent perturbations with the error decaying proportional to −1 as has been found previously for unconfined gravity currents (

Vertically varying permeability
The results obtained in previous sections generalise to layers in which the permeability varies vertically with = ( ). The flow is axisymmetric and self-similar after buoyancy segregation. Here, we find the analytical solutions that arise in the small regime with a less viscous input fluid ( < 1). One key difference to the uniform case is that shock-like regions in which the interface is relatively steep can occur even for < 1 (see §6.1). Despite this, we show that at late times the axisymmetric solutions for any continuous permeability profile, ( ), are stable. We define the dimensionless permeability as where the denominator is the mean permeability. We also define the dimensionless depthintegrated permeability (Hinton & Woods 2018), This model applies provided that the fluids have become segregated by buoyancy. In a uniform layer this applies at times given by (2.41) with representing the vertical permeability. In a layer with vertically varying permeability and less viscous input fluid ( < 1), the dimensional time for buoyancy segregation is adjusted to where ( ) is the vertical permeability. In an isotropic, vertically heterogeneous layer, ( ) = ( ). The time (6.5) becomes singular if the permeability is zero at any height in the layer (and in this case our model does not apply). The time for the shallow approximation to apply is given by (2.42).
As for a uniform layer, there are self-similar axisymmetric solutions, H 0 = H 0 ( ), P = P 0 ( ) to (6.3), (6.4), with the pressure gradient at the top boundary given by and H 0 satisfies the following ordinary differential equation The interface gradient and pressure gradient at the leading contact line are (6.8) respectively. Whilst at the trailing contact line The small-analytic solutions are given by neglecting the right-hand side of (6.7) from which we obtain the following implicit solution , 0 <H 0 < 1, (6.11) (6.12) where the contact lines are The associated pressure gradient in the interface region is given by . (6.14) In a uniform layer, the interface (6.11) is monotonic whenever < 1 but in a layer with a vertical permeability variation, turning points can arise even when the input fluid is less viscous than the ambient (Hinton & Woods 2018). The solution (6.11) is only valid when the interface is monotonic. Otherwise, the interface has a turning point with buoyant input fluid lying below denser ambient fluid. This would invalidate the model, which assumed that the fluids have segregated owing to buoyancy. For example, for a linear permeability structure with dimensionless variation, Δ , (6.16) Provided (6.11) is monotonic, the interface forms an axisymmetric self-similar intrusion, which extends along the upper boundary in a qualitatively identical fashion to the case of a uniform layer (see figure 9a). In the case that (6.11) has a turning point, a shock must be introduced; this regime is studied in §6.1. We apply our numerical method to layers with vertical variations in permeability for a wide range of initial conditions and find that the self-similar axisymmetric solutions are stable to both -dependent and axisymmetric perturbations for any > 0. Next, we generalise §4 to show how the pressure gradients suppress instabilities in a layer with vertically varying permeability (but no shock-like regions; for that case, see §6.1). Between the contact lines, the magnitude of the pressure gradient in the ambient fluid is less than or equal to that in the input fluid owing to the contribution of the term − H 0 / (as in §4). The interface gradient at the contact lines is discontinuous and the stability there is treated separately.
We can bound the location of the contact lines of the self-similar axisymmetric solution. In the case that the small-analytic solution (6.11) is monotonic, the contacts lines for > 0 satisfy 0 2K (0) , 0 2 K (1), (6.17) with equality as → 0 + . These inequalities reflect the action of buoyancy to extend the interface. The gradient of the interface and the pressure at the upper boundary at the contact lines are given by (6.8) and (6.9). Thus, the pressure gradient in the ambient fluid just ahead of the leading contact line is given by whilst in the input fluid just upstream of the leading contact line, the pressure gradient is 19) and the inequality (6.17)a ensures that the magnitude of the pressure gradient is larger in the input fluid, which stabilises the interface, with equality as → 0 + . The pressure gradient in the ambient fluid just ahead of the trailing contact line is given by whilst in the input fluid just upstream of the trailing contact line, the pressure gradient is (6.21) and the inequality (6.17)b ensures that the magnitude of the pressure gradient is larger in the input fluid, which stabilises the interface.

Shock-like regions
In the case in which the small-interface (6.11) has a turning point, a shock must be introduced to vertically segregate the fluids. This shock may occupy part or all of the layer (see figure 9b). By way of an example, we consider a linear variation in permeability with Δ = 1.5 and viscosity ratio, = 0.3. The axisymmetric self-similar profile has a shock-like region near the top of the layer, even though < 1 (red line in figure 9b). The shock at the top of the layer is at location = and is of thickness H , which satisfy These arise from continuity with the solution (6.11), which is valid in H > H , and mass conservation, respectively. The shock-like region may also extend across the entire layer (e.g. black line in figure 9b). For more complicated permeability profiles, there can be multiple shocklike regions separated by rarefaction-like regions (Hinton & Woods 2018). We apply our numerical method to layers with shock-like regions and find that the selfsimilar axisymmetric solutions are stable to both -dependent and axisymmetric perturbations. An example is shown in figure 10 with a linear variation in permeability with Δ = 1.5, viscosity ratio = 0.3, and = 0.05, demonstrating that mode five fingers are suppressed. The smallanalytic solution is shown as a dashed line in in figure 10c; there is a shock-like region in the top part of the layer. For > 0, the shock-like regions are smoothed by buoyancy so that they are not vertical (Hinton & Woods 2018). The smoothed shock-like region has radial extent proportional to in the similarity coordinate . Hence, in dimensional terms, the aspect ratio of the shock and eventually the shock-like region of the flow is long and thin (when the right-hand side becomes large) and the shallow flow approximation applies at sufficiently late times. This smoothing of the shock-like region also occurs when the shock encompasses the entire thickness of the layer (e.g. figure 9b). Thus the fluid-fluid interface is not cylindrical in the case. Instead, due to buoyancy, the interface becomes gradually shallower over time, which suppresses the classical instability. The stability analysis generalises to the case in which there are shock-like regions. At the top of the layer, the leading contact line for > 0 satisfies 0 2K (0) . (6.24) The first inequality arises from buoyancy smoothing the interface beyond the shock location (see figure 10c) whilst the second arises because the profile (6.11) has a turning point in the shock-like region and the shock conserves mass so it must be ahead of this profile at = 0. Hence (6.17)a applies even when there is a shock-like region at the top of the layer. A similar argument applies when there is a shock-like region at the bottom of the layer and (6.17)b remains correct. The gradient of the interface and the pressure at the upper boundary at the contact lines are given by (6.8) and (6.9), even when there are shock-like regions. Hence the stability argument of §6 applies. In summary, the self-similar axisymmetric solutions in a layer with vertically varying permeability are stable owing to the same mechanism as in a uniform layer. The less viscous input fluid intrudes into the ambient fluid sufficiently far that the destabilising increase in the magnitude of the pressure gradient associated with the viscosity contrast is smoothed over the extent of the interface. Buoyancy segregation corresponds to a monotonic interface and the input fluid may be forced into the low permeability region (see also Debbabi et al. 2018). We note that there may be significant intermingling of the fluids in the transient evolution to the buoyancy-segregated flow (e.g. Huppert et al. 2013). There may also be a competition between viscous fingering and heterogeneity-driven fingering prior to buoyancy segregation (De Wit & Homsy 1997).

Discussion and conclusion
We have examined the stability of the horizontal flow of one fluid injected into another fluid of greater viscosity and density. When there is a sharp interface and the two fluids have segregated owing to buoyancy, the flow evolves in an axisymmetric self-similar fashion. Whilst Saffman-Taylor fingers may occur prior to buoyancy segregation, we have demonstrated that the buoyancy-segregated self-similar flow is stable to both axisymmetric and angular-dependent perturbations. The input fluid forms an intrusion with large radial extent neighbouring the upper boundary. This disperses the destabilising pressure gradient jump associated with the viscosity contrast between the fluids that drives the classical instability. The mechanism generalises to layers with a vertical variation in permeability and anisotropic layers with different horizontal and vertical permeability. After buoyancy segregation, the flow in a heterogeneous layer is always stable even though the interface may contain steep shock-like regions owing to slow flow in the low permeability zone. The time for buoyancy segregation can be long if the porous layers has zones of very low permeability.
In future work, these ideas could be extended to viscous displacements in horizontal channels as is relevant to mantle plumes (Schoonman et al. 2017). The base self-similar flow was found for a line source and a point source by Zheng et al. (2015b) and Hinton (2020) but the stability has not yet been investigated. The analysis would be more complicated than the present work owing to the effect of the no-slip boundaries and the shear flow (Snyder & Tait 1998;John et al. 2013).
In the context of CO 2 sequestration, viscous fingering is undesirable as it may reduce storage efficiency. Our study could be usefully extended to consider strategies to vary the injection rate to avoid viscous fingering. Given that the interface is stable after buoyancy segregation (even for relatively high injection rates) one could analyse how to vary the input flux during the presegregation transient to ensure that fingers never occur.