Foam front propagation in anisotropic oil reservoirs

The pressure-driven growth model is considered, describing the motion of a foam front through an oil reservoir during foam improved oil recovery, foam being formed as gas advances into an initially liquid-filled reservoir. In the model, the foam front is represented by a set of so-called “material points” that track the advance of gas into the liquid-filled region. According to the model, the shape of the foam front is prone to develop concave sharply curved concavities, where the orientation of the front changes rapidly over a small spatial distance: these are referred to as “concave corners”. These concave corners need to be propagated differently from the material points on the foam front itself. Typically the corner must move faster than those material points, otherwise spurious numerical artifacts develop in the computed shape of the front. A propagation rule or “speed up” rule is derived for the concave corners, which is shown to be sensitive to the level of anisotropy in the permeability of the reservoir and also sensitive to the orientation of the corners themselves. In particular if a corner in an anisotropic reservoir were to be propagated according to an isotropic speed up rule, this might not be sufficient to suppress spurious numerical artifacts, at least for certain orientations of the corner. On the other hand, systems that are both heterogeneous and anisotropic tend to be well behaved numerically, regardless of whether one uses the isotropic or anisotropic speed up rule for corners. This comes about because, in the heterogeneous and anisotropic case, the orientation of the corner is such that the “correct” anisotropic speed is just very slightly less than the “incorrect” isotropic one. The anisotropic rule does however manage to keep the corner very slightly sharper than the isotropic rule does.


Introduction
During petroleum production, only a fraction of the oil present in an oil reservoir is extracted via primary production (i.e. driven out of the reservoir under the reservoir's own pressure). Subsequently [1] fluids must be injected into the reservoir with the aim of raising pressure and displacing the oil that is still present. Foam offers advantages [2,3] over other possible choices of driving fluids (e.g. water or gas) for various reasons. Foams have comparatively low mobility [4,5] within porous media such as oil reservoirs. This low mobility arises [6] not only from some bubbles in the foam becoming trapped within pores, but also from high levels of viscous dissipation associated with those bubbles that are moving. Thus, whereas gas in the absence of foam is highly mobile, gas that is within foam has much lower mobility, and this suppresses fingering or channelling phenomena (which otherwise drive flow a e-mail: paul.grassia@strath.ac.uk along preferential flow paths bypassing much of the oil). As a result [6], foams give a comparatively uniform sweep through a reservoir reaching parts of the reservoir that other displacing fluids might not reach.
Injection of low mobility fluid into an oil reservoir would ordinarily demand very high driving pressures. Producing foam via co-injection of gas and surfactant solution therefore requires high pressures. This situation can be mitigated by injecting first surfactant solution and then gas (the so-called surfactant-alternating-gas process [7][8][9][10]), and this is the situation that we study here. Foam (i.e. bubbles of gas separated by thin liquid films, typically with the bubbles being longer than the pore is wide) is now formed in situ at the interface of the surfactant solution and gas, and the resulting zone of low mobility is then confined close to this interface, sufficient to ensure uniform sweep without requiring excessive overall driving pressure [6]. The degree to which the low mobility zone is localised at the foam front depends on how readily the foam collapses upstream of the front: less surfac- The evolution of the foam front can be described by discretising it into segments and following the evolution (in an X, Y coordinate system) of material points which are moving to the right and downwards, with the normal to a front segment being oriented at an angle α from the horizontal. Although the front shape is predominantly convex (viewed from downstream) it is possible for concave corners to form at discrete points on the front. (b) In the case where the permeability is anisotropic (typically being larger in the horizontal than in the vertical), the direction in which the front material points are moving is not necessarily the same as the direction of the front normal. As drawn here, t h en o r m a li sa ta na n g l eα from the horizontal while the dir e c t i o no fm o t i o ni sa ta na n g l eβ α from the horizontal, with β α being a function of α, but not equal to α. tant implies faster collapse and hence a more localised low mobility zone [6]. If the surfactant concentration evolves over time (due e.g. to downward migration of surfactant solution under gravity [11] or perhaps due to surfactant adsorbing onto pore walls), then the extent of the low mobility zone should be smaller than in the case of constant surfactant concentration: such complications are not however considered here. The presence of oil also helps to favour foam collapse [12], again tending to imply a low mobility zone of smaller extent.
In order to exploit foams in improved oil recovery processes, is it useful to be able to predict how the foam moves through the oil reservoir. Towards this end, computationally intensive foam simulation techniques have been developed [13][14][15][16][17][18]. A simpler alternative however is to use a "pressure-driven growth" model [6,19], which is a very simple description of the evolving shape of a foam front moving through an oil reservoir. Specifically the model makes it possible to represent the front as a set of material segments and to track the trajectories of material points on these segments (thereby tracking the advance of gas into the liquid-filled region): see fig. 1(a). The directions in which material points move depend on the orientation of the material segment on which they are located, whilst the speeds of material points turn out to depend on their instantaneous location and their path history (details to be given later). By considering a sufficient number of such material points it is thereby possible to reconstruct a discretised representation of the resulting foam front shape.
The physics behind the model is very simple [6,19]: the foam front advances according to the difference between an injection pressure and a hydrostatic pressure in the reservoir. Since the hydrostatic pressure grows with depth, the top of the front advances the furthest, while points lower down advance less [19]. The net driving pressure (injection pressure less hydrostatic) is balanced by dissipative losses (the pressure losses associated with both trapping films and dragging them through the reservoir being lumped together). These dissipative losses are assumed to be concentrated at the foam front itself, in the aforementioned low mobility zone, where the number of foam films is greatest. Upstream of the foam front [6,20] there are relatively few films (as films tend to have already ruptured upstream) and downstream of the front there are no foam films (only surfactant solution). Hence, away from the foam front, mobility tends to be high and dissipative losses are negligible there.
In a typical situation at oil production scale, the dissipative foam front region might be on the order of a hundred metres across [19] (containing a multitude of individual bubbles). The front itself advances much more than this, e.g. as much as several kilometres horizontally [19,21]. This reiterates that the foam front itself is an exceedingly thin region compared to the distance that the front has advanced overall. The front shape can then be idealised as a 1-d curve in either 2-d space or more generally in a 3-d axisymmetric space. We focus here on the 2-d case (which is slightly simpler to formulate than the 3-d axisymmetric one, but which gives remarkably similar results [6,19]). Useful data that can be extracted from the pressure-driven growth model include not only rates of advance of foam fronts (and hence times for foam to migrate from an injection well to a production well), but also the front shapes [6,19]. Based on these front shapes (i.e. 1-d curves) it is possible to determine at any given time, how much of the reservoir has been swept by foam, and how much is left unswept, and more particularly where those unswept regions are.
The front shapes resulting from the model can be classed into two broad categories: wholly convex shapes (viewed from downstream to upstream) and shapes containing concavities. For the most part, the pressure-driven growth model predicts convex front shapes and the model handles these with ease [19]. Nonetheless there are certain situations [11,19,22] where concavities are known to arise in the front shape during at least part of its time evolution.
These include gravity-driven drainage of the surfactant solution used to make the foam [11], a sudden increase in the injection pressure [22], and the case of a heterogeneous stratified reservoir (mentioned by [19] but never studied in detail).
Concavities have an interesting dynamics in the pressure-driven growth model, focussing down from a smooth concave curve at early times to an "arbitrarily sharp" concave corner later on [11,19,22]. "Arbitrarily sharp" in this context means that the orientation of the foam front changes on a length scale comparable with the front thickness, which (as mentioned above) is much smaller than the overall distance that the front propagates. In the pressuredriven growth model, these concave corners where two otherwise convex segments of front meet require special treatment. Specifically if, between two adjacent segments of front, the tangent to the front turns through an angle θ at a concave corner, then the corner must move more quickly than surrounding material points [19] by a factor 1/ cos(θ/2). This reflects the fact that the concave corner is not itself a material point but instead a geometric feature of the front shape: indeed material points are continually being consumed by the concave corner. In the language of hyperbolic partial differential equations [23], the material points propagate along characteristics, whilst the concave corner moves as a "shock" (i.e. a discontinuity in the tangent to the front, albeit with the front location itself being continuous), and the speed of the discontinuity differs from that of the characteristics.
If the concave corner is propagated without including the 1/ cos(θ/2) speed correction factor which was mentioned above, serious consequences can occur in a numerical implementation of the pressure-driven growth model [19]. Specifically, points on the front neighbouring the corner are found to cross over one another, forming a spurious loop that is topologically infeasible. Worse still, as a result of the changed topology, the numerical scheme treats the spurious loop as if it had higher pressure than its surroundings causing it to grow indefinitely [19]. It must be emphasised that the formation and growth of these spurious loops are unphysical artifacts of the numerical scheme, and not part of the pressure-driven growth model per se. Including the 1/ cos(θ/2) correction factor in the speed of the corner removes these artifacts [11,19,22].
It was mentioned earlier that the case of a heterogeneous reservoir was one of the cases of interest for which concavities might arise in a foam front. Oil reservoirs are heterogeneous as a result of stratification, some strata having different permeability than others. Another effect of stratification however is that reservoirs might be anisotropic: permeability is higher along a stratum than between strata. The pressure-driven growth model has been reformulated for the case of an anisotropic but nonetheless homogeneous reservoir [21]. This led to a convex front shape, similar to that for the heterogeneous reservoir, albeit with increased "gravity override", i.e. for a given propagation distance of the leading edge of the foam front, the area left unswept beneath the front is greater.
However, in the particular case of a heterogeneous, anisotropic reservoir (as opposed to a homogeneous, anisotropic one), concavities and hence concave corners might be expected to arise. These concave corners presumably then need to be propagated rather more rapidly than surrounding material points, but the question then arises of precisely how rapidly they should move. The 1/ cos(θ/2) speed up rule for concavities, that was mentioned earlier [19], strictly speaking only applies to the case where permeability is isotropic. Turning to anisotropic systems however, using the incorrect speed up factor (e.g. employing the "isotropic" 1/ cos(θ/2) factor instead of a more general formula for anisotropic systems) can potentially lead to the same problematic numerical issues as omitting the speed up factor altogether would, i.e. it can lead to points on the curve crossing over one another giving spurious loops as topologically infeasible numerical artifacts. Thus a speed up factor for concave corners in the presence of anisotropic permeability needs to be determined: this is the main purpose of the current study.
The rest of this work is laid out as follows. Section 2 gives a general formulation of the pressure-driven growth model. Then sect. 3 derives the "speed up" rule for a concave corner in a pressure-driven growth system, accounting for anisotropy. Next sect. 4 analyses the consequences of the anisotropic "speed up" rule for the motion of a concave corner. After that sect. 5 describes a numerical scheme for pressure-driven growth in heterogeneous and anisotropic systems incorporating the motion of concave corners. Results of the numerical scheme are described in sect. 6 demonstrating that good numerical behaviour is obtained provided a corner "speed up" rule is employed: surprisingly however the system is not exceedingly sensitive to which particular "speed up" rule is used, a result which follows owing to the particular orientation of the corners that develop in heterogeneous and anisotropic systems. Conclusions are offered in sect. 7.
2 Formulation: Pressure-driven growth model As alluded to above, foam improved oil recovery generally proceeds by flooding an oil reservoir with surfactant solution and then injecting gas under a specified driving pressure (the so-called "surfactant alternating gas" process [7][8][9][10]21]). A foam front forms at the boundary between the surfactant solution and the gas, and this advances into the reservoir. The foam films encounter resistance as they advance, and this is reflected in a low mobility region near the foam front itself [6] which is the zone where the number density of foam films is highest.
In what follows, equations governing how the foam front moves under these circumstances are given (sect. 2.1) and then are cast in dimensionless form (sect. 2.2). Modified equations incorporating explicit effects of heterogeneity and anisotropy are introduced in sect. 2.3.

