## Abstract

We consider the displacement of one fluid by a second immiscible fluid through a long, thin permeable channel whose thickness and permeability decrease away from the axis of the channel. We build a model that illustrates how the shape of the fluid–fluid interface evolves in time. We find that if the injected fluid is of the same viscosity as the original fluid, then the cross-channel variations in permeability and thickness tend to focus the flow along the centre of the channel. If the viscosity of the injected fluid is smaller than the original fluid, then this flow focusing intensifies, leading to very poor sweep of the original fluid in the system, with the injected fluid bypassing much of the channel. We also show that if the viscosity ratio of the injected fluid to the original fluid is sufficiently large, then a blunt nose may develop at the leading edge of the injected fluid, whereas the remainder of the fluid–fluid interface becomes stretched out along the edges of the channel. This leads to a much more efficient sweep of the original fluid from the channel. We generalize the model to illustrate how buoyancy forces and capillary pressure affect the evolution of the system and compare our model predictions with some simple laboratory experiments. This partial stabilization of a fluid interface in a channel of non-uniform width represents a generalization of the classical Saffman–Taylor instability, and our nonlinear solutions for the evolution of the interface highlight the importance of cross-channel variations in permeability and thickness in modelling flow in channelled reservoirs.

This article is part of the themed issue ‘Energy and the subsurface’.

## 1. Introduction

There has been considerable interest in the dynamics of two-phase flow in porous media with important applications for oil recovery and also carbon sequestration [1,2]. During the enhanced recovery of oil, water is often injected into a porous reservoir rock to displace the oil [3]. In such flows, the injected fluid may be less viscous than the original fluid in the formation, and this can lead to a viscous instability in which the injected fluid develops fingers that migrate through the original fluid [4,5].

Saffman & Taylor [4] modelled the viscous instability of a fluid–fluid interface migrating through a long, thin Hele-Shaw cell of finite and constant thickness. Their analysis has been developed by many subsequent researchers for both Hele-Shaw geometries and porous layers [5]. In the case of a gap of constant width, there have been detailed numerical, analytical and experimental works exploring the nonlinear development of the instability.

A problem of particular importance relates to the mathematical description of the nonlinear fingers, which develop as the instability grows. In the original work of Saffman and Taylor, a solution was presented for the nonlinear shape of the viscous finger, using potential theory and assuming a width for the finger as a fraction of that of the channel. Subsequent work has identified a selection mechanism for the finger width in terms of the flow far upstream [6,7].

In situations such as enhanced oil recovery, prediction of the nonlinear evolution of the instability is key for quantifying the areal sweep of the reservoir, and hence the fraction of the oil which is recovered through the injection of water. In some situations, such as enhanced oil recovery, viscosifiers are added to the fluid to suppress the instability, so that the less viscous fluid is displaced by a more viscous one [5]. Several analyses have been carried out to minimize the growth rate of the instability, often assuming a simplified geometry for the system [8,9]. More recently, Al-Housseiny *et al.* [10] demonstrated that if the thickness of the Hele-Shaw cell decreases in the direction of the flow, then the gradual increase in the surface tension across an advancing immiscible front, which results from the gradual decrease in the radius of curvature, can stabilize the interface.

In the related problem area of carbon sequestration, CO_{2} may be injected into deep saline aquifers, displacing the original water in place [11]. During the displacement of water by CO_{2} injected into a porous layer, it is important to understand the evolution of the fluid–fluid interface as this informs estimates of the fraction of reservoir pore space which is invaded by CO_{2}. As the flow spreads from the injection well, the flow is initially controlled by the source pressure, and again may be hostage to the viscous fingering instability mentioned earlier. In porous layers of large vertical extent, the flow typically transitions to a buoyancy-dominated flow, and the motion is often described using models of porous gravity currents [12,13]. As such a flow spreads through the porous layer, the CO_{2} may disperse around heterogeneities in the formation [14], and, ultimately, it may dissolve into the ambient water in which it is weakly soluble [15–18]. Over long time scales, in excess of 10 000 years, any background hydrological flows may enhance this dissolution through supply of CO_{2}-unsaturated aquifer water [19]. Much of the analysis of these flows has been developed for homogeneous isotropic porous media, and tested experimentally using a Hele-Shaw cell, which consists of a narrow constant gap between two parallel plates [2].

