Pressure-driven growth in strongly heterogeneous systems

The pressure-driven growth model for advance of a foam front through an oil reservoir during foam improved oil recovery is considered: specifically the limit of strong heterogeneity in the reservoir permeability is treated, such that permeability variation with depth more than outweighs the tendency of the net pressure driving the front to decay with depth. This means that the fastest moving part of the front is not at the top of the solution domain, but rather somewhere in the interior. Moreover the location of the foam front on the top boundary of the system can no longer be specified as a boundary condition, but instead must be determined as part of the solution of the problem. Numerical solutions obtained from the pressure-driven growth model under these circumstances are compared with approximate analytic solutions. An early-time approximate solution is found to break down remarkably quickly (far more quickly than breakdown would occur in the analogous homogeneous system). Numerical solutions agree much better with local quasi-static solutions centred about local maxima in the front shape, each local maximum corresponding to a depth within the reservoir at which a high permeability stratum is found. These individual local solutions meet together at sharp concave corners to cover the entire depth of the foam front. As time continues to progress however, the system evolves towards a long-time, global quasi-static solution, corresponding to the fastest moving of the aforementioned local maxima. Additional key features of the predicted front shapes are elucidated. The foam front is found to meet the top boundary obliquely despite an established convention in pressure-driven growth that the front and top boundary should meet at right angles. In addition, at each sharp concave corner, discontinuous jumps are predicted in the path length that material points travel to reach either side of the corner. Moreover the long-time, global quasi-static solution is found to admit smooth concavities, as opposed to the aforementioned sharp concave corners, which only tend to be prominent earlier on.


Introduction
During oil production, only a fraction of the oil originally in place within an oil reservoir is driven out under the reservoir's own pressure. Once the reservoir's pressure is depleted, additional oil can be extracted by injecting driving fluids into the reservoir. Amongst the candidate fluids for driving out the oil are gas-liquid foams [1][2][3], employed in the process of "foam improved oil recovery".
The reasons why foam is a good candidate driving fluid for improved oil recovery stem from a number of properties of foam, particularly the properties it exhibits when present within porous media such as oil reservoirs [4,5]. Using foam as a driving fluid, rather than using injected gas alone, tends to prevent gas (within the foam) from simply rising to the top of the reservoir and "overriding" the oil present [6,7]. Moreover since foam is a fluid with complex rheology, such that the tendency of foam films either a e-mail: paul.grassia@strath.ac.uk to block or to flow through pores [8] can be sensitive to pore size, foam can potentially reach oil within pores that other injection fluids might not have reached. Yet another key property of foam is that its mobility within porous media is surprisingly low, particularly when the foam is finely-textured [6]. This low mobility enables a foam front within a porous medium to advance in a uniform fashion (rather than developing fingering instabilities [9], which would involve flow then occurring upon preferential flow paths, namely the fingers, but bypassing the oil still in place).
Multiphase flow in porous media is classically described by Darcy-type relationships [10], but with different fluids in the multiphase flow having different transport properties [8,11] (e.g. different relative permeabilities and/or different effective viscosities, with the relative mobility being the ratio between relative permeability and effective viscosity). Recognising, however, that foams can be much less mobile than other fluids in the reservoir,  Fig. 1. A front of finely-textured foam, which is initially vertical but which moves from left to right as time evolves. The front (shown here as a thick curve) separates liquid (surfactant solution and oil) downstream (towards the right) from coarse foam upstream (towards the left). The front can be described at any time t by a curve x vs. y, with each front material point having displaced through a distance s to reach its current location. The instantaneous direction of motion of front material points is along the front normal n, with α denoting the angle between the normal and the horizontal. For a system with strongly heterogeneous permeability, x can exhibit several local maxima, x max,1, xmax,2, xmax, 3, whilst x at the top (denoted x top) might lag behind one or more of these. Sharp concave corners can appear between the various maxima.
Darcy's law can be simplified to a case in which essentially all the pressure drop in the flow is realised across the least mobile fluid present [12].
Specifically, we shall consider the case of foam production in situ within the reservoir, with a slug of aqueous surfactant solution injected first, followed by a slug of injected gas: this process is known as surfactant alternating gas [13,14]. A foam front is formed at the interface between surfactant and gas and this interface migrates through the reservoir as injection proceeds. Under these circumstances, the finest textured foam (i.e., the smallest bubbles) is found at the front itself, whereas upstream the foam is much coarser (since bubbles formed upstream have already collapsed [15][16][17]), and downstream there is no foam (only surfactant and oil). The least mobile part of the flow by far is the finely-textured foam at the front [12], so as we have noted earlier, this is where the bulk of the Darcy pressure drop must be realised.
In typical applications, this finely-textured foam front is a region which is perhaps tens or hundreds of metres thick containing a myriad of foam bubbles [18]. Nonetheless we are dealing with a front that can itself propagate over a distance of kilometres. Relative to the distance scale over which the front propagates, the front thickness is considered to be comparatively small. To a reasonable approximation then, the foam front can be represented (see, e.g. fig. 1) by a 1-D curve in a 2-D domain (or more generally by a 2-D curved surface in a 3-D domain [12,18], although we do not treat that generalisation here). If the entire Darcy pressure drop is realised across this effectively 1-D curve, then Darcy's law reduces to a simpler mathematical model known as "pressure-driven growth" (see [12,18,19]). Following a number of recent studies in the literature [20][21][22][23][24][25], this is the model we will analyse here.
Pressure-driven growth solely involves computing how the locus of the 1-D foam front evolves over time: it does not require a solution of the full 2-D domain. According to the model [12,18], any given element of the 1-D front advances due to the net pressure driving it, the net pressure here being the difference between the injection pressure (behind the front) and a hydrostatic pressure (ahead of it). The hydrostatic pressure increases with depth, meaning that the net driving pressure decays with depth: eventually one reaches a maximum depth at which hydrostatic pressure balances injection pressure, and the front cannot advance beyond this. At depths less than the maximum, the net driving pressure is balanced by the dissipation associated with moving along the finely-textured, low mobility foam front [12,18]: this balance sets the speed of each element of foam front, with the direction of motion of each element being normal to the front.
Although conceptually the pressure-driven growth model is relatively easy to explain, it admits a surprisingly complex set of mathematical behaviours. The reason why the mathematical behaviour manages to be so rich is that the equation governing the model turns out to be a hyperbolic partial differential equation, and it is well known that such equations admit solutions that are not necessarily smooth [26]. In other words, pressure-driven growth, as formulated, lacks any regularising "diffusive"type term that would require the solutions to be smooth. The most striking manifestation of "non-smooth" solutions in pressure-driven growth is that concave front shapes can potentially focus down to sharp concave corners: again see fig. 1. Physically what these sharp corners represent is that the foam front can reorient itself on a length scale comparable with the thickness of the finelytextured foam front (rather than on a length scale comparable with the much longer distance over which the front propagates).
In spite of the tendency of the model to develop solutions that are obviously not smooth, for many problems of interest, including the advance of a foam front through a homogeneous reservoir, these sharp concave corners tend not to manifest [12,18]. Instead the front tends to be convex (or strictly speaking so close to being convex that tiny departures from convexity can be easily overlooked [25]): convex shapes do not lead to the same mathematical illbehaviour in the pressure-driven growth model that concave shapes do [18].
Nonetheless there are important classes of problems in which concavities can and do arise [20,21,23]. Amongst these (see [23]) is the case of a heterogeneous reservoir in which the permeability of the porous medium varies with depth. The reservoir might for instance be stratified, with different strata in the reservoir at different depths. Concavities in the front shape then arise naturally in those strata for which the front motion lags behind the motion in neighbouring strata, and as [23] found, these concavities can indeed focus down into sharp corners, such as those shown schematically in fig. 1.
That prior study [23] however left a number of questions still unanswered. One unanswered question was whether there might be a simple way to represent the front shape away from the concave corners. Another unanswered question was how the front shape evolved at long times, and moreover whether or not the concave corners persisted at these long times. Another key unanswered question concerned the behaviour of cases with comparatively strong heterogeneity, and in particular how foam fronts in such cases behave near the top boundary of the reservoir. We explain the issue near the top boundary as follows.
We have already stated that the net driving pressure difference (injection pressure less hydrostatic) is highest at the top [18]. Ordinarily this leads to the top of the front advancing further than the remainder of the front lower down. Typically moreover, a boundary condition is imposed such that the tangent to the front meets the top boundary of the solution domain at right angles. Since the direction of motion of front material points is normal to the front itself, this boundary condition serves to keep the top of the front moving along the top boundary of the domain, as one might expect [12]. As we will explain in more detail later on, a consequence of this boundary condition turns out to be that the location of the top of the front is known a priori for all times [18], and this can be imposed on the solution for the front shape. A potential issue however arises in a heterogeneous reservoir for which the reservoir permeability happens to be decreasing on the approach to the top boundary. The permeability affects the rate of advance of the foam front, and if the decrease in permeability near the top is sufficiently strong, it can overcome the tendency of material points higher up to have a higher rate of advance owing to their higher net driving pressure. The system is then said to be "strongly heterogeneous", and in effect, front material points on the top boundary could lag behind points slightly lower down: see, e.g., fig. 1. This makes it impossible, at least within the pressure-driven growth model, to impose a condition that the front tangent meets the top boundary at right angles. Failure to impose this particular boundary condition also implies that the location of the top of the front is no longer specified as part of the problem formulation but instead is a priori unknown. In the language of hyperbolic partial differential equations, it is possible to show that the trajectories of front material points correspond to characteristic curves with information being propagated along the characteristics [18]: instead of characteristics entering the solution domain from the top boundary, both they and the information that propagates along them are now being directed outwards from the solution domain towards that boundary. We emphasise that this situation only occurs when the heterogeneity in the permeability is sufficiently strong as to outweigh the depth variation of the net driving pressure. When the heterogeneity is weaker, the usual top boundary condition applies, as per what occurs in homogeneous systems.
The purpose of the present work is to address some of the issues and unanswered questions that were identified above specifically for pressure-driven growth in systems having heterogeneous permeability, with a focus on strongly heterogeneous systems in particular. The remainder of this work is laid out as follows. Section 2 outlines the governing equations of pressure-driven growth, firstly for a homogeneous system, and subsequently extended for a (strongly) heterogeneous one. Numerical solution techniques are discussed in sect. 3, with a number of approximate analytical solutions outlined in sect. 4. Results are presented in sect. 5, which also compares data from numerical simulations with approximate analytical predictions. Conclusions are offered in sect. 6.