Governing equations for the foam front
Finding the time evolution of a foam front (see e.g. the definition sketch in fig. 1(a)) is equivalent to finding the trajectory of a set of material points that give a discretised representation of the front. Suppose that coordinates X, Y represent a front material point at a given time t, with X being measured to the right, and Y being measured upwards. The pressure-driven growth model, which is described in detail by [6,19], and which, in the first instance, is considered in a homogeneous and isotropic system, predicts the rate of change of X and Y with time t to be where k is a reservoir permeability (homogeneous and isotropic in the first instance), λ r is a relative mobility (the reciprocal of an effective viscosity) applicable specifically to the foam front region (motion of the foam front being highly dissipative [6], but mobility elsewhere being exceedingly high), ∆P is a pressure difference (i.e. the difference between an injection pressure P inj and a background hydrostatic pressure P hyd , the latter growing with depth), S g is the gas fraction in the foam specifically at the front where the foam meets the liquid-filled region downstream (S g can be estimated via fractional flow theory [6]), φ is the porosity of the reservoir (which again we assume to be homogeneous), s is the distance that the point on the foam front has propagated to date (since injection), andτ is the ratio between the front thickness and the distance that the front has propagated (which is a constant much smaller than unity; the fact thatτ can be taken to be constant follows from fractional flow theory [6,15,20], which demonstrates that the thickness of the low mobility zone at the foam front grows proportionally to the distance that the front travels). Finally α is the angle that the normal to the front segment makes to the horizontal: assuming that the foam front at any given instant is represented by a curve X vs. Y , then α can be determined as α = arctan(dX/dY ), where the derivative dX/dY is taken at a particular time t, moving along the front from material point to material point. As mentioned above, the hydrostatic pressure P hyd grows with depth, the gradient of hydrostatic pressure being ∆ρ g (where ∆ρ is the density difference between foam and surfactant solution, and g is acceleration due to gravity). The implication is that, for a given injection pressure P inj , there is a maximum depth d max to which foam can penetrate (equal to P inj /(∆ρ g)): this is the point at which injection pressure is exactly balanced by hydrostatic pressure.
Equations (1)-(2) correspond to Darcy's law taken across the low mobility zone at the foam front. Although it might appear that capillary effects have been discarded entirely from the model, the fact that gas mobility within foam is much lower than gas mobility in the absence of foam actually arises due to capillary effects. Remember that what the model is describing here is not an individual foam film in a channel (for which capillary effects might well be included explicitly in the model [24,25]), but rather the evolution of an entire foam front, itself containing a multitude of bubbles (and with the entire foam front now being represented simply by a low mobility zone).

Governing equations in dimensionless form
We make distances dimensionless with respect to the maximum depth d max , and we also place the origin of the spatial coordinate system at this maximum depth. Pressures are likewise made dimensionless with respect to P inj . We also make times dimensionless 1 with respect to S g φd 2 maxτ /(kλ r P inj ). In dimensionless form then, now treating X, Y , s and t as dimensionless variables, we deduce On the right-hand side of eqs. (3)-(4), the term Y represents the dimensionless analogue of the net pressure difference ∆P , with Y varying from 0 (at the bottom of the front) to 1 (at the top surface). Note also that (since s represents the total distance travelled by a material point) it follows ds/dt = Y/s.
Equations (3)-(4) are solved with suitable initial conditions, typically with the foam front being initially vertical. However the front begins to tilt over with time, as material points higher up (which experience a higher net driving pressure) migrate faster than those lower down. As the front tilts over however, points must also begin to move downwards (the angle α grows as the tilt grows, leading to vertical motion in eq. (4)). Hence the generic motion of a typical material point in fig. 1(a) is to the right and downwards. A boundary condition is however also required at the top surface, namely that the direction of travel must be along that surface. This implies α =0 at the top surface and hence dY/dt = 0 there, which also requires distance s travelled at the top equals Cartesian coordinate X. Equation (3) can then be solved for X at the top surface (Y = 1) as a function of time t, the solution (assuming an initial condition X = 0 when t =0 ) being simply X = √ 2t.