However, for both of the above-mentioned processes of enhanced oil recovery and also carbon sequestration, the flow of the fluid–fluid interfaces may be strongly influenced by the structure and boundaries of the flow domain. Geological deposits in many sedimentary systems include relatively long and narrow channelled deposits associated with specific geological events; these channels are often confined by layers of low-permeability shale and they may exhibit variations in thickness, grain size and permeability across the channel, as the deposit typically thins and becomes more fine-grained on the flanks or wings of the channel. Such channelled deposits may become stacked vertically, owing to multiple phases of deposition over geological time, and this can lead to preferential pathways for fluid to migrate through the formation ([20]; figure 1). Seismic studies of a number of turbidite deposits have illustrated that, in the channelled region of the flow, the thickness of the deposit wanes towards the outer edge of channels [21–23] and geological mapping shows that, as well as this decrease in thickness, there can be a trend of increasing interbedding of sandstones and mud-stones, as at the Ross Formation, County Clare, Ireland [24], and the vertically stacked channel complexes in the Karoo, South Africa [25]. There may also be a gradual decrease in the grain size towards the flanks perhaps associated with stratification in the flow [26] or hydrodynamic fractionation [27], as reported at the turbidite deposit at San Clemente State Beach, CA, USA [28]. These structural features can lead to a reduction in the permeability and depth of the channel from the axis to the wings of the deposit.

When water is injected into such a formation to displace oil, if the fluid invades one of the channels, we expect that it will first migrate along the central section of the channel, in the region of highest permeability and thickness. Indeed, there are numerous examples of flow through heterogeneous permeable media where the flow tends to follow the most permeable pathway [2,29,30] and, as a result, the injection of a viscous fluid then leads to an effective stabilization of the flow as the more heterogeneous part of the formation tends to be flooded first. In the present situation, of a channel with non-uniform permeability in the cross-channel direction, a wave of injected fluid of increasing width will develop, gradually displacing more of the original fluid from the flanks of the channel. Similarly, if CO_{2} is injected as part of a sequestration scheme, the variation of properties from the centre to the flank of a channel will tend to focus the highly mobile CO_{2} along the axis of the channel if this is the most permeable pathway. Such localization of the flood front can limit the effective storage potential of a reservoir as the flow bypasses the wings of the channel and therefore only accesses a fraction of the pore space.

The objective of this paper is to develop a model of such topographically controlled viscous fingering in a channelled formation and assess its impact for the injection of both relatively viscous and relatively inviscid fluids, with application to enhanced oil recovery and CO_{2} sequestration. To this end, we build a model to describe the displacement flow that develops in a porous channel with permeability and thickness that vary in the cross-channel direction. In this situation, we expect the flow speed to be greater in the centre of the channel. This leads to an evolution of the shape of the interface along the channel (figure 2) and informs a calculation of the fraction of the cross-sectional area of a channel that is swept by the displacing fluid as a function of time. On initial consideration, one might infer that, even if the displacing fluid is more viscous, the interface between the two fluids will become stretched out, as the injected fluid bypasses much of the original fluid and runs along the most permeable region of the channel. However, in this work, we demonstrate that a fluid–fluid interface can, in fact, be partially stabilized if the displacing fluid is sufficiently viscous. We extend the analysis to account for buoyancy effects in an inclined system, and show that if there is a sufficiently large density contrast between the fluids then the interface may be stabilized or destabilized by the buoyancy force. In summary, the model provides a series of idealized solutions for the development of a topographically controlled viscous finger.

Once we have built the model, we illustrate the phenomena with a series of simple analogue experiments in which one fluid is displaced by a second immiscible fluid in a long, narrow channel in which the gap thickness varies in the direction perpendicular to the flow. Finally, we apply the model to provide insight into the displacement pattern of CO_{2} or polymer-laden water injected into a channelled reservoir, inferring estimates for the fraction of the channel that may be invaded by the advancing fluid and hence providing insights into the storage capacity of the system for CO_{2} sequestration and the sweep efficiency of a polymer flood for enhanced oil recovery.

## 2. Model for a topographically controlled viscous finger