Model and governing equations
The formulation of the pressure-driven growth model can be found in several publications that are already in the literature (see, e.g., [12,18,20,21,23]). Hence we will not derive the model again in full here, opting instead merely to introduce the governing equations and describe what they mean physically.
The key governing equation is Here x represents the location of a material point on the foam front, t represents time, k represents the reservoir permeability (assumed to be homogeneous in the first instance), λ r is the relative mobility of the finely-textured foam front, S w is the liquid saturation in the foam (i.e., the volume fraction of aqueous surfactant solution in the foam), φ is the reservoir porosity, ΔP is the net driving pressure (injection pressure less hydrostatic), τ thick is the front thickness, and n is the front normal. That this is a hyperbolic partial differential equation can be appreciated since on the left-hand side we have a time derivative of x, whereas on the right-hand side we have the front normal n, which is simply a rotation of the front tangent through a right angle, with the front tangent itself being a derivative of x with respect to distance along the front. The net driving pressure ΔP is the difference between a driving injection pressure P drive and a hydrostatic pressure, the latter increasing with depth. If ρ is the liquidto-gas density difference and g is acceleration due to gravity, there is a maximum depth to which foam can reach, d max ≡ P drive /(ρg) injection pressure being balanced to hydrostatic pressure at this particular depth. Furthermore, although the front thickness τ thick is much smaller than the distance over which the front propagates, it turns out that the front thickness is proportional to the distance over which front material elements displace. This is expressed as τ thick = τ s where s is the distance through which a front material point has displaced, and τ is a dimensionless constant much smaller than unity. Typically (see [18]) we choose τ = 0.01, although the analysis we present works for any τ value, provided τ 1.
We define an (x, y) coordinate system, such that x = 0 is the location of an injection well and y = 0 is the aforementioned depth at which injection pressure and hydrostatic pressure balance. We make the governing equations dimensionless by scaling distances by d max and times by a scale Values of d max and t scale are of course sensitive to the parameter values we set, but typical values reported by [20] suggest d max values from several hundred metres up to a couple of kilometres, and t scale roughly on the order of a week. Continuing in what follows to employ the symbols x ≡ (x, y), s and t, to represent dimensionless variables (the entire discussion from here on being in dimensionless form), the governing equation becomes where in addition ds/dt = n · dx/dt = y/s.
Here the term in y on the right-hand side represents the net driving pressure: this vanishes at the bottom of the domain, and is highest at the top (y = 1 in our dimensionless system). This means that points higher up on the front tend to move faster and further than those lower down. The initial condition is ideally x = 0 and s = 0 at t = 0 for all values of y. This however leads to arbitrarily large initial values of dx/dt and ds/dt. To avoid this, in numerical computations we tend to set x = s = s 0 at t = 0 where s 0 is a small dimensionless parameter, typically s 0 = 0.01. As we will discuss later, the solutions that develop over time are really very insensitive to s 0 , so the exact choice of s 0 has little bearing on the model predictions.
A boundary condition is also usually imposed on the top boundary. To achieve this we first define an angle α such that the subscript "f " here representing a derivative along the foam front. According to this definition α is the angle that the front normal makes to the horizontal; the sign convention is that a positive α corresponds to a front normal pointing below the horizontal (see fig. 1), which is what we expect to develop if upper parts of the front move further and faster than parts lower down.
Having now defined α, the usual boundary condition that is imposed is α = 0 at the top. This ensures that the front normal is along the top boundary, and hence (according to eq. (3)) the direction of motion of front material points is along that boundary.
If points are moving along the top boundary, the distance through which they have displaced s must be the same as their horizontal coordinate x, which we now denote by x top . Thus for the top boundary y = 1, eq. (3) reduces to dx top /dt = 1/x top of which a solution is More generally the solution can be written although the difference between eqs. (6) and (7) is usually negligible: for the value s 0 = 0.01 mentioned earlier, the difference is negligible for any t value greater than order 10 −4 . All of the above development has envisaged a homogeneous permeability. A heterogeneous permeability can easily be incorporated into the model by defining a dimensionless function J(y) which describes how permeability is modulated with depth. The governing equations then be- Following [23] we select a sinusoidally varying function for J(y), representing the fact that a reservoir is stratified with a well-defined wavelength for the strata In what follows, in the interests of maintaining consistency with [23], we set k het = 0.3 and n het = 3. Notice that as a result of this definition J(0) = J(1) = 1. Notice also that the derivative J (y) ≡ dJ/dy is J (y) = −2πk het n het cos(2πn het y).
When evaluated on the approach to y = 1 using the k het and n het values given above, this turns out to be negative with the value −2πk het n het . Suppose we define a function j(y) to be the product of J(y) (the heterogeneity modulation) and y (the net driving pressure), this product being what appears in eqs. (8) and (9). Hence which we have plotted in fig. 2 for our chosen k het and n het values. It follows that Evaluating this at y = 1, we deduce that the value of j at the top (denoted j top ) satisfies j top = 1 − 2πk het n het . This implies that for our chosen k het and n het , not only is J (y) negative at the top boundary, it is actually sufficiently negative that j (y) is also negative there, as fig. 2 indeed indicates. This is what we call a strongly heterogeneous system. The conventional picture whereby points higher up invariably move further and faster than those lower down no longer applies, and the boundary condition given by eq. (6) (or more generally by (7)) ceases to apply. Instead the value of x at the top boundary x top must be obtained as part of the solution of the problem. How the solution can be obtained is discussed in the following sections.  10) and (12)). Parameters describing the heterogeneity are k het = 0.3 and n het = 3. Local maxima (×) and local minima (•) of j are indicated.