Accounting for heterogeneity and anisotropy
The above formulation applies to a homogeneous and isotropic case. Changes to the model accounting for het- (a) In an anisotropic system, two segments of foam front oriented at angles α ± θ/2 from the horizontal, and moving in directions β ± ψ/2 from the horizontal, join up at a concave corner. The apparent speed v app of the concave corner is expected to be greater than that of a hypothetical material segment oriented at angle α, and similarly is expected to be greater than that of a hypothetical material segment moving in direction β. Moreover the apparent motion of the concave corner is in a direction β app, which is the direction of motion corresponding to a hypothetical material element oriented in a direction α app. In general βapp differs from β,a n dαapp differs from α.(b)Inatimeτ , a material point originally at C reaches position B, moving a distance v aniso,C τ along a line with slope − tan β + ≡ tan(β + ψ/2). Likewise in a time τ ,am a t e r i a l point originally at D reaches position B, moving a distance v aniso,D τ along a line with slope − tan β− ≡ tan(β − ψ/2). The concave corner moves from location A to B in this same time period, covering a distance v app τ along a line with slope − tan β app. The original location A of the concave corner can be found given the slopes of the material segments AC and AD, which are respectively 1/ tan(α +) ≡ 1/ tan(α + θ/2) and 1/ tan(α −) ≡ 1/ tan(α − θ/2).
erogeneity [19] and anisotropy [21] are introduced generalising eqs. (3)-(4) as follows: where we continue to present the governing equations in dimensionless form, and where α is (as before) the orientation of the front normal measured with respect to the horizontal, β α is the direction of motion of the material point (again measured with respect to the horizontal, with β α being a function of α, but as shown in fig. 1(b), not identical to α)a n dk v < 1 is an anisotropy factor, i.e. the ratio of vertical to horizontal permeability. Note that β α can be evaluated (see [21]) as arctan(k v tan α), which is derived by considering the ratio of −dY/dt to dX/dt. The term cos(α − β α ) in the denominator represents [21] the fact that fractional flow theory [6,15,20] gives the thickness of the front measured along the direction of motion β α , but the thickness of the front measured along the normal direction α is less than this. Moreover J(Y )i sa heterogeneity factor, the precise form of which depends on how the reservoir is stratified into layers. To illustrate the model, we choose the following form for J(Y ) where k het < 1 is the relative variation of permeability from layer to layer, and n het (which for simplicity we constrain to be an integer) is the number of high and low permeability layers. For simplicity we assume that the porosity of the reservoir (which enters our definition of dimensionless time) is fixed, despite the variations of permeability. However variations of porosity could readily be absorbed into the function J(Y ) if required. In lieu of eq. (5), the distance s that the foam front has propagated now satisfies (remembering that β α equals arctan(k v tan α) as mentioned above, and also employing a trigonometric identity [26]), We also assume (again for simplicity) that the ratio between foam front thickness (measured along the propagation direction) and distance that the front has propagated is independent of permeability variations. This now completes the formulation of the pressuredriven growth model in the case of a heterogeneous and anisotropic system. The important thing to note is that, owing to heterogeneity, material points in low permeability regions (with J(Y ) < 1) might lag behind those in higher permeability regions (with J(Y ) > 1). This can cause concavities to form in the front shape, and those concavities might subsequently focus down into sharp concave corners (such a corner is shown in the definition sketch fig. 1). These corners are not material points, and so propagate differently from material points. We need to know how to propagate those concavities, and in particular (if the system is anisotropic in addition to being heterogeneous) we are faced with the challenge of propagating corners with anisotropy. This is discussed in the next section.

Configuration of a concave corner
Consider a case ( fig. 2(a)) where a section of front oriented at an angle α + ≡ α + θ/2 joins up with a section of front oriented at an angle α − ≡ α−θ/2, both angles being measured from the horizontal to the front normal. Clearly the front turns through a (concave) angle θ at the junction, whilst the definition of α is a little more general than in sect. 2 (previously it was the orientation of a particular section of front, now it is the orientation of the bisector between adjacent segments). Suppose that the directions of motion of the two adjacent segments of front are oriented at angles β + ≡ β + ψ/2a n dβ − ≡ β − ψ/2 (again measured from the horizontal). Note that β as defined here is the bisector of the two propagation directions, which is not necessarily the same as β α (given, as in sect. 2, by the formula arctan(k v tan α)) although β and β α are usually quite close (as we will see). The model described by [21] makes it possible to determine β and ψ in terms of α and θ. In particular where (recall from sect. 2), k v is the ratio between vertical and horizontal permeability, which for a stratified reservoir is taken to be a constant somewhere between zero and unity. Cases in which θ is a small value (i.e. θ ≪ π)a r e common, because the concave corner once it has begun to form, tends to be propagated in such a way as to prevent it sharpening any further. In that case, we obtain β ≈ arctan(k v tan α) (small θ limit). (12) This implies that β α (as defined in sect. 2) and β tend to coincide provided θ is small. Moreover ψ can be obtained (by differentiating this β vs. α relation, and using a relevant trigonometric identity [26])