We consider the injection of fluid into a long, relatively narrow channel, of small but finite thickness, which is confined within an impermeable seal rock. We assume the injected fluid displaces the original fluid down the channel, as sketched in figure 2, and that the flow extends far down the channel. In practice, the channel may be sinuous, but we work with the along-channel coordinate, *x* say, and assume that the wavelength of the sinuous fluctuations is long compared with the amplitude of these fluctuations.

### (a) Model equations

We consider the flow in a long channel, 0<*x*<*L*, of thickness *b*(*y*) and width −*d*<*y*<*d*, where *b*(*y*)≪*d*≪*L*, so that the pressure gradient is nearly uniform across the thickness of the channel, and the flow is essentially two dimensional, varying with position across (*y*) and along (*x*) the channel (*u*(*x*,*y*), *v*(*x*,*y*)) (figure 2). We also assume that the channel permeability, *k*(*y*), and thickness, *b*(*y*), decrease smoothly from the centre to the edges of the channel. In the original reservoir fluid, of viscosity *μ*_{o}, the depth-averaged two-dimensional flow follows Darcy’s Law,
2.1in terms of the pressure gradient along and across the channel where *k*_{o} is the relative permeability of the original fluid in the formation. Continuity requires that
2.2Given that *L*≫*d*, equations (2.1) and (2.2) imply that *u*/*v*=*O*(*L*/*d*)≫1. We infer that to leading order the flow is parallel to the long axis of the channel and the pressure gradient across the channel is much smaller than that along the channel.

We now consider the injection of a second fluid, *i*, with viscosity *μ*_{i}, into the cell to displace the original fluid. If the fluids are immiscible, then the surface tension, *σ*, will lead to a pressure jump between the two fluids given by a capillary pressure function *p*_{c}(*σ*,*δ*), where *δ* is a function of the grain size and order [1,3]. On the pore scale, two-phase flow theory suggests that there will be a fractional flow of the displacing and the displaced fluid, but as the flow develops this leads to formation of a saturation jump, across which there is a transition from flow of the original fluid to a mixed flow of the injected and original fluid. Behind this saturation jump, there is typically an increasing flux of injected fluid, and a decreasing flux of the original fluid, until a point upstream where the original fluid reaches residual saturation, and all the flow is composed of injected fluid [1,3]. In the present analysis, for simplicity, we assume that downstream of the saturation jump all the flow consists of the original reservoir fluid, whereas upstream all the flow consists of injected fluid with a fraction of the pore space being host to the residual saturation of the original fluid. We now aim to model the advance of this saturation front on the scale of the channel. Strictly the effective permeabilities seen by the injected fluid upstream and the original fluid downstream of the saturation jump equal the product of the rock permeability *k* and the relative permeability of each phase, *k*_{o} and *k*_{i}, and this leads to the mobility ratio given by *V* =*μ*_{i}*k*_{o}/*μ*_{o}*k*_{i}. Provided the effective capillary number is small, *kp*_{c}/*uμ*_{o}*d*≪1, where *p*_{c}/*d* is the capillary pressure gradient scaled with the width of the channel, then this front across which there is the saturation jump will remain localized in the cross-channel direction [5].

We assume that the displacing fluid advances through the centre of the channel, displacing the original fluid as a single topography-controlled finger, of width *w*(*x*,*t*). In general, we expect this finger to become progressively narrower with distance down the channel, especially if the displacing fluid is less viscous than the original fluid in the formation. Assuming there is a well-defined interface such that the displacing fluid is located in the region 0<*y*<*w*(*x*,*t*), we now seek a relation governing the evolution of *w*(*x*,*t*). First, the mass flux of the displacing and the original fluid are given by
2.3and
2.4where
2.5If the width *w*(*x*,*t*) of the displacing fluid varies slowly along the channel such that ∂*w*/∂*x*≪*w*/*d*, then, locally, the pressure in the original fluid may be related to the pressure in the displacing fluid by the approximate relation
2.6

At each point along the channel, the total mass flux along the channel, *Q*, may be expressed in the form
2.7where *V* is the mobility ratio of the original to injected fluid, which is linearly proportional to the viscosity ratio of the injected to the original fluid. The local conservation of displacing fluid may be expressed by the integral expression
2.8Combining this with equation (2.6) to eliminate the pressure gradient, we find the equation for the evolution of *w*(*x*,*t*), the width of the displacing fluid at distance *x* along the channel,
2.9where *F*(*w*), the fractional flow of the displacing fluid along the channel in the absence of surface tension, is given by
2.10Equation (2.9) is an advection–diffusion equation in which the advection speed of surfaces of constant *w* is given by
2.11This speed varies with *w*, leading to a tendency to focus the flow along the central part of the channel where the permeability is expected to be the greatest. The capillary pressure, which manifests itself in the expression on the right-hand side of equation (2.9), is diffusive in character. If the interface, *w*(*x*,*t*), becomes stretched in the along-channel direction, the surface tension term becomes progressively smaller, and to leading order the evolution of the interface will be described by the advection speed (equation (2.11)), except in a small region near the nose where the interfacial tension will smoothly connect the solution in the region *y*<0 with the solution in the region *y*>0. Equation (2.10) is reminiscent of the equation for Buckley–Leverett theory for the migration of saturation surfaces in two-phase flow based on analysis of the pore-scale flows, but the present theory applies to the macroscopic flow on the scale of the channel.

### (b) Shock-like solutions

Provided that (1/*b*)(d*F*/d*w*) varies monotonically with *w* then the interface will stretch smoothly in time and the evolution of *w*(*x*,*t*) is of a dispersive nature. However, if (1/*b*)(d*F*/d*w*) has a local maximum, then the model predicts that at some point *x* along the channel the interface position relative to the centre of the channel will become multi-valued. We expect that such multi-valued behaviour in the speed of the saturation wave can occur if (i) the displacing fluid is sufficiently viscous that the rate of the increase in the fractional flux with width of the displacing fluid initially increases as the width of the displacing fluid increases from zero, but provided that (ii) the thickness or permeability of the outer part of the channel is sufficiently small that the fractional flux only converges towards unity slowly as the width of the displacing fluid approaches the width of the channel. In this case, we expect that the interface will develop an inner region, 0<*w*<*w**, which advances with constant speed, and for interface widths *w*>*w** the speed gradually decreases with interface width, leading to an outer region where the interface gradually spreads out and disperses. We note that in the nose region, 0<*w*<*w**, moving with constant speed, there will be a detailed readjustment of the flow just behind the front. This process is somewhat analogous to the phenomena of shock formation in a capillary tube when fluid of one viscosity displaces a second fluid of different viscosity [31].