Numerical scheme
Equations (8), (9) can be readily solved numerically. Details of the numerical scheme for systems with homogeneous permeability have already been published in the literature [18]. Here therefore we focus mostly on how we adapted that existing scheme to deal with heterogeneous permeability.
Basically the existing numerical scheme is a finite difference scheme in which the front normals at any given point on the front are computed to second order spatial accuracy by examining the locations of neighbouring spatial points either side of the given point but giving more weight to whichever neighbour is closer by. Point positions are then updated in time to second order temporal accuracy using Heun's method.
In the original scheme (for a homogeneous permeability system) the numerical time step was 5 × 10 −5 dimensionless units and the initial spatial separation between material points was 0.025 dimensionless units, data for these choices being found to converge with data generated with much smaller time steps and finer spatial separations [18]. However since the front shape predicted was convex, and since for convex front shapes material points move apart as the front advances, it was necessary to impose a maximum permitted spacing between material points: this was 0.05, i.e. twice the initial spacing. Once the maximum permitted spacing between front material points was exceeded, a new material point was inserted between the existing ones.
Given however that in the present heterogeneous system the permeability is modulated (with n het = 3 wavelengths of the modulation fitting into the vertical extent of the domain) we chose a finer spatial grid here. Specifically we selected an initial point spacing 0.01 and a maximum permitted spacing 0.02. The time step was also chosen to be correspondingly smaller 10 −5 .
As we are now dealing with a front which has concave parts in addition to convex ones, some of the elements actually shrank over time. We therefore also had to assign a minimum permitted point spacing, below which we re-moved points from the front: this minimum spacing was chosen as 0.005. Eventually the concavities shrank down into sharp concave corners. As has been explained in [18] these need to be specially handled in numerical schemes: a corner that turns through an angle δθ needs to be sped up by a factor 1/ cos(δθ/2) relative to the speed of a material point, as otherwise spurious loops develop in the front shapes. This speed up relative to material points is justified, as strictly speaking the concave corners are not themselves material points, but rather geometric points at which material points are being consumed. We considered a corner to be sharp once it turned through an angle that exceeded the square root of the minimum point spacing, ensuring that the curvature at the corner was typically on the order of the inverse square root of the point spacing (and was therefore large). This follows the recommendation of [18] which suggests that the critical corner angle required to apply speed up should be simultaneously small compared to unity but large compared to the point spacing.
Results from computations using the above numerical scheme are presented later on in sect. 5. There are however a number of analytical approximations that we can make which are useful for analysing the front shape: these are considered in sect. 4.

Analytical approximations for front shape
Some approximate analytical solutions are now presented. Section 4.1 considers a small time t 1 approximation. Section 4.2 discusses instead longer time, quasi-static solutions. There turn out to be several local quasi-static solutions on various parts of the front (see sect. 4.3), eventually giving way to a global quasi-static solution (sect. 4.4) over the entire front depth. Section 4.5 meanwhile contains a detailed analysis of the front shape not considering the entire depth, but rather in the neighbourhood of the top boundary.