Motion of a concave corner
Returning to the case of finite θ, despite knowing β and ψ via eqs. (10)-(11), we do not know a priori either the speed or direction of motion of the concave corner. Its speed may be greater than the speed of adjacent material points, and it need not move along the direction β that bisects the motion of adjacent material points, nor for that matter along the direction β α .
The actual motion of the concave corner can however be deduced with reference to fig. 2(b). That figure shows the concave corner instantaneously at location A, and two nearby material points at locations C and D. The material points C and D follow well defined trajectories (in the directions β + ψ/2a n dβ − ψ/2 respectively) and, after a certain elapsed time τ , they coincide at location B.T h e concave corner must in this same time τ move from A to B.
Finding the locations of points A, B, C and D relative to one another is an exercise in coordinate geometry. We elect to work in Cartesian coordinates here, although note that a mathematically elegant vector formulation is also available (see appendix A). In order to proceed, we recall the shorthand notation introduced at the start of this section, symbols α + and α − being defined such that α ± ≡ α ± θ/2, and symbols β + and β − being defined such that β ± ≡ β ± ψ/2. The slopes of the lines BC and BD in an x, y coordinate system (x measured to the right, and y measured upwards, but with the origin of the coordinate system placed at B, and hence translated with respect to the X, Y coordinate system defined in fig. 1) are respectively − tan β + and − tan β − . The lengths of the lines BC and BD now become τ times the speed of each material point, either point C or point D.
A complication is that, in an anisotropic system such as this one, the material point speed depends on the orientation of material segments. Specifically there is a material point speed that we denote v aniso,+ , corresponding to segments of front upon which the front normal is oriented at angle α + from the horizontal and moving in the direction β + . There is a different material point speed, denoted v aniso,− , corresponding to segments of front upon which the front normal is oriented at angle α − and moving in direction β − . These speeds v aniso,± can be related to the speed of a material point in an analogous isotropic system (denoted v iso ) via a formula (which can be derived from the work of [21], using also some trigonometric identities [26]) where v iso is a known and well-defined function of the material point position and of its path history. Via the theory presented in sect. 2, this is equal in fact to where Y is the coordinate defined in fig. 1, and s is path length executed to date, and J(Y ) represents the heterogeneity (if any), and α ± and β ± are (as above) shorthand for α ± θ/2a n dβ ± ψ/2. Equation (14) can be analysed as follows. In the limit as k v → 1, then β ± → α ± , and hence v aniso,± → v iso . Moreover in the opposite limit as k v → 0, then β ± → 0, and v aniso,± is again the same as v iso . However in a special case where θ and ψ are small (so that α ± ≈ α and β ± ≈ β and moreover v aniso,+ ≈ v aniso,− ), and supposing also 2 β ≈ α/2, the ratio v aniso,± /v iso is less than unity, specifically it is cos(α)/ cos 2 (α/2) which (owing to a trigonometric identity [26]) is 1 − tan 2 (α/2).
According to eq. (14), v aniso,± depends not only on α ± and β ± , but also on v iso (i.e. the speed of material points in an "analogous" isotropic system). Remembering that material points tend to move not only to the right, but also downwards, observe that v iso (which equals Ys −1 J(Y )as we have said) tends to fall with falling Y as a material point moves downwards (since the difference between the injection pressure and the reservoir hydrostatic pressure then falls). Moreover v iso falls as a material point executes longer and longer paths thereby increasing s: this is because the dissipative "low mobility" region at the foam front becomes thicker the further the point travels [6]. Finally in the case of a layered heterogeneous reservoir v iso is sensitive to any permeability variations J(Y )fromlayer to layer in the reservoir.
We focus here on a concave corner and material points very close to it. Although there will be some variation in v iso from material point to material point moving away from the corner in either direction, locally this will be insignificant compared to the significant change in front orientation at the corner itself. In spite of this, it is actually possible to contemplate a discontinuous jump in v iso at the corner, since v iso is sensitive to path lengths and front material points may have different path histories depending upon which side of the corner they are located. In the interests of keeping the presentation simple however, these possible discontinuities in v iso at the corner will be neglected (numerical evidence suggests they are less significant than the changes in orientation angle [27]). In any case, the motion of the corner even with discontinuities in v iso can still be computed (see e.g. appendix A), but we focus for now on the simpler "continuous" case. Time variation in v iso (associated not only with the time evolution of path lengths s of material points but also with changes in Y coordinate) is also ignored for the purposes of computing the corner's instantaneous velocity, because points closely neighbouring the corner are liable to be consumed by that corner within a very short time, and v iso barely changes within that short time frame. To summarise, even though v iso does actually vary in space and time, in order to compute the local dynamics of the corner itself (which is our main focus here), it is possible to treat v iso as being effectively uniform and constant, at least for the purposes of computing the dynamics of the corner relative to nearby material points. Our immediate aim then is to compute how much the speed of the corner differs from this assumed locally uniform and constant v iso (which is the speed corresponding to "isotropic" material points).
Returning now to fig. 2(b), if we multiply the speed of material point C (denoted v aniso,C say, and corresponding to segment orientation α + and motion direction β + )b y the time interval τ , we obtain the length of BC. Likewise, if we multiply the speed of material point D (denoted v aniso,D say and corresponding to segment orientation α − and motion direction β − ) by the time interval τ ,w eo btain the length of BD. Since we have the slopes of BC and BD already (− tan β + and − tan β − respectively), we know (for any given τ ) exactly where the points C and D are located relative to B.
The slopes of the lines AC and AD are also known. Since the normals to these lines are oriented at α + θ/2 and α − θ/2 from the horizontal, their tangents are at π/2−α−θ/2andπ/2−α+θ/2 from the horizontal, and the slopes are 1/ tan(α+θ/2) and 1/ tan(α−θ/2) respectively. Given the slopes of these lines AC and AD, and given they pass through known coordinate locations C and D, their intersection, namely point A, can be determined (relative to point B). After some algebra, the coordinates of point A relative to B are found as: We define a speed v app (the apparent corner speed) as We are interested in the ratio v app /v iso .A l s oo fi n t e r e s ti s the ratio v app /v aniso,α where v aniso,α denotes the speed of a hypothetical material segment oriented at angle α from the horizontal: this segment of orientation α is referred to as "hypothetical" because, in the configuration of interest here, the orientation actually jumps discontinuously from α−θ/2toα+θ/2. The speed v aniso,α is given (analogously to eq. (14)) by where recall β α is defined by β α ≡ arctan(k v tan α). This is the direction in which this hypothetical material segment moves, and it turns out to be close to β, albeit 3 not exactly β.
We define a direction β app (the apparent direction of propagation of the corner) as and we also define (the orientation of a hypothetical material element that would move in the direction β app ). We are interested in how β app (the apparent direction of motion of the concave corner) compares with β (the bisector of the direction of motion of the material elements that are adjacent to the corner). To the extent that β is close to β α , this is essentially equivalent to asking how α app compares with α, but (to avoid repetition in what follows) we will focus on results for β app rather than for α app .
Another quantity of interest is the ratio v app /v aniso,hyp where v aniso,hyp is the speed of a hypothetical material point moving along the direction β app (a contrast with eq. (19), which indicates the speed of a hypothetical material point oriented in the direction α and hence moving along the direction β α ). We deduce (again completely analogously to eq. (14)) This completes our set of formulae for the motion of a concave corner. Before we can proceed to analyse these formulae however, we need to decide on the domain of α values to be analysed. This is discussed in the next subsection.

Domain of α under consideration
Note that fig. 2 considered a case where α>0, i.e. the components of the front normal point downward and to the right. As already mentioned earlier, this situation is believed to be typical because the upper parts of the foam front tend to advance more rapidly than the lower parts (a larger hydrostatic pressure opposes the front motion the deeper down in the reservoir that one moves). However in a heterogeneous reservoir it is possible that parts of the front in low permeability regions may be left behind relative to other parts of the front, and this may produce regions with α<0 (front normal pointing upwards and to the right, see e.g. sect. 6 later on). It is thus possible to contemplate a concave corner that occurs in a region with α<0.
It turns out however that it is possible to transform an α<0 configuration into an α>0 configuration simply by reflecting vertically about a horizontal plane. This transformation replaces the original α by −α, but preserves the original value of θ (indeed θ>0 is a requirement for the corner to be concave). It is now possible to proceed to compute values of β, ψ, x AB , y AB , v app , β app ,e t c .u s i n g the geometry corresponding to α>0 as discussed above. We can reflect back from the transformed configuration into the original one, which changes the sign of β, y AB and β app , but leaves ψ, x AB and v app unchanged. In effect this means that, for the purposes of computing the motion of the concave corner, it is sufficient to analyse the case α ≥ 0 which is what we consider below.
There is one further subtlety regarding the domain of α here: fig. 2 shows a case for which both angles α + θ/2 and α − θ/2 are less than π/2. For any given finite θ,i t is possible to contemplate a situation such that α + θ/2 exceeds π/2 but α − θ/2 remains less than π/2. Consider, e.g., the particular case for which π/2 − θ/2 <α<π / 2. Equations (10)-(11) still apply, but with the slight complication that now tan(α + θ/2) < 0. When we multiply this by k v and take the arctan in eqs. (10)-(11) we must select a solution branch between π/2a n dπ. This is of course only a minor complication, but in what follows we avoid it altogether by truncating the α domain to exclude values very close to π/2. In other words we consider 0 ≤ α ≤ π/2 − θ/2 instead of 0 ≤ α ≤ π/2. In practice truncating the domain by an amount θ/2 in this fashion is not at all restrictive because θ tends to be quite small (the corner tends to propagate in such a way as to stop θ from increasing, as we already mentioned earlier) and moreover in systems of practical interest (e.g. the case of anisotropic and heterogeneous reservoirs) concave corners are never observed for α values exceedingly close to π/2, but instead tend to develop at rather smaller α values. This is because (as we will see later on in sect. 6) the corners tend to develop in regions that lag behind the rest of the front, implying X is typically not too far away from a local minimum, so that segment orientations (given in sect. 2 as arctan(dX/dY )) and hence α itself are seldom exceedingly large.