The location down the channel, *x**(*t*), of the point at which *w*=*w** may be found by noting that (*a*) for mass conservation
2.12whereas (*b*) the continuity of the speed of the interface requires that equation (2.11) applies at *x*=*x**
2.13Combining (2.12) and (2.13), we see that *w** is given by the solution of the relation
2.14

### (c) Illustration of the flow structure

To make progress, it is useful to have a specific form for the channel thickness as a function of the position relative to the centre of the channel. For convenience, we now work with the dimensionless channel position, , with , and the dimensionless along-channel position . Given the decrease in thickness and permeability away from the axis of the channel, expected from geological deposits [20], as a simple parametric model, we take
2.15corresponding to a gap of width *b*_{o} at , which then decreases to width *b*_{o}(1−λ)^{α} at . We also set the permeability to decay with width as
2.16Substituting these functional forms into equation (2.11) leads to the relation
2.17where the dimensionless time is defined as *τ*=*Qt*/*Ldb*_{o}.

In figure 3, we illustrate the variation of the speed of the interface as a function of the cross-channel width of the displacing fluid, *w*, for three values of *V* in the case that the displacing fluid is (*a*) more mobile, as given by *V* =1.00 (green), 0.75 (red) and 0.50 (blue), and (*b*) less mobile, as given by *V* =1.25 (blue), 2.00 (green) and 5.00 (red). In both cases, in this example, we assume *α*=1, *γ*=2 and λ=1. In figure 3*a*, it is seen how the finger of fluid which advances along the centre of the channel becomes increasingly focused as the mobility of the injected fluid increases above that of the original fluid in the channel (*V* <1). We describe this as topographically controlled fingering, and this is of key importance for modelling the migration of CO_{2} along a channel (§5). In figure 3*b*, it is seen that when the mobility of the injected fluid decreases below that of the original fluid in the aquifer, *V* >1, the speed of the interface decreases. With a mobility ratio *V* =1.25, the interface speed still decreases monotonically with distance from the centre of the channel and so this again will lead to a dominant focusing of the flow towards the centre of the channel, but the flow now becomes somewhat broader. This illustrates the importance of the topography of the channel in regulating the flow focusing. However, as the mobility of the injected fluid decreases, so that *V* increases to value 2.00, a maximum in the speed of the interface, (*Q*/*b*) d*F*/d*w*, is predicted to develop away from the centre of the channel (figure 3). In §2b, we proposed that this leads to the formation of a shock near the centre of the channel, which advances with uniform speed in the region , whereas a dispersed interface spreads out along the channel in the region beyond. A similar current structure develops in the case *V* =5.00, and in figure 3*c* we illustrate the location of the shock (red line) in the centre of the channel, and indicate the region where there is a dispersing tail in the outer parts of the channel. In this calculation, we have taken λ=1, so that the thickness of the channel falls to zero when .