Small time asymptotic approximation
An asymptotic solution for pressure-driven growth that is valid in the limit of small times, assuming homogeneous permeability, was presented by de Velde Harsenhorst and co-workers [19]: in the literature [25] this has been referred to as the Velde solution. The solution is where the effect of the parameter s 0 (see sect. 2) is treated as negligible. The basis of the Velde solution is that, at early times, material points move primarily in the horizontal direction, and the Velde solution is what one obtains if the vertical motion is neglected in the first instance.
Using the Velde solution it is also possible to estimate the extent of the vertical motion that has been neglected. Via eq. (14), we estimate the derivative of x with respect to y along the front (dx/dy) f ∼ t/(2y) (15) noting that this equals tan α according to eq. (5). Given the way that α has been defined, the vertical component of n is − sin α, but for small time, the value of α is also small and tan α ∼ sin α ∼ α. Substituting eqs. (14) and (15) into (3), we find that front material points drift vertically according to This equation, for the case of a homogeneous permeability system, gives a uniform vertical velocity regardless of y (or more correctly, as will be explained below, a vertical velocity that is independent of y for the overwhelming majority of the y domain). It is possible to employ the vertical motion in eq. (16) to improve the estimate of the horizontal motion: this then leads to an improved Velde solution. The technique for generating this improved solution is explained in [25] (see also an analogous technique in [24]), but here we merely note that in view of eq. (16), a point currently at location y at time t, historically was higher up in the domain, and historically was moving faster than it is currently moving, meaning it manages to reach a slightly higher x, namely In [25] it is shown that this formula, albeit formally applicable for t 1, manages to fit numerical data even for times as large as t = 0.5. Note however that the formula only applies for points which have spent the entire period up to time t drifting downwards with vertical velocity component given by (16): it does not therefore apply for points that find themselves either very close to the bottom boundary (since all motion is arrested there) or very close to the top (specifically the formula does not apply when y > 1 − t/2, and indeed (17) is not compatible with the boundary condition (6)). In the small t limit however, the domains in which eq. (17) breaks down account for just a very small part of the overall y domain.
The above analysis applies to a homogeneous permeability system. In a heterogeneous system, we need to replace y by a function j(y) given by eq. (12). The Velde solution becomes x ∼ s ∼ 2 j(y) t (18) and the vertical drift associated with it becomes The improved Velde solution (after working through the heterogeneous analogue of the technique presented in [24,25]) becomes If we are particularly interested in using eq. (20) to estimate the (a priori unknown) trajectory of x on the top boundary, we could estimate where we have used the fact that in our system j = 1 at the top (by construction), and where also the value of j at the top (j top ) is known: see sect. 2. Despite the similarities between the homogeneous and heterogeneous analyses, there are some significant differences. In the heterogeneous case, for our chosen j (see eqs. (10)- (13) and see also fig. 2), the function j (y) changes in sign according to the y value. Unlike in the homogeneous case in which points on the front had a uniform downward drift velocity, in the heterogeneous case, the vertical drift velocity given by (19) here is non-uniform and can in fact change sign according to the value of y.
Since the rate at which front points are moving horizontally is governed by j (see eq. (18)), it follows that eq. (19) which drives vertical drift opposite to the direction of j is driving front material points from regions of higher speed (high j) to regions of lower speed (low j). This is why eq. (20) always predicts an x value larger than (18) does: regardless of whether points are migrating upwards or downwards, historically they experienced faster horizontal motion than their current y location would indicate.
The assumption used to derive (20) is that the current vertical drift velocity (for a point currently at y) is representative of the drift velocity that a material point has seen throughout its entire evolution up to the current time t. Clearly in a homogeneous system, with uniform vertical drift, that is a good assumption. However, in a strongly heterogeneous system it is a poor assumption, since the rate of vertical motion is y dependent. This suggests that the time domain for which the heterogeneous equation (20) is applicable might be rather less than that for which the homogeneous analogue (17) applies, and consequently (21) might be a comparatively poor approximation to x top in a heterogeneous system.
It is possible moreover to argue that eq. (20) (from which eq. (21) is itself derived) will break down far sooner towards the top of the domain than towards the bottom. The reason is because of the differing behaviour of j and j near the bottom and near the top, as we now explain. Recall that j and j are given by eqs. (12)-(13), where the values of J and J that appear in these equations are in turn given by eqs. (10)- (11). Both J and J are oscillatory, but for our chosen k het and n het values (k het = 0.3, n het = 3) the value of J is typically an order of magnitude larger than that of J: the function J oscillates between ±2πk het n het whereas J oscillates only between 1±k het . As fig. 2 indeed shows, the largest oscillations in j are therefore seen in the upper part of the domain (where the term in y J (y) within eq. (13) is significant), with less oscillation seen in the lower part of the domain (where y J (y) must be smaller). It follows also that in the upper part of the domain j can be up to an order of magnitude larger than J and hence up to an order of magnitude larger than j. Applying this result to the early-time formula for x, whereas in the homogeneous formula (17) the term in t 2 /6 represents, even for times as large as t = 0.5, just a small perturbation to the Velde solution, in the strongly heterogeneous case the analogous term (j ) 2 t 2 /6 in (20) represents a very significant perturbation, particularly in the upper part of the front where j tends to be largest. This again supports the view that eqs. (20)-(21) might only be applicable over just a very limited time domain.
Yet another way of drawing the same conclusion comes from recalling that the Velde solution (18) and the "improved" Velde solution (20) obtained from it are based on the idea that motion is primarily horizontal with only weak vertical displacements. At leading order, the vertical drift velocity is −j /2 (eq. (19)) for a heterogeneous system, compared with − 1 2 (eq. (16)) in a homogeneous case. Given the way that j behaves, the vertical drift velocity can be up to an order of magnitude larger in the heterogeneous system considered here than in the analogous homogeneous system, suggesting that the vertical displacements in the heterogeneous system are not necessarily weak. This then explains why the Velde solution should break down sooner in the heterogeneous case.

Local quasi-static behaviour
Notwithstanding the arguments that we have just advanced concerning the Velde solution possibly being unreliable in heterogeneous systems, there are certain y locations at which it can still be applied to estimate x. Consider the locations at which j is a local maximum (we denote these locations by y max with the corresponding j values being j max ). According to the small time asymptotic expansion (see eq. (19)), these same locations should have no vertical drift (since j = 0 there). At (or near) these locations therefore, the location of the front should match well with the Velde solution (even if the rest of the front does not). In particular the x coordinate of the front should be at a location that we denote x max satisfying In the homogeneous system the maximum value of x always occurred at the top of the front y = 1, although in a strongly heterogeneous system the local maxima x max are expected to be at smaller values of y. What [18] demonstrated, at least for a homogeneous system, was that the shape of the front behind the maximum x value developed over time a quasi-static shape, in which the "apparent" horizontal motion of the front was uniform with height. The actual motion of front elements in this case was normal to the front, but since the front normal was oriented obliquely to the horizontal, the "apparent" motion (in the sense indicated in fig. 3) differed from the actual motion. Returning to consider the heterogeneous case, it is valid to ask by analogy, whether the front shape in the neighbourhood of x max ever attains a quasi-static shape, and whether the resulting "apparent" horizontal motion is likewise uniform with height.
The equations governing the (heterogeneous) quasistatic shape can be derived as follows. Consider (see fig. 3) a point at a y value slightly below y max , such that the normal to the front at y is at an angle α from the horizontal (with α ≡ 0 at y max itself, but not at other y values). In a small interval of time δt, the point displaces through a distance j(y) δt/x max , assuming that up to time Fig. 3. At a given time t, a local maximum of the front is at location (x max, ymax), and in a small time δt, this maximum displaces horizontally by j max δt/xmax. Although a material point at a nearby but otherwise arbitrary vertical location (denoted y), actually displaces normally to the front by an amount j(y) δt/x max (moving at an angle α to the horizontal), a new material point arrives at this location y starting from both higher up and further ahead horizontally (motion of this point indicated here by a dashed arrow). The net result, if the front shape is to be quasi-static, is that the entire section of front apparently displaces in the horizontal uniformly, with the same velocity j max/xmax as the local maximum.
t all points in the neighbourhood of (x max , y max ) have displaced by a comparable path length namely x max . If the apparent horizontal motion of the quasi-static shape is δx max = j max δt/x max uniformly for all y in the neighbourhood, then we require that the projection of δx max from the horizontal onto the front normal at location y matches the normal displacement j(y) δt/x max , and hence This is entirely analogous to the quasi-static solution for the homogeneous case as presented in [18], which actually satisfies cos α = y (even though that is not obvious from the way in which the solution is written in [18]; noting however that cot α = y/ 1 − y 2 , equivalence with the results of [18] can be obtained).
Returning to consider the heterogeneous case, and using eq. (23), we can derive quasi-static formulae for x vs. y as follows: (25) Using the j values plotted in fig. 2, the predicted quasistatic front shapes x vs. y via eq. (25) are plotted in fig. 4. The function j vs. y in fig. 2 has three local maxima (which we can denote from top to bottom as j max,1 , j max,2 , j max,3 at respective locations y max,1 , y max,2 , y max,3 ). Hence there are likewise three branches of local quasi-static solutions in fig. 4. To simplify the plot in fig. 4 we have introduced an arbitrary horizontal shift into the solutions corresponding to the two lower local maxima, to ensure that all the solutions join up continuously at concave corners, the corners being placed at the y locations corresponding to the local minima of j.