Analysing the motion of a concave corner
We now analyse the behaviour of the above equations for various combinations of the parameters k v , α and θ (the values of β and ψ then being computed in terms of k v , α and θ).
In sect. 4.1 we examine specifically the quantities β and ψ. Then we look at the apparent speed v app of the concave corner starting with a simple symmetric case (sect. 4.2) and then more general cases (sect. 4.3). Next we consider the apparent direction of motion β app (sect. 4.4) with a special focus on the limit of small k v (sect. 4.5), giving also a geometrical interpretation of the implications of the results (sect. 4.6). Based on the geometrical interpretation, sect. 4.7 considers how a concave corner interacts with adjacent material points on the front, and we re-analyse the velocity fields (sect. 4.8) in light of this geometry.

Behaviour of β and ψ
Before we can analyse the motion of the corner (i.e. the behaviour of v app and β app ), we first need to understand the behaviour of β and ψ. The general behaviour of β and In a system where two segments of front meet at a concave corner, a plot of β (bisector of the direction of motion of the two adjacent segments) as a function of α (bisector of the orientation of the normals of the two adjacent segments). Both a weakly anisotropic case k v =0 .9 and a strongly anisotropic case k v =0.1 are shown. Cases are moreover considered where the angle θ through which the front turns at the corner is arbitrarily small, as well as a situation with finite θ = π/18: this latter case is only plotted for k v =0.1 (as for kv =0.9, any dependence on θ is exceedingly weak). For isotropic systems (with k v → 1), β is identical to α (regardless of θ).
ψ can be readily appreciated by focussing on the expressions in the small θ limit (given in eqs. (12) and (13)). These show that β and ψ are increasing functions of α, with β changing from 0 to π/2a sα itself changes from 0 to π/2. Meanwhile ψ over this same domain changes from k v θ to θ/k v (remember that k v is less than or equal to unity here). When k v is close to unity (weak anisotropy), the ratios β/α (see e.g. fig.3)andψ/θ (see e.g. fig. 4(a)) are roughly uniform (with very little sensitivity to α). However for k v ≪ 1 (strong anisotropy), significant deviations from uniformity of β/α and ψ/θ are observed (again see fig. 3 and also fig. 4(b)). In the domain α ≪ 1, linear approximations apply, giving β ≈ k v α,a n dψ/θ ≈ k v . In fact, when k v ≪ 1, we have β ≈ k v tan α (an order k v quantity) even for α values of order unity. The value of β only starts to increase significantly when tan α becomes order 1/k v which requires α values within an amount O(k v )o f π/2. In other words, when k v is small, the value of β remains small for almost all values of α, with significant increases in α only seen when α is very close to π/2.
The above analysis is restricted to the limit of arbitrarily small θ. There is a slight modification in the case of small but finite θ particularly in the case of strongly anisotropic systems, with k v ≪ 1. Recall that (as mentioned in sect. 3.3) we restrict consideration to values of α ≤ π/2 − θ/2 (implying that the maximum permitted value of α + is π/2). In this case we expect the maximum possible α − to be π/2 − θ. According to eqs. (10)- (11), and approximating tan α − for small values of θ, we deduce the maximum possible β is π/4+ 1 2 arctan(k v /θ), and the corresponding ψ is π/2 − arctan(k v /θ). In cases where In the case of two segments of front meeting at a concave corner, the ratio ψ/θ plotted against α,whereψ is the change in direction of motion from segment to segment, θ is the change in orientation from segment to segment (data here correspond in fact to the limit of arbitrarily small θ), and α is the orientation of the bisector of the two segments. Two weakly anisotropic cases are shown with k v =0 .9a n dkv =0 .98. (b) Analogous data for ψ/θ vs. α but now for strongly anisotropic systems k v =0 .01 and kv =0 .1: note the logarithmic scale. In addition to presenting data in the limit of arbitrarily small θ,ac a s ew i t hk v =0 .01 and finite θ = π/18 is also plotted. Nonetheless data for k v =0 .1w i t hθ = π/18 are not shown here since (on the scale of this graph), they would be difficult to distinguish from the situation where k v =0 .1 with arbitrarily small θ.
Clearly, based on the above discussion, there are some quantitative differences in the behaviour of β and ψ in the case of finite θ compared to the limiting case of arbitrarily small θ. However the qualitative behaviour is much the same, namely (regardless of the value of θ) both β and ψ/θ are increasing functions of α.
Moreover, when k v ≪ 1, there are strong non-uniformities in the values of β/α and ψ/θ. Specifically β/α and ψ/θ are much smaller than unity for most of the α range of interest, but show sudden increases near the largest values of α within the domain. Restricting the domain of α a little (i.e. considering not 0 ≤ α ≤ π/2 but instead 0 ≤ α ≤ π/2 − θ/2, as we do here) tends to truncate this region of sudden increases. Ordinarily θ is comparatively small (since the "speed up" rule for the concave corner is specifically designed to prevent the corner itself from sharpening and hence to prevent increases in the value of θ). Truncating the domain from 0 ≤ α ≤ π/2t o0≤ α ≤ π/2 − θ/2, has little bearing however on computing the motion any concave corners that might appear, as these typically occur at rather smaller α values where there is very little difference in β/α and/or ψ/θ between the case with arbitrarily small θ and that with small but finite θ. Hence the intuition gained from the small θ limiting case is actually very useful in practice.

Symmetric case α − = −α +
Having elucidated the general behaviour of β and ψ with respect to α, θ and k v , we now begin to consider the motion of the concave corner. To begin, we consider a symmetric case for which α = 0 and hence α − = −α + : two elements of front therefore join at a concave corner such that the bisector of the corner is horizontal. On symmetry grounds it follows that β = β app = α app =0a n d moreover β − = −β + . We deduce y AB = 0 (on symmetry grounds) and it follows from eqs. (16) and (18) that x AB = v iso τ/cos β + , and hence v app = v iso / cos β + .
Note that (since β =0)w eha v eβ + = ψ/2, and in cases of interest here with small values of θ, eq. (13) implies ψ ≈ k v θ. It follows therefore that the "speed up" factor v app /v iso is no longer 1/ cos(θ/2) (as per the isotropic case) but instead 1/ cos(k v θ/2). Provided θ is small, the speed up factor approximates to 1 + k 2 v θ 2 /8, which is closer to unity than was the case for isotropic systems (where the speed up factor approximates to 1 + θ 2 /8). Thus already we can see that a different "speed up" rule is needed in anisotropic systems compared to isotropic ones.
We continue to investigate how this "speed up" rule changes as a result of anisotropy in the following subsection.