In figure 4*a*,*b*, we illustrate the evolution of the shape of the interface as a function of time according to this simplified advection model, for the cases (i) *V* =0.75 and (ii) *V* =5.00 corresponding to (i) a very unstable flow and (ii) a flow with a strong viscous stabilization effect. In the first case, a topographic viscous finger develops along the axis of the channel, as time increases, whereas in the latter case the blunt nose structure develops in the centre of the channel as expected.

In the case that a finger develops the pure advection model predicts that there is a discontinuity in the gradient of the fluid–fluid interface at the axis of the channel. The capillary pressure will smoothly connect the solution in the region with that in the region . From a scaling perspective, equation (2.9) illustrates that the length scale of this capillary adjustment region at the nose of the flow is given by
2.18and, in order that the present modelling applies, we require that Δ*x*≪*d*, so that the adjustment occurs in a local region near the nose of the flow (cf. §2a). Upstream of this adjustment region in the centre of the channel, we expect the flow to follow the model described by figure 4*a*.

### (d) Condition for a shock

Given the parametric model for the channel thickness and permeability, it is possible to develop some expressions for the onset of formation of a shock, and also for the condition when the shock fills the whole width of the channel. These are useful predictions, because they indicate how large the viscosity of an injected fluid should be in order to stabilize a topographically induced fingering flow. From the above analysis, we recognized that a shock-type solution forms if the right-hand side of equation (2.17) has a maximum in the region . If we differentiate expression (2.17), then it follows that there is a maximum if there is a root of the relation
2.19for some in the range . Because increases with , then a root of equation (2.19) exists if *V* *(0)<*V* <*V* *(1). In the limit , equation (2.19) leads to the result that
2.20Furthermore, in the case , we find that
2.21

As a specific example, we consider the limit λ=1 in which the gap width falls to zero at the outer edge of the channel. We also impose the conditions that *α*=1 and *γ*=2, which also corresponds to flow in a Hele-Shaw cell with a linearly decreasing gap width as a function of the cross-channel position. In this limit, we find that in order that there is a maximum wave speed in the region , and hence the formation of a shock, we require that , because . We also find in this case that so that there is no mobility contrast which is sufficient that the shock extends to the outer edge of the channel, and so the flow will also lead to a dispersive tail, in the absence of capillary effects. More generally, if λ<1, then *V* *(1) remains finite, and so for *V* >*V* *(1) the shock is predicted to extend across the whole width of the channel.

## 3. Effects of buoyancy

In the event that, if there is a density difference between the original fluid in the channel and the injected fluid, there is a component of gravity along the main axis of the channel, say, where *θ* is the downward angle of tilt to the horizontal, the motion of the interface will also be influenced by the buoyancy force, according to the relations
3.1and
3.2By including these additional gravitational terms in the model, then, in the absence of capillary pressure effects, the governing equation becomes
3.3where
3.4corresponds to the ratio of the buoyancy-driven flow along the channel to the supply flux *Q* and Δ*ρ*=*ρ*_{i}−*ρ*_{o}. The speed of surfaces of constant *w* is now given by
3.5