Comparison between local quasi-static solutions
Although it is possible to construct local quasi-static solutions x vs. y about each local maximum, y max,1 , y max,2 , y max,3 , it is not possible in general to match together these local solutions in a uniform fashion to construct a global quasi-static solution. The reason why such matching is impossible is that each local maximum has a different x max value: we denote the values by x max,1 , x max,2 , x max,3 etc. and moreover (since they correspond to different j max values) they move over time at different rates (as eq. (22) predicts). Thus, even though it might be possible to match both x and y values for adjacent local solutions at one particular instant in time (see, e.g., fig. 4 where the solutions have been matched via an arbitrary horizontal shift), it will not be possible to retain such matching at later times. Instead the fastest moving of the local maxima (which typically is the one that is higher up than any of the others, since the function j tends to be larger in the upper part of the domain than lower down) will therefore run ahead of the others. As it runs ahead, this fastest moving maximum will start to invade the domain occupied by the slower moving local maxima underneath it. A global quasi-static solution must therefore be based on the fastest moving local maximum and not on any of the others: the implications of this are explored below.

Global quasi-static solution
A global quasi-static solution for the front shape is now described. Suppose the global maximum of j is j max,1 oc-curring at y = y max,1 . Equation (25) is now assumed to hold over the entire solution domain (instead of just a local region centred about y max,1 ): the relevant solution for the function j obtained via eqs. (10) and (12) is plotted in fig. 4. Extending the domain in this fashion however predicts an interesting feature of the global quasi-static solution, described as follows. As y decreases below y max,1 , the value of j(y) decreases, cos α decreases according to eq. (23), and hence α increases: the front reorients as we expect. Eventually however a local minimum of j(y) is reached (see, e.g., fig. 2). Moving downwards below this, according to eq. (23), α starts to decrease again. This then signals a change in the front shape in fig. 4 from convex to concave.
It is known, however (see sect. 1 and also [18]) that concave shapes are problematic in the pressure-driven growth model, since they have a tendency to focus down into sharp concave corners. It is valid to ask therefore how it is possible for a concavity to be sustained indefinitely in a quasi-static solution. The reason turns out to be that individual elements of front only spend a finite amount of time in the concave region: they migrate into the concavity at one vertical location and then migrate out of it again at another vertical location before they have shrunk down to a sharp corner. The proof of this result is given below.
We use the symbol S to denote distance measured along a foam front downward from the top: this is not to be confused with s which is the distance that front material points displace during their trajectory. Consider an element on the front of size δS. Following a foam front element, the rate at which δS evolves, denotedδ S, is known to satisfy (see [27] where κ is front curvature, and v ⊥ denotes the normal velocity of the front (which, in the case of interest here, can be obtained from (9)). Given that we are following a foam front element, and since the system is assumed quasi-static, we can replace the time derivative here by a convected derivative, where z is a coordinate measured down from the top (as opposed to y which is measured up from the bottom), and where v z is the velocity component in the z direction. Based on the definition of the orientation angle α, we know that v z = v ⊥ sin α (28) and also (based on the definition of curvature) where recall S is a coordinate along the front measured downwards from the top, and where we have used the fact that dz/dS = cos α. Substituting eqs. (28)-(29) into (27) it follows that sin α d(δS)/dz = δS cos α dα/dz (30) and it is easy to check that the solution of this is Consider a front element that starts off at some y location slightly underneath y max,1 . As this front element migrates downward, its α value will grow and hence according to eq. (31) its δS value will likewise grow (the element will stretch). Eventually however the element will move from a convex to a concave region of the front: the value of α now starts to decrease, so the element's δS value must then begin to decrease also (the element shrinks). However δS never shrinks down to nothing, since sin α never falls all the way back to zero. The smallest value of α or equivalently the largest value of cos α that we attain, can be readily computed. Via eq. (23), the cos α in question corresponds to a local maximum of j/j max,1 : this is always strictly less than unity, since none of the other local maxima for j (j max,2 , j max,3 , etc.) are ever as large as j max,1 . Thus even though front elements are indeed shrinking when they are on concave regions of the front, they manage to migrate reasonably quickly onto the next convex region of the front where they manage to start growing again.
To summarise, even though the front shape admits sharp concave corners at finite times, the long-time, quasistatic front shape admits smooth concavities (but not sharp concave corners). As we will see later on, whilst sharp concave corners are in fact transiently produced during the evolution of this heterogeneous system, those sharp corners are driven to the bottom boundary of the system where motion of the front ceases, meaning these corners play no role in the long-time, quasi-static behaviour.

Front orientation at the top boundary
The analysis that we have presented to date has been focussed mostly upon what happens below the global maximum, i.e., below y max,1 . Here however we switch focus to what happens above y max,1 and in particular what happens adjacent to the top boundary at y = 1. In particular we want to know, at the top of the domain, what angle (denoted α top ) the front normal makes to the top boundary. We know that initially α top ≡ 0, since the front starts off as a vertical line. Over time however α top should evolve, and we wish to establish whether α top will ever attain a long-time, steady-state value, denoted α top,ss say. In fact, according to eq. (23), the long-time, quasi-static value of α top should satisfy cos α top,ss = j top /j max,1 , where j top is the value of j at the top surface, with in fact j top ≡ 1 in our system (see eqs. (10) and (12)). It follows therefore that α top,ss = − arccos(1/j max,1 ) the sign convention being that α is negative when the front normal has a component pointing upwards.  Fig. 5. Sketch of a foam front shape immediately adjacent to the location x top at the top boundary of a reservoir, the shape consisting of a sharply curved "inner solution" that is realised over a dimensionless distance on the order of γ (with γ 1) matched onto a much lower curvature "outer solution". The inner solution manages to meet the top boundary at a right angle, whereas the outer solution is oriented obliquely, with the normal n to the front being at an angle |α top| above the horizontal.
A non-zero α top is counterintuitive since, as has already been mentioned in sect. 2, conventionally in pressure-driven growth [18,23] (either in homogeneous systems or in weakly heterogeneous systems which still have j > 0 at y = 1) we impose a boundary condition α top = 0. It is only when j top < 0 that the pressure-driven growth system no longer satisfies the boundary condition α top = 0. This leads to a paradox since the reason for imposing this condition is to prevent the foam front from penetrating the top boundary, and this "no penetration" condition now appears to be violated.
The way to resolve the paradox is to recall from [18] that the pressure-driven growth model is in fact a singular limit of a more general model known as the viscous froth, which has been extensively studied in the literature [27][28][29][30][31][32][33][34][35]. It is found that there needs to be a "boundary layer" immediately adjacent to the top surface within which all terms in the viscous froth model must be retained, even though certain terms would ordinarily be neglected in the pressure-driven growth limit. As indicated in fig. 5, eq. (25) represents an "outer solution", but there is also an "inner solution" or "boundary layer" across which the front reorients abruptly to avoid penetrating the top boundary. As fig. 5 shows, the front shape in the boundary layer is necessarily concave (since the front orientation angle α must change from zero right at the top to a value α top given, e.g., by eq. (33) moving downwards through the boundary layer). There is however no risk of this concavity focussing down to a sharp corner (because a sharp corner can only arise in the context of pressuredriven growth, but not for a viscous froth [18]).
As per the discussion in [18], the viscous model here (assumed to be in dimensionless form) satisfies where x top is the location of the front on the top boundary, Δp is the (dimensionless) net driving pressure difference, γ is a (weak) dimensionless surface tension term and K is the front curvature (but defined now with a sign convection such that K is positive in concave regions, i.e. K = −κ).
Since we are interested primarily in a long-time, quasistatic state, we suppose that x top remains a fixed distance behind x max,1 . Since however x max,1 is growing with time, it follows that x top in eq. (34) can be well approximated (at least in relative terms) by x max,1 , and hence Moreover, near the top boundary, the dimensionless net driving pressure becomes Δp ≡ 1. Finally it is recalled that γ is a small dimensionless parameter. As is explained in [18], in the current context, γ should not be thought of as "surface tension" in the conventional sense. Rather it is a parameter that limits the growth of the front curvature: although represented here by a 1-D curve, in reality the foam front has a finite thickness and it cannot curve over length scales less than its thickness. To generate the "inner solution" at the top boundary suppose we look for a quasi-static front shape of the form where ξ is a function we must determine. Differentiating with respect to time gives the apparent horizontal speed of the front as j max,1 /(2t), and it then follows: Substituting eqs. (35) and (37) into (34), we find j max,1 cos α = 1 + γK, an equation which admits a quasi-static solution (in which α and K depend on y but not on t). Moreover where S is a distance measured downward along the front, where we recall our sign convention is such that α becomes increasingly negative as we move downward from the top boundary, and where we have used the fact that dy/dS = − cos α. Thus j max,1 cos α = 1 + γ cos α dα/dy.
Equations of this form have been studied by [36] and implicit solutions (giving y as a function of α) are easily obtained. Moreover dξ/dy = tan α, so once y vs. α is known, the function y vs. ξ can also be obtained. We will not give the detailed formulae for these "inner solutions" here as they are already published in literature [36]. The key features of these solutions are however as follows. At the top boundary y = 1, the "inner solution" has α = 0 and hence has a large curvature K = (j max,1 − 1)/γ, remembering that γ is a small parameter. If we move a vertical distance on the order of γ below the top surface, the front reorients and the curvature decays. At a distance such that the curvature term has fallen to a negligible value, α attains the value α top,ss as given by (33), corresponding to the aforementioned "outer solution".
It is already known in the context of foam improved oil recovery within homogeneous systems that "outer" and "inner" solutions are actually required to represent the front shape [18]. However the outer solution is rather different in the homogeneous case, as it already manages to satisfy α top = 0. As a result, the inner solution is comparatively unimportant, serving not so much to reorient the foam front (as must happen for a strongly heterogeneous case), but rather just to relax a weak singularity in the front curvature. Moreover, in the homogeneous system, the front shapes predicted by both "inner" and "outer" solutions are convex [18]. Here, for a strongly heterogeneous system, we have a quite distinct inner solution which is concave, and which matches on to an outer solution with non-zero α top . Whereas in [18] individual material elements of the front could migrate from an "inner" region to an "outer" one, here with strong heterogeneity the migration is in the opposite sense (from the "outer" to the "inner" region), and once in that "inner" region, elements reorient their motion, so as ultimately to travel parallel to the boundary. The paradox of the pressuredriven growth model in strongly heterogeneous systems appearing to violate the no penetration condition at the top boundary is thereby resolved.
In what follows we will focus exclusively upon the "outer solution" namely pressure-driven growth. Important results in the case of pressure-driven growth are (33) for α top and (35) for x top both in the long-time, quasistatic limit. According to eqs. (9) and (12), the actual speed of material points near the top boundary is less than that of material points at the location y = y max,1 because j at the top is smaller than j max,1 . Nonetheless the fact that the front meets the top surface obliquely means that the apparent horizontal motion along the top boundary exceeds the actual speed of material points at the top. The value of x top then grows like 2j max,1 t (as per eq. (35)) which exceeds √ 2t (see eq. (6), being the value that x top would have had if the front had met the top boundary at right angles instead of obliquely).
This completes the discussion of the various analytical approximations for the front shape. In the next section we will ascertain how the analytical approximations compare and contrast with numerical data.

Results
In what follows we present numerical data obtained via the method outlined in sect. 3, comparing these data with analytical predictions from sect. 4 where relevant. Section 5.1 examines data for front shapes (x vs. y for various t), whereas sect. 5.2 focusses on the situation at the top boundary. Finally, sect. 5.3 examines data for the path length s through which material points on the foam front displace.

Front shapes
Numerically computed front shapes at a selection of times from between t = 0.125 and t = 8 are plotted in fig. 6. It is clear that the front evolves from a configuration that is near vertical at early times to one that is much more tilted over at later times. The front has a tendency to form concave corners: these are very evident for t = 0.125 through t = 1, but are rather less evident for t = 2 through t = 8. The concave corners do not remain at a fixed y location but instead are seen to migrate downwards. Indeed if one looks closely at the t = 2 and t = 4 curves, very small concave corners can be detected towards the bottom of the domain. The more striking feature of the t = 2 through t = 8 curves, however, is that they admit (smooth) concavities which are clearly distinct from concave corners: this is consistent with the predictions of sect. 4.4.
In fig. 7 some early-time numerical data are compared with the predictions of the Velde solution (18) and the improved Velde solution (20). In fig. 7(a) for t = 0.0625 we see very clearly that the improved Velde solution is a better fit. In fig. 7(b) for t = 0.125 however, the "improvement" is less clear. On average the improved Velde performs better, since the Velde solution is consistently behind the numerical curve, yet the improved Velde solution is sometimes behind and sometimes ahead. However pointwise at t = 0.125 the improved Velde solution is no better, in the sense that there are points (particularly those in the upper part of the domain) at which the "improved" Velde solution is a similar distance away from the numerical data as the Velde solution is. As anticipated in sect. 4.1, even though the improved Velde solution is known to work well for homogeneous systems even out to t = 0.5, for a heterogeneous system it breaks down much sooner, and that breakdown manifests itself primarily in the upper part of the solution domain.
A far better fit to the t = 0.125 data can in fact be obtained via the quasi-static solutions given via eq. (25). As is seen in fig. 8(a), locally a quasi-static solution seems to develop centred about each of the local maxima in x, growing outwards from there. The agreement with the nu- merical data is good especially given that there are no adjustable fitting parameters here, as each local maximum in the x direction is given by eq. (22) with the rest of each curve then being determined from eq. (25). There are however mismatches where the quasi-static solutions associated with different local maxima fail to join up continuously: in the numerical data, these manifest themselves in the form of concave corners that run ahead of the quasistatic predictions.
The situation at t = 0.5 has evolved somewhat, as fig. 8(b) shows. There are still recognizable local quasistatic solutions about each of the local maxima in x, but there is a clear tendency of the quasi-static solutions higher up to advance further and thereby invade parts of the domain formerly occupied by the quasi-static solutions associated with maxima lower down. This is associated with the concave corners at the junctions between the various local quasi-static solutions having migrated downwards relative to their earlier locations in fig. 8(a).
By t = 4 (see fig. 9) there is a recognisable tendency for the foam front to converge to a global quasi-static solution. However it is clear that convergence is not uniform: agreement with the global quasi-static solution is achieved quickly in the upper part of the domain, with convergence in the lower part of the domain being rather slower. One of the issues faced here is that for the numerical data y necessarily falls to zero when x = 0 but for the long-time, quasi-static state (see eq. (25)), y only falls to zero as It is only as t increases from t = 4 to t = 8 to t = 16 that we see the numerical curve in the lower part of the domain starting to migrate towards the long-time quasi-static state curve.
Despite the slow convergence, the fact that the solution does indeed evolve towards a long-time, quasi-static state is a significant one from the point of view of improved oil recovery operations. By construction (see eq. (25)), over the overwhelming majority of the depth of the domain (with the exception of y 1), the location x of the front is no more than an O(1) distance behind x max , i.e. typically a point on the front is no further behind the leading x value than the domain itself is deep. This implies that the advancing foam front manages to avoid so called "gravity override" an undesirable situation in which isolated parts of the front, typically those higher up in the domain, far outrun the majority of the front.

Front motion along the top boundary
Another observation that we can make from the curve shapes in, e.g., fig. 6 through fig. 9, is that the curves do not quite meet the top boundary at right angles (i.e., α top = 0). As was discussed in sect. 4.5, this is associated with x top growing slightly faster than it would for a right angle configuration. Here we study in detail how x top vs. t behaves.
In fig. 10(a) the numerical x top vs. t data are compared against two formulae: √ 2t (eq. (6), the formula that would apply if α top = 0 at all times) and 2j max,1 t (eq. (35), the formula that would apply if α top = α top,ss ≡ − arccos(1/j max,1 ) at all times). By time t = 1 it is clear that the numerical data fit the latter formula far better than they fit the former. However it is difficult to detect which formula fits better for small times t 1, since all the curves converge together on the plot as t → 0. In order to visualise this better, in fig. 10 Over time this is seen to migrate from 1 to close to j max,1 (corresponding to the angle α top reorienting from zero at t = 0 to a final steady-state value α top,ss ≡ − arccos(1/j max,1 )).
The evolution of α top with time is shown explicitly in fig. 10(c). Clearly there is some noise in the numerical data: this is actually a numerical artifact arising from the fact that we quite often have to delete front material points near the top boundary. As explained in sect. 3, we delete material points whenever the point spacing falls below a minimum permitted value. Near the top boundary, material points are predicted to have an upwards velocity component, so the spacing between these material points and the top boundary is continually decreasing. In spite of the noise, it is clear that the data match the predicted steady state α top,ss .
It is possible to use the Velde and improved Velde formulae (eqs. (18) and (20)) to determine (via eq. (5)) a value of α top : the resulting formulae are also plotted on fig. 10(c). Clearly the Velde and improved Velde formulae for α top only apply at very early times: for any time t greater than about 0.05 the numerical data tend to be closer to α top,ss than to the Velde data.
The data we have considered here have a well defined value of j top (with j top being negative as is clear from fig. 2): based on eqs. (10)- (13), j top evaluates to 1 − 2πk het n het as already mentioned in sect. 2. It is interesting to speculate what might happen in a situation where the magnitude of j top is decreased. This could be achieved for instance by decreasing the value of the parameter k het , which would maintain the value of j at the top, denoted j top , equal to unity. Reducing k het from the original choice k het = 0.3 down to k het = 1/(2πn het ), would make j top tend to zero. Since vanishing j corresponds to having a maximum in j, the more that the magnitude of j top is decreased, the closer j max,1 must move to j top ≡ 1. The value of α top,ss given by eq. (33) must then migrate towards zero, implying that the x top vs. t formula from eq. (6) will be recovered.
The discussion presented above concerned the case j top ≤ 0, but with the magnitude of j top decreasing towards zero. A boundary condition α top = 0 is also typically imposed on systems for which j top is strictly greater than zero, e.g. in the case of homogeneous permeability we have j ≡ y and hence j top = 1. It is known however that imposing α top = 0 when j top is strictly greater than zero leads to considerable complexity in the structure of the foam front at early times and near the top [25]: this arises due to an incompatibility between those material points forced to remain on the top boundary and points immediately below the top which must drift downwards with velocity component −j top /2 according to the Velde solution (see eq. (19)). This incompatibility is removed in the special case when j top → 0, suggesting that the complicated structure near the top of the foam front needed to resolve that incompatibility (as presented by [25] for a homogeneous system) is not actually required in the special case of a heterogeneous system with vanishing j top .
In what follows however we return to the case of nonvanishing j top and more particularly strongly heterogeneous systems with j top < 0.

Behaviour of path length variable s
So far we have focussed attention primarily upon front shapes and quantities (e.g., x top , α top ) which can be obtained directly from them. Here we switch the focus instead to the path length variable s. Although this must be computed in order to determine the front shape, the evolution eq. (8) for the front shape being dependent upon it, we have not yet examined its behaviour in detail.
We know that s is necessarily greater than x, because to reach a given x location, front material points generally must displace both horizontally and vertically. However since front motion is initially predominately horizontal, s should be only slightly greater than x (at least initially). Figure 11(a) shows a plot of s vs. y superposed on a plot of x vs. y for various times up to t = 0.5. As expected, the s values are just very slightly larger than x.
The biggest discrepancies between s and x tend to be seen in the regions immediately above concave corners in the front shape. This is unsurprising. We know (via fig. 8) that in the time regime considered here, the front shape is well described by local quasi-static solutions satisfying eq. (23). For s to grow above x we need a significant vertical component of motion, i.e. we need α different from zero, and this is achieved via eq. (23) when j is significantly smaller than j max . Note that, the angles α are not expected to be symmetric either side of a concave corner: selecting the y value corresponding to the corner between the two uppermost local quasi-static solutions, j(y)/j max,1 will be less than j(y)/j max,2 for instance, these values then determining α according to eq. (23). The front is therefore tilted over to a greater extent above the corner than below, which also contributes to the gradual downward migration of the corner itself, since the algorithm proposed by [18] sets the direction of motion of the corner to be along the corner bisector. Another consequence of having a larger α above the corner than below is that s is larger above the corner (because as explained above, larger α implies more vertical motion, which in turn implies a larger path length s travelled to reach a given horizontal location x).  It is clear therefore that there is actually a jump in s at each of the concave corners which is what fig. 11(a) shows. This is the first time that a jump in the value of path length s has been clearly demonstrated in numerical results of the pressure-driven growth model. The result is actually a fairly significant one for the numerical implementation of the model, since the rule we use for evolving the corners i.e. a speed up factor relative to material points of 1/ cos(δθ/2) where δθ is the angle through which the corner turns, strictly speaking only applies when material points either side of the corner are moving at the same speed (but in different directions). In reality, the speed of material points is sensitive to s (see eq. (9)), so if there is a jump in s at the corner, there must also be a jump in that material point speed.
Strictly speaking then we ought to be using a more general technique for determining how the concave corners move: indeed a general technique for determining how sharp corners must evolve in pressure-driven growth when the speed of material points undergoes a jump at a concave corner has already been derived by [22]. The derivation specifically envisaged the case of an anisotropic permeability, so that a jump in orientation of the front (as necessarily happens at a concave corner) implies a jump in speed of material points, regardless of whether or not path length s undergoes a jump at the corner. In the present system we have a heterogeneous but isotropic permeability, but the technique derived by [22] for propagating the concave corners can still be employed, as all that matters in the technique is that there is a jump in the speed of material points regardless of how it is caused.
We leave the implementation of this more general technique for computing the evolution of the locations of the concave corners for future work. We note however from fig. 11(a) that the jumps in s compared to the value of s itself are fairly small, so the error involved in propagating the concave corners assuming equal s and hence (according to eq. (9)) equal speed either side of the corner should likewise be small. Moreover the mathematical structure of the pressure-driven growth equations should help to maintain the integrity of the solution in spite of the exact details of how the corner is propagated: information about the front shape tends to propagate outwards in convex regions but inwards in concave regions (towards the concave corner). As front material points migrate closer and closer together at the concavity and thereby the spacing between them diminishes, we eventually need to remove points from the front altogether: in the neighbourhood of the corner then we destroy information, meaning that slight inaccuracies in that information cease to be relevant. Sufficiently large inaccuracies can and do produce spurious behaviour (see [18]), but slight inaccuracies can be tolerated (since the inaccurate information itself is eventually destroyed).
Yet another consideration is that as time proceeds the concave corners themselves become less prominent so the algorithm for propagating them should likewise become less critical: fig. 11(b) for instance shows x and s vs. y at longer times up to t = 4. The concave corners on the x vs. y plots are now difficult to see, although jumps in s at certain y values remain evident: these jumps in s are however gradually being pushed to the bottom of the solution domain, beyond which no further advance of the foam front is possible.

Conclusions
We have used the pressure-driven growth model in order to study the advance of a foam front through an oil reservoir in the context of foam improved oil recovery. Whereas previous studies have focussed on homogeneous and/or weakly heterogeneous systems, here we have considered a strongly heterogeneous case: strong heterogeneity in this context means that the reservoir permeability varies significantly from stratum to stratum, and in particular decreases sufficiently in the direction moving toward the top boundary, that it overcomes the tendency of the net driving pressure (injection pressure less hydrostatic pressure) to grow moving upwards. Thus the part of the foam front that moves the furthest and the fastest is not that right at the top of the reservoir, but rather somewhere in the interior, albeit typically relatively close to the top.
An early-time solution for the shape of the front is available that assumes that the motion of front material points is predominately in the horizontal. This solution predicts that the horizontal displacement of the foam front develops several local maxima, each maximum corresponding to a high permeability stratum. This solution however breaks down comparatively quickly (much more quickly than the analogous early-time solution for a homogeneous system does) since it turns out that vertical motion of front material points can be rather more significant in a strongly heterogeneous system than in a homogeneous one.
A better match to the front shape can be found by assuming a local quasi-static solution about each of the local maxima of the front displacement. Each local quasistatic solution assumes a locally uniform "apparent" horizontal motion matching that of the corresponding local maximum. The various local solutions cannot however be matched together to give a global quasi-static solution since it turns out that the different local maxima move at different rates from one another. Over time then, the faster moving local quasi-static solutions tend to invade the domains occupied by their slower moving neighbours with the result that the global long-time, quasi-static solution is that corresponding to the fastest moving out of all the local maxima: this is typically the maximum that is located highest up in the domain. All points on the front are now predicted to have an apparent horizontal motion equal to that of the leading local maximum.
One interesting prediction is that the global, long-time solution admits concave regions in the foam front shape. Unlike what is normally seen in pressure-driven growth, these are smooth concavities rather than sharp concave corners: individual material elements spend only a limited time migrating through any given concavity, implying insufficient time for the concavity to sharpen. Sharp concave corners do in fact appear in the front shape early on, but they turn out to be transient being driven to the bottom boundary of the front where they can advance no further: the concavities that persist at long times are smooth ones.
The motion of the uppermost of the local maxima on the front (associated with the uppermost of the high permeability strata) also influences the solution for the front at the top boundary of the domain. The result is that the front along the top boundary moves faster in a strongly heterogeneous system than in a homogeneous or weakly heterogeneous one. The foam front also appears to meet the top boundary obliquely, rather than (as is usual for pressure-driven growth) at a right angle: this reflects a rapid spatial re-orientation of the foam front near the top boundary, which can be captured by a more general model (the viscous froth, an "inner solution" valid very near the boundary) albeit not via pressure-driven growth (an "outer solution").
The prediction of a long-time, global quasi-static solution for pressure-driven growth even in a strongly heterogeneous system such as this one is significant for petroleum engineering operations. Since the global quasi-static solution predicts that all points move with the same apparent horizontal motion as the forwardmost point (i.e., as the leading local maximum), this means that the front itself remains just a limited distance behind that forwardmost point. This then suggests that override (an undesirable situation in which certain points on the front can advance arbitrarily far ahead of others) should not occur in this system.

Open Access
This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.