Variation of v app with respect to α
In fig. 5(a) we show a graph, plotted as a function of α, of the apparent speed of the concave corner, v app (normalised by v iso which recall, in the neighbourhood of the  Here v app is the apparent speed of the concave corner (which is not itself a material point), v iso is speed of material points in an analogous isotropic system, and v aniso,α is the speed corresponding to a (hypothetical) material segment oriented along the bisector direction. Also shown is the ratio of speeds or "speed up" factor for the isotropic system 1/ cos(θ/2). (b) A zoomed view of the same data.
corner, is taken to be constant and uniform in the first instance) in the case k v =0 .9 with θ = π/18: these parameter values have been chosen arbitrarily in order to illustrate the model. For comparison, on fig. 5(a), we also show the constant value 1/ cos(θ/2), which is the "speed up" factor for the isotropic case. The main observation is that v app /v iso is quite close to unity (and hence quite close a l s ot o1 / cos(θ/2) with θ = π/18 here) in this k v =0 .9 case. This is unsurprising, given that in the limit k v → 1, we find that v app /v iso is identically 1/ cos(θ/2) regardless of α, a fact which can be demonstrated directly from eq. (18) by taking the k v → 1 limit (although we omit the proof here). Notice however that (for k v =0 .9) as α increases, v app /v iso falls. For very small α, v app /v iso starts off in fact greater than unity (consistent with the predictions of sect. 4.2) but it falls below unity for α greater than about 0.2 (see the zoomed view in fig. 5(b)). Interestingly however the function v app /v aniso,α is invariably in excess of unity, and in fact for α values greater than about 0.25, it even exceeds 1/ cos(θ/2). This means that were we to treat the concave corner as if it were oriented in the direction α, but with a velocity "sped up" relative to a material point by a factor 1/ cos(θ/2) (the isotropic speed up factor) we would actually be moving the point too slowly, at least for α greater than 0.25. The discrepancies in velocity are however never very large in this weakly anisotropic k v =0 .9 case: the largest value that v app /v aniso,α attains is around 1.008, compared with 1/ cos(θ/2) which is around 1.004.
For comparison, in fig. 6 we show data for k v = 0.1 (again the value is chosen arbitrarily here, but the anisotropy is now much stronger than before) still with θ = π/18 as above. Qualitatively we see similar features to what is seen in fig. 5(a) namely v app /v iso falls as α increases, whereas v app /v aniso,α grows to exceed 1/ cos(θ/2). However the overall scale of the graph is very different: in particular, as α grows, v app /v aniso,α reaches values in excess of 3. This implies that were we to move a concave corner with a speed up factor of just 1/ cos(θ/2) relative to a material point oriented in direction α we would in fact be moving the corner far too slowly, at least for larger α values.
The details of what is happening near α ≈ 0a r ed i fficult to resolve on the scale of fig. 6: in the interests of brevity, we have not provided zoomed views of this plot. W ea s s e r th o w e v e rt h a tf o rα → 0, v app /v iso is much closer to unity than 1/ cos(θ/2) is: it is in fact found to be 1+k 2 v θ 2 /8 with k v =0.1 here. As a result of lim α→0 v app being so close to unity, we only need comparatively small α values for v app /v iso to fall below unity: this happens already by α =0 .1. Moreover as α grows, v app /v aniso,α rises above 1/ cos(θ/2) comparatively quickly: this happens around α =0.2.
These fine details for comparatively small α will become relevant later on (see sect. 6). Our main point for now is that for certain α values, v app /v aniso,α is potentially significantly larger than unity, and likewise significantly larger than 1/ cos(θ/2). In other words, were we to speed up the concave corner by an amount of just 1/ cos(θ/2) relative to material points (instead of the correct amount v app /v aniso,α ), we would actually be moving the corner too slowly. A possible consequence is that the angle through which the corner turns continues to sharpen, and eventually the front shape could develop a spurious loop along the lines already discussed in the introduction (these spurious loops being a consequence of moving the corner too slowly).
To date we have considered just the speed of motion of the concave corner (compared to that of material points oriented in direction α). There is however another relevant question concerning whether the direction of motion of the concave corner is necessarily the same as that of a hypothetical material segment oriented in the direction α. In the isotropic case (k v → 1) this is certainly true, as the direction of motion of the concave corner bisects directions of motion of points on the adjacent material segments. The anisotropic case however is more subtle, as we shall see below. Note that β α is nearly halfway between β + and β − . The value of β app is however slightly less than β α .I td oe s however always fall between β + and β − at least for this weakly anisotropic case k v =0.9. As a first approximation then, it seems reasonable to assume that the concave corner moves in a direction roughly corresponding to β α (but moving at a slightly higher speed than a material point would move as per fig. 5(a)-(b)).

Variation of β app with respect to α
Analogous data for k v =0 .1( i.e. strong anisotropy) still with θ = π/18 are however somewhat surprising: see fig. 7(b). Over much of the domain of α, the value of β app is actually significantly below β α ( a n de v e nb e l o wβ − ). Hence, not only are there significant adjustments to the speed of the concave corner due to anisotropy (see fig. 6), but there are also significant adjustments to its direction of motion. The reasons for this are considered below.

Direction of concave corner motion for small k v
In the limit of small k v , the values of β, β + and β − ,a r e expected to be much smaller than unity over almost the By definition however β + − β − = ψ. Moreover to a good approximation, β ± ≈ k v tan α ± , and if in addition θ is small, we can approximate β + + β − by 2k v tan α which is essentially 2β. Likewise 2(tan α + − tan α − ) can be approximated by 2θ sec 2 α so that where eq. (13) has been used. Equation (20) then implies Since we are restricting consideration to small values of β (and hence to small values of k v tan α) we can approximate further In other words β app is (in this small k v limit), a factor k v times smaller than β (with β itself being on the order of k v times smaller than α). As a result, the concave corner is moving very nearly horizontally. It follows that β app is not only less than β, but also can be less than β − . Note in particular that β − becomes greater than zero once α − itself becomes greater than zero (or equivalently when α exceeds θ/2). It follows that β app given by β app ≈ k v β (with k v ≪ 1) can be smaller than β − for α values just very slightly in excess of θ/2 (and θ/2 is itself typically quite small, e.g. it is θ/2=π/36 for our data).

Geometrical interpretation for small k v case
The above observations (in particular the direction of apparent corner motion β app being smaller than β − )c a nbe interpreted geometrically as follows. Consider a concave corner where the segment orientation jumps from α + to α − , and the direction of motion of material points jumps from β + to β − . Were we to "round off" this concave corner into a region of large but finite curvature, we would encounter material segments with all possible orientations between α + and α − and moving in all possible directions between β + and β − . However the motion of the concave corner itself (prior to being "rounded off") is not in any of these directions between β + and β − , and does not correspond to any orientation between α + and α − .
There is no contradiction here: the concave corner is not a material point, and so need not move in a direction that corresponds to any of the material elements in the "rounded off" corner. The result is nonetheless surprising when seen in the context of the isotropic system, for which the direction of motion of the corner always corresponds to the direction of motion of the material points at the centre of the rounded off corner.
Another consequence of β app being smaller than β − is that the definition sketch of fig. 2 must be modified: see fig. 8. In fig. 8, material point D (referred to more correctly in this case, for reasons to be explained shortly, as "virtual" material point D) is no longer "below" the concave corner A (as was formerly the case in fig. 2) but now is higher up. In other words, the slope of the trajectory AB of the corner is smaller in magnitude than the slope of the trajectory DB of material point D (these magnitudes being tan β app and tan β − , respectively).
The fact that the system adopts a configuration such as this has some additional important implications (as we consider in the next subsection).  Fig. 8. Material segment AC meets material segment AE at the concave corner A, the segment orientation turning through an angle θ at the corner. Over a time τ , material point C migrates to a new location B along a line with slope − tan β +. During this same time τ , the concave corner which is not itself a material point, also migrates to B along a line with slope − tan β app. All the material points between A and C are thereby consumed. Meanwhile material point E migrates t oan e wp o s i t i o nF along a line with slope − tan β −.T h e distance BF is however greater than AE. Even though material point E has a small velocity component in the direction from E to A, the concave corner has an even larger velocity component in this direction, and so "outruns" the motion of point E. This means that new material points must continually be "extracted" from the concave corner. In particular the new material point which is "extracted" or created at time τ , can be considered to have originated from a "virtual" material p o i n ta tl o c a t i o nD (on an extrapolation of the segment AE), with the additional property that DB is parallel to EF. In fact all the "virtual" material points on the segment AD between A and D need to be extracted from the corner at some stage during the time interval τ .