If the buoyancy force is stabilizing, this may lead to the formation of a shock in the centre of the channel, thereby acting to suppress the focusing of the flow along the centre of the channel (figure 5*a*); in contrast, if the buoyancy is destabilizing, it may accelerate the focusing of the interface along the centre of the channel by increasing the flow speed in the centre of the channel, as seen in figure 5*b*.

## 4. An exploratory experimental model

We have carried out a series of very simple exploratory experiments in a long thin Hele-Shaw type of cell composed of a series of plates of Perspex that form a gap of variable thickness in the cross-channel direction. The channel was 15 cm wide, 180 cm long, and the gap increased from 0 to 3 mm from each edge to the centre of the channel. The cell was initially filled with vegetable oil, and this was displaced by supplying a mixture of water and glycerol at a constant flow rate. The end of the cell in which the source fluid was added was sealed, except for a series of injection ports that were uniformly distributed across the width of the channel, as in a simple model of a line well. The other end of the cell was open, to allow fluid to spill out. During the setting-up phase of each experiment, the tank had a small inclination to prevent draining of fluid from the tank prior to the start of the experiment. Just before the pump was turned on to start an experiment, the inclination of the tank was set to a range of values from being horizontal to being vertical. We note, for reference, that, in a Hele-Shaw cell, the permeability is given in terms of the gap thickness *b* as *b*^{2}/12, and so the Hele-Shaw cell of variable width represents a good analogue of a channel of variable permeability in the cross-channel direction.

In the experimental system, *Q* was set to values of order 1–10 ml s^{−1}, while (*ρ*_{glycerol}−*ρ*_{oil})≈200 kg m^{−3}, so that the dimensionless parameter, which compares the gravitational and forced flow speeds (cf. equation (3.3)), is given by *Gk*(0)*b*(0)≈1−10. This shows that, with the cell tilted upwards, the gravitational stabilization of the front associated with the greater density of glycerol compared with the oil is of comparable magnitude to the destabilization associated with the applied flux *Q*. We now present some results with different angles of slope and different viscosity ratios. The surface tension across the interface, *σ*, leads to an effective capillary pressure. With *σ*∼40–50 Nm^{−2}, the ratio of the flow associated with the surface tension gradient and the supply flux, as given by (cf. equation (2.9)) *σb*^{2}(0)*w*(0)*d*/12*μQ* is very small, of order 10^{−3} to 10^{−4} for these flow rates.

Figure 6 presents a series of results from three experiments in which a water–glycerol mixture displaced oil in the Hele-Shaw cell, for cases with (i) viscosity ratio 2 : 1 and with angle of slope 10°; (ii) viscosity ratio 2 : 1 and with angle of slope 58°; and (iii) viscosity ratio 10 : 1 and with the cell being horizontal. In each case, in figure 6*a*, we show a photograph of the red finger of glycerol displacing the clear oil. Figure 6*b* shows a series of lines corresponding to the position of the interface at a series of times during the experiment, to illustrate the time evolution of the front. Figure 6*c* illustrates the superposition of these individual interface profiles, , but with each rescaled by the time at which it was taken after the start of the experiment, . The interface profiles suggest that with the smallest viscosity ratio, case (i) and a small angle of tilt (10°), a narrow finger migrates along the central part of the channel. As the cell is tilted to an angle of 58°, the gravitational force stabilizes the interface, leading to a much flatter interface (figure 6*b*), whereas with the largest viscosity ratio, case (iii), the fluid nearly fills the width of the cell, as would be anticipated from the model, because the viscosity ratio suppresses the tendency for the flow to focus in the centre of the channel. In each experiment, the rescaled data showing the cross-channel position of the interface as a function of the scaled along-channel position suggest that the data collapse to a universal curve that depends on the viscosity ratio and angle of tilt.