Consumption vs. creation of material points
In the case of an isotropic system, whenever two material segments turn through an angle θ at a concave corner, it is known [19] that the motion of the concave corner consumes nearby material points. Specifically the motion of the concave corner bisects the angle at which the two material segments meet. This means that the motion of the concave corner has a component that is tangent to each of the material segments. Material points themselves have no tangential motion in an isotropic case, moving instead purely along segment normals. Thus the tangentially moving corner invariably manages to consume the tangentially fixed material points.
The case of a weakly anisotropic system as shown in fig. 2 is just a slight generalization of what happens in the above mentioned isotropic case. The main complication is that the velocities of the material points now also have a component along the tangent to the material element. For example, in fig. 2(b), material point C (which is moving in the direction CB) has a component of its motion parallel to the segment AC. Nonetheless the tangential motion of the concave corner A (which is moving in the direction AB) actually exceeds that of the material point C implying that the corner still manages to consume ma-terial point C as time proceeds (an event which occurs at location B).
The case of material point D is more subtle. The motion of concave corner A (moving along AB as mentioned above) still has a component tangent to the segment AD but this component can be either in the direction A to D or from D to A, depending on the exact level of anisotropy. The motion of material point D (moving in the direction DB) however very definitely has a tangential component directed from D to A: this material point therefore migrates towards the corner, allowing the corner to consume it (again an event which occurs at location B).
Contrast this with the case of a strongly anisotropic system in fig. 8. Whilst the concave corner A still manages to consume material point C in the same manner as described above, the situation on the segment denoted AE is rather different. Material point E indicated on fig. 8 has a tangential component of its motion directed from E to A, but the concave corner at A has an even larger tangential component in this same direction, and so "outruns" the material point. By the time the concave corner reaches location B, the material point E has attained location F , even further away from the corner than where it started. As a consequence therefore, new material points need to be created or "extracted" from the corner.
These newly created points can be considered to have originated from "virtual" material points which are distributed along a segment AD s h o w ni nfi g .8w h i c hi s initially behind the foam front. The point D for instance ceases to be "virtual" only once the corner reaches location B. Once extracted from the corner, the material point in question then migrates away from the corner. There is an implication for a numerical simulation of pressure-driven growth that new numerical grid points continually need to be created in the neighbourhood of the corner. Whilst the need to add material points to simulations of pressuredriven growth is not in itself unusual (indeed it is standard practice in the case of convex foam front shapes [19]), this is the first situation where we have encountered a need to add new material points to a pressure-driven growth system for a concave shape. This would moreover never happen in isotropic systems, where concavities always consume material points, but never create them.
The need to create new material points in particular on the segment that is oriented with its normal closer to the horizontal (but not on the segment oriented with its normal further from the horizontal) seems to be associated with the former segment being the faster moving of the two. In an extreme case for instance in which material points on segment AE move at the maximum allowed speed v iso and material points on segment AC do not move at all, the direction of motion of the concave corner is well defined: the corner is confined to the existing segment AC. Given that the material elements turn through an angle θ at the concave corner, it is then an exercise in elementary trigonometry to compute the speed of the concave corner (which is v iso cosec θ in this special case) and the speed at which new material points are extracted (which is v iso cot θ). We emphasise however that these results only  v app /v iso v aniso,hyp /v iso Fig. 9. Case of a strongly anisotropic system kv =0 .1w i t h the front turning through θ = π/18 at a concave corner. The velocity ratios v app/viso and v aniso,hyp /viso are plotted against orientation α of the corner bisector, v app being the apparent speed of the concave corner, v aniso,hyp being the material point speed on a hypothetical material segment travelling in the same direction as the concave corner and v iso is the material point speed in an analogous isotropic system. apply in this extreme (and somewhat artificial) case in which one material segment moves and the other is fixed in position. Generally speaking both material segments move, albeit at different rates, and the trajectory of the corner must reflect that.