In figure 7, we compare the experimental results shown in figure 6 with the predictions of the model accounting for both the capillary pressure and also the gravitational effects (cf. equations (2.9) and (3.3)). To this end, we have solved the model equation
4.1numerically using a finite-difference code in which we adopt an upwind differencing scheme for the advection terms and a central differencing scheme for the diffusion-type terms [32]. The advection terms, which correspond to the second term on the left-hand side of the equation, are associated with the applied flow *Q* and the buoyancy force, whereas the diffusion-type term, on the right-hand side of the equation, is associated with the surface tension, *p*_{c}=*σ*/*b*. It is seen that in each case there is reasonable agreement between the model prediction and the experimental observation of the width of the advancing finger of glycerol except for some discrepancy near the head of the finger. Near the front of the advancing flow of glycerol, there is an adjustment from the velocity profile in the advancing current of glycerol behind the head to the flow in the head. This involves a two-dimensional adjustment zone that has not been described in the present model. In order to visualize the adjustment of the flow profile in the main current and in the head, in figure 8, we present four photographs at successive times for the case *V* =5.00 in which we illustrate the motion of four lines of dye advancing along the channel. These lines of dye were released at successive times in the experiment, and it can be seen how the flow is more rapid in the centre of the cell, as anticipated by the theory. It may also be seen how the fluid adjusts in the near-head region. It would be interesting to develop a two-dimensional model for the structure of this flow in the head region, but this is beyond the scope of the present investigation. Nonetheless, the present model provides a reasonable representation of the evolving shape of the interface behind the head region. We also note that any deformation of the cell walls associated with the pressure of the flow could lead to some adjustment in the position and shape of the front, and we plan to explore such interactions in more detail in future work.

## 5. Applications

The main purpose of this work has been to explore how variations in the cross-channel permeability and thickness influence the development of the flow along the channel. We have shown that the flow within an individual channel is very strongly controlled by the boundaries and also the variation in permeability across the channel. We have developed a simple model that leads to some idealized analytical solutions for the shape of the nonlinear fingers as they migrate down the channel, based on a long-wave approximation for the flow. When the injected fluid is less viscous than the original fluid in the channel, a localized topographic focusing of the flow along the centre of the channel develops. In the opposite limit relating to the injection of a relatively viscous fluid, the stabilizing effect of the viscosity contrast tends to suppress the topographic channelling of the flow. Now, the simple long-wave model predicts the development of a multi-valued interface shape that leads to formation of a shock-like adjustment at the head of the flow and the formation of a relatively blunt nose. Although the models are highly idealized, and focus in particular on the sweep efficiency of a flood front, they can be used to provide insights into two problems of interest, namely the storage potential of a channel system as part of a CO_{2} sequestration project and the effectiveness of polymer flooding for enhanced oil recovery.

### (a) CO_{2} sequestration

With CO_{2} sequestration, supercritical CO_{2} is typically injected to depths of 1.5–2.0 km below the surface. In such conditions, we expect that the injected fluid is less dense and less viscous than the original water in the aquifer, with a viscosity ratio in the range of 1 : 10–1 : 20 and a density ratio of order 0.6–0.9 [11,17,33]. As a result, we envisage that, if the CO_{2} is injected into a channel system, both the viscosity ratio and the buoyancy ratio may lead to flow focusing along the centre of a channel, and so the flow will bypass the outer parts of the channels. This may lead to an inefficient sweep of the channel by the injected CO_{2}, limiting the storage potential of the system. To illustrate the idealized sweep efficiency of an injection, in figure 9*a*, we illustrate the fractional flow of CO_{2} at a point *L* downstream of the injection point (cf. equation (2.11)) and also the fraction of the volume of the channelled reservoir that has been swept by the injected fluid, as a function of the time after the start of injection. In figure 9, time is measured in terms of the number of pore volumes in that region of the channel located between the injection well and the point a distance *L* downstream which have been injected into the formation (figure 9*a*). For the calculation, we assume that the channel has a geometry given by the idealized model with *α*=1 and *γ*=2 and λ=1. Illustrative results are given for fluid with mobility ratios of 0.05 (blue), 0.1 (red) and 0.25 (black). The dashed lines represent the fraction of the pore space that has been swept as a function of time, whereas the solid lines represent the fraction of the injected CO_{2} flux that passes the point *L* downstream of the injection point as a function of time. It is seen that the CO_{2} breaks through to the point a distance *L* downstream very rapidly (blue solid line), because the flow tends to migrate through a very localized region in the centre of the channel. As a result, only a very small fraction of the reservoir volume (<0.1) is swept by the injected fluid (blue dashed line). As the injection process continues, there is only a slow increase in the fraction of the reservoir pore space swept by the CO_{2}. The model illustrates the challenge of accessing a large localized volume of pore space through injection of CO_{2} for sequestration. Because the CO_{2} appears to short-circuit much of the pore space in the channel, the strategy for sequestration of a large volume of fluid would then benefit, for example, from there being a structural trap downstream of the channel [34].

### (b) Enhanced oil recovery

For the process of enhanced oil recovery through injection of polymer-laden water, the viscosity of the water may be increased through addition of a range of polymers, with the purpose to increase the fraction of the reservoir swept by the flow. The model provides us with an idealized approach to exploring the impact of different viscosities of the polymer solution on this sweep efficiency. In particular, our model suggests that, as the viscosity of the injected fluid increases above a critical value in excess of the viscosity of the original fluid, there is a transition to a shock-type structure that forms a blunt nose near the leading edge of the fluid–fluid front. We have solved the relationship for the width of this shock as a function of the flow rate, for the case *α*=1, *γ*=2, and in figure 9*b* we illustrate how this varies with the ratio of the viscosity of the injected fluid compared with the original oil in place. As the viscosity ratio increases to values of 2–3, there is an initial significant uplift in the swept volume of oil from the system. However, there is a law of diminishing returns for larger polymer concentrations in that the increase in the size of the shock becomes more gradual with further increases in viscosity ratio. In addition with substantially larger values of the viscosity of the injected fluid, the mechanical power required to inject progressively more viscous treatment fluid increases and eventually the large injection pressures that may be required have the potential to damage the formation.

## 6. Conclusion

We have developed a model for the migration of CO_{2} or polymer-laden water injected into a channelled reservoir, illustrating how the shape of the fluid–fluid interface evolves with time. The model is based on an integral approach for the flow along the channel, and uses the typical separation of scales between the thickness, the width and the length of the channel to generate an analytical description of the flow. It is shown how the viscosity ratio of the fluids and the magnitude of any buoyancy forces can substantially change the sweep efficiency of the flow. In particular, we illustrate how the topographic focusing of the flow along the high-permeability central axis of the channel can dominate, leading to a very poor sweep of the system. Even if the injected fluid is of the same viscosity as the original reservoir fluid, the flow tends to focus as it migrates along the reservoir.

Application of the model suggests that, owing to the relatively low viscosity of CO_{2}, the sweep efficiency of a CO_{2} injection may be very low, so that the CO_{2} may access only a very small fraction of the pore space. Such modelling may be valuable for assessing the overall storage potential of a site. The model has also shown that the use of a viscosifier in the injected water used in enhanced oil recovery can substantially enhance the sweep along a channel, especially if the injected fluid is more viscous than the oil in place, with the viscosity contrast suppressing the natural tendency for the flow to focus along the high-permeability central axis of the channel.

In planning application of these technologies in a real system, there is a significant element of uncertainty in the structure of the formation and the detailed properties of the rock. The models we have presented in this work do not account for the details of the pore-scale displacement processes, but assume that there is a well-defined interface between the original and the injected fluid. This enables us to focus on the macroscopic sweep efficiency of the flow, and has identified how the geometry of the formation and the structure of the permeability within such a formation can be key for determining the fraction of the pore space accessible to CO_{2} during injection, as part of a sequestration project. Analytical models, such as those presented herein, can also be used to inform assessments of the sensitivity of model predictions to uncertainties in the geological description of the reservoir.

In conclusion, we note that the model provides an interesting new lens on the process of viscous instability; in particular, the development of a simplified model of the nonlinear growth of the ‘topographic’ viscous finger in a channel of non-uniform thickness, through use of long-wave approximations, complements the body of work on modelling viscous instability in a Hele-Shaw cell of constant thickness. We plan to report on a more detailed suite of experimental models of the process in the near future, but the initial experiments presented in §4 illustrate the general trends identified in the modelling.

## Data accessibility

Data presented in this paper may be obtained from Andrew Woods via email: andy{at}bpi.cam.ac.uk.

## Competing interests

The authors declare that they have no competing interests.

## Funding

We are grateful to bp for funding.

## Authors' contributions

A.W.W. carried out all the mathematical analysis, wrote the paper and carried out the comparison of the experimental data with the model. N.M. carried out the experiments. Both authors read and approved the manuscript.

## Footnotes

One contribution of 12 to a theme issue ‘Energy and the subsurface’.

- Accepted June 10, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.