Re-analysis of velocities at the corner
There is one further consequence of the angle β app being outside the range β − through β + in the case of strongly anisotropic systems: this consequence is related to speeds of motion. Consider a comparison between v app (the apparent speed of the corner) with v aniso,hyp (defined in sect. 3 as the speed of a hypothetical material point moving in the same direction β app as the corner). Results are shown for the case k v =0 .1 in fig. 9. Observe that (counter-intuitively) there are α values such that the corner can move more slowly than the (hypothetical) material point. For k v =0.1, this actually happens for α values greater than about 0.1, but on the scale of the figure, it is easiest to see for α values greater than about 0.7.
At first sight this appears problematic, since intuition from the isotropic case suggests that a sharp concave corner moves faster than any material point contained in a "rounded off" version of the corner (where the two material segments are joined by an arc of large but finite curvature). There is however no contradiction here: in the case of a strongly anisotropic system with k v much smaller than unity, the "rounded off" version of the corner only contains material points moving in directions between β − and β + : it does not contain any faster-moving material points travelling in the direction β app . The corner (moving at speed v app ) is actually travelling faster than material points moving in directions between β − and β + (see e.g. the curve labelled v app /v aniso,α in fig. 6 which corresponds specifically to motion in the direction β α , with β α itself between β − and β + according to fig. 7(b)).
This completes our analysis of how the "speed up" rule of a concave corner is affected by anisotropy (and in particular how the corner's motion differs from the isotropic case). In what follows, we aim to implement this within a numerical code to simulate a propagating foam front via the pressure-driven growth model. However, rather than analysing a concave corner oriented at a specified predetermined angle (as has been done above), what we will do instead is study a system that is both heterogeneous and anisotropic, allowing the concave corners to develop naturally as a result of the heterogeneity, with whatever orientation that the dynamics producing the corner demands.

Numerical scheme: Pressure-driven growth
The numerical scheme for simulating isotropic systems via pressure-driven growth is already well established in the literature [11,19,22], so we give here just a very brief description of how the scheme generalizes to a heterogeneous and anisotropic case. For the purposes of simulating a heterogeneous and anisotropic reservoir we chose parameter values k het =0.3a n dn het = 3 (heterogeneity; see eq. (8)) and k v =0 .1 (anisotropy; see eq. (7)). The above values are chosen arbitrarily for the purposes of illustrating the model. We adopted the X, Y coordinate system from sect. 2. We discretised the foam front initially into 200 material points. We assumed that these were initially equally spaced vertically over the interval 0 ≤ Y ≤ 1. Initially X and s values were set to a uniform value s 0 with s 0 =0.001. Choosing a small but non-zero s 0 ensures that velocities (which are proportional to s −1 ) are finite at time t =0.
Choosing a time step δt =1 0 −5 , we propagated the material points according to eqs. (6)- (7): the distances between adjacent material points tends to change over time. We added new material points whenever segments grew to be greater than 0.02 and removed material points whenever segments shrunk to be smaller than 0.002. We also moved the top boundary point denoted X top according to X top =(2t + s 2 0 ) 1/2 which turns out to satisfy our boundary condition that the front motion at the top is parallel to the boundary itself.
As the front propagated, concavities were found to develop upon it, typically in the low permeability regions for this heterogeneous system. These concavities became increasingly tightly curved over time. We considered that a sharp concave corner had formed once the angle between adjacent front elements reached the value θ = π/18. As before, this value is chosen arbitrarily for the purposes of illustrating the model, but similar results are expected [19] for any "small" θ value provided θ ≪ π.
We could determine v iso specifically at the sharp corner as it formed (v iso is equal to Ys −1 J(Y ) in this case). Moreover the angle α at the concave corner is defined as  Fig. 10. The predicted shape of a foam front for a heterogeneous, anisotropic system: k het =0 .3, n het = 3 and kv =0 .1. Concave corners are formed, but points in the concavity continue to be propagated as if they were material points without any speed up factor, which leads to spurious loops (see inset) after time approximately t ≈ 8.
the angle below the horizontal of the bisector of the two material segments that meet there.
Using the theory in sect. 3, we could then determine v aniso,α and v app (via the ratios v aniso,α /v iso and v app /v iso ), as well as the angles β α and β app .
The concave corner once formed of course might need to be propagated differently from neighbouring material points. We looked at several different options here.
1. The corner was propagated (incorrectly) exactly as if it were a material point, i.e. at speed v aniso,α and at an angle β α from the horizontal.
2. The corner was propagated (again incorrectly) at a speed 1/ cos(θ/2) times v aniso,α at an angle β α from the horizontal (corresponding to what would happen in an isotropic system).
3. The corner was propagated at the correct speed v app and at the correct angle β app from the horizontal for an anisotropic system.
In the following section, we consider the numerical consequences of moving the concave corner according to the various different rules above.

Numerical results
Numerical data for the computed front shape at selected times are shown in fig. 10. As time proceeds, a concave corner is seen to form. In this particular figure, we do not apply any speed up rule whatsoever once the concave corner forms. Spurious loops develop as a result: already by time t = 8 they have begun to form (see the zoomed inset), and they grow as time proceeds.
In fig. 11 we show data where the concave corner is sped up after it forms (i.e. once the angle θ through which the front turns at the corner exceeds π/18). In   fig. 10, but moving the concavity faster than material points, as soon as the angle θ between adjacent material segments exceeds π/18. Specifically the corner is propagated as it would move in an isotropic system using a speed 1/ cos(θ/2) higher than the material point speed, albeit with the direction of motion unchanged. Spurious loops do not now appear. (b) Foam front shape propagating the concave corner using the correct anisotropic propagation rules, i.e. with the correct apparent speed v app and in the correct apparent direction βapp (see formulae in sect. 3). Again spurious loops do not appear, but the corners are maintained comparatively sharp. fig. 11(a) the speed up factor relative to material points is 1/ cos(θ/2), which is the speed up rule applicable to an isotropic system rather than an anisotropic one. The corner still moves in the same direction as material points do, i.e. it moves in direction β α when its bisector is oriented in direction α. This is sufficient to eliminate the formation of any spurious loops. In fig. 11(b) we show data for the case where the concave corner propagates instead at speed v app in a direction β app : these are the correct anisotropic speed up rules. Again no spurious loops appear.
Comparing fig. 11(a) and fig. 11(b) we can however see subtle differences between them. The concavities in the latter figure are sharp, but in the former figure they are rounded off. This is a reflection of the fact that the speed assigned to the corner in fig. 11(b) is actually less than that in fig. 11(a), and the lower the speed, the sharper the corner. At first sight fig. 11(b) having a lower propagation speed appears surprising since fig. 5 and fig. 6 for instance show that over most of the domain of α values the anisotropic speed up factor v app /v aniso,α exceeds the isotropic one 1/ cos(θ/2). The resolution of this apparent discrepancy however arises from the fact that all the concavities in figs. 10-11 form at small values of α. Specifically they appear in regions where the permeability is low, meaning those parts of the front lag behind the rest. At these points which lag behind the rest of the front, the value of X as a function of Y is a local minimum, so the front tangent is vertical and the front normal is horizontal, implying α is low (α ≪ 1). The values of β α and β app are likewise very low, so the corner is moving nearly horizontally regardless of which speed we assign to it. For small α however, based on the analysis in sect. 4.2, it follows that the anisotropic speed up factor approximates to 1/ cos(ψ/2) in contrast with the isotropic factor 1/ cos(θ/2). The anisotropic speed up is then slightly smaller since ψ ≈ k v θ with k v < 1.
We have shown therefore that in a heterogeneous and anisotropic system, it is necessary to take specific account of the anisotropy in order to propagate a concave corner in such a way as to keep it sharp (but still not so sharp as to form spurious loops). The results shown here are of course just preliminary ones demonstrating the feasibility of our algorithm for propagating concavities. A thorough parametric study of the behaviour of the heterogeneous case, both isotropic and anisotropic, in terms of the various system parameters k het , n het and k v still remains for future work.

Conclusions
In the context of a pressure-driven growth model, we have considered the rule for propagating a concave corner in the shape of a foam front, such as might occur during improved oil recovery within a heterogeneous and anisotropic oil reservoir. The concave corner needs to propagate more quickly than neighbouring material points on the foam front. If it were propagated (incorrectly) at the same speed as those neighbouring material points, (undesired) spurious loops would develop in the numerically computed front shape.
Whilst this situation also occurs in the case of a concave corner that develops in isotropic systems, the corner "speed up" factor for anisotropic cases may need to be rather different from that for isotropic ones. Depending on the orientation of the concave corner, the anisotropic "speed up" factor determined here could be significantly larger than the isotropic one. In that case, applying just the isotropic "speed up" factor to the anisotropic systems, still implies a risk of producing spurious loops. Moreover switching from isotropic systems to anisotropic ones affects not just the speed of the concave corner, but also the direction in which it moves. Specifically it can propagate in a direction that is much closer to the horizontal than the direction of motion of any neighbouring material points. For a numerical scheme computing the front shape, this also means that new material points may need to be created in the neighbourhood of the concave corner (rather than material points being consumed by it, as is the norm for isotropic systems).
One way that fronts are expected to develop concavities is by having a permeability that is heterogeneous varying with depth. Concavities then tend to develop in low permeability regions, such that the foam front in these regions is lagging behind the remainder of the front. The concave corners are then oriented such that the corner bisector is nearly horizontal, and the motion of the corner is likewise nearly horizontal. Some of the more counterintuitive behaviours alluded to above (e.g. concave corners causing material points to be created rather than consumed) are now avoided. Indeed this is a special case in which the anisotropic speed up factor is actually less than the isotropic one, so spurious loops are avoided using either of the two speed up factors. Using the slightly smaller anisotropic speed up factor however keeps the concave corners sharp.
Having now developed here a numerical scheme for pressure-driven growth that can handle systems which are anisotropic (and more generally systems which are both heterogeneous and anisotropic), a thorough parametric study of the foam front evolution in terms of the parameters controlling heterogeneity and anisotropy remains for future work. We anticipate that anisotropy will lead to increased gravity "override": when foam first arrives at a production well, more oil will be left in place in anisotropic systems than in analogous isotropic ones. Heterogeneity may compound this, since (as mentioned above) it is the introduction of heterogeneity that causes the formation of concave corners. These corners then lag well behind the leading part of the front (again implying significant amounts of oil having been left in place between the concave corner and the leading part of the front).
PG acknowledges funding from CONICYT Chile (folio 80140040). PG and CTU acknowledge funding from FONDE-CYT project number 1120587. EMH acknowledges scholarship funding from the EPS-CONACyT programme.

Appendix A. Vector formula for motion of the concave corner
This appendix presents an alternative derivation of the formula for the motion of a concave corner that is equivalent to the presentation considered in the main text. Whereas the presentation in the main text was convenient because it was particularly amenable to geometric interpretation, the presentation here is arguably less easy to interpret geometrically, but has a compact vector representation, making it easy to implement in a computer program.
Consider the points A, B, C and D as defined in either fig. 2 or fig. 8: the derivation works for either configuration. Suppose that the velocity vector of the concave corner denoted v A can be written in vector form as v A = v AC n AC + v AD n AD , (A.1) where n AC and n AD are normals to front elements, and v AC and v AD are coefficients that we must find. Projecting onto n AC , and equating to n AC ·v C (v C being the velocity of material point C) v aniso,C cos(α + − β + )=v AC + v AD cos θ, where v aniso,C denotes the speed of the material point C, and the factor cos(α + − β + ) represents projection, the motion of point C being in direction β + and the normal to AC being in direction α + . The angle θ represents the jump in direction between segments AC and AD.I nt h e above equation, all terms except v AC and v AD are known.
Geometrically the equation ensures that (despite the corner A not actually being a material point) the segment from A to C propagates in the same way as a material segment would, without the motion of the corner introducing any artificial rotation. Similarly projecting onto n AD , and equating to n AD · v D (v D being the velocity of material point D) where eq. (14) has been used with v iso being the local material point speed in an analogous isotropic system. Equation (A.4) as written remains valid even if v iso exhibits a discontinuous jump 4 at A, leading to different v iso values at points C and D, although that is a complication we do not consider any further here, treating instead v iso as being uniform at least locally.