Issue 
Sci. Tech. Energ. Transition
Volume 78, 2023
Characterization and Modeling of the Subsurface in the Context of Ecological Transition



Article Number  42  
Number of page(s)  40  
DOI  https://doi.org/10.2516/stet/2023036  
Published online  22 December 2023 
Regular Article
Foam placement for soil remediation: scaling foam flow models in heterogeneous porous media for a better improvement of sweep efficiency^{}
^{1}
IFP Energies nouvelles, Earth Sciences and Environmental Technologies Division, 1 et 4 avenue de Bois Préau, 92852 RueilMalmaison Cedex, France
^{2}
IFP Energies nouvelles, Digital Sciences and Technologies Division, 1 et 4 avenue de Bois Préau, 92852 RueilMalmaison Cedex, France
^{*} Corresponding author: frederic.douarche@ifpen.fr
Received:
2
September
2022
Accepted:
10
November
2023
If injected with a large gas fraction, foam reduces mobility more in highpermeability layers and diverts flow to lowpermeability layers. Here is a qualitative statement that has been claimed many times in many works related to environmental remediation and oil recovery. It is so true and relevant for foam flow in porous media and yet so little quantified and even less exploited in Darcyscale numerical simulation. After briefly reviewing opportunities and challenges related to the use of foams in porous media and its Darcyscale implicittexture and populationbalance modelling, we make a detour out of the strict framework of mathematical models by revisiting with a fresh eye the physics of foams on the large scale of heterogeneous natural porous media in terms of scaling laws. Specifically, it has been recently shown experimentally and theoretically that foam mobility reduction scales approximately as the square root of rock permeability within the framework of Darcytype implicittexture foam flow models [Douarche et al. (2020) Scaling foam flow models in heterogeneous reservoirs for a better improvement of sweep efficiency (Paper ThB04), in:17th European Conference of the Mathematics of Oil Recovery (ECMOR), Edinburgh, Scotland, 14–17 September; Gassara et al. (2020) Trans. Porous Media 131, 1, 193–221]. This also appears to hold for populationbalance models under the local steady state assumption. This quantitative scaling law for the effect of permeability on foam properties was inferred from an analogy between foam flow in porous media and foam flow in capillary tubes and was found consistent with the modelling of available experimental data. We demonstrate by varying the permeability contrast and anisotropy of a stack of porous layers how foam regulates the flow and leads to flow diversion from high to lowpermeability layers. The threshold in permeability heterogeneity for which such a foamdriven diversion becomes effective is quantified through a sensitivity study accounting for foam injection type, permeability heterogeneity and anisotropy, heterogeneity structure, and scaling procedure. The simulations carried out clearly indicate that ignoring mobility reduction dependence on permeability in the foam process assessment of heterogeneous formations leads to an underestimation of mobility reduction benefits to improve flow conformance. The question of cores selection, as this rocktyping strategy in relation to the porous medium characterization may guide a smart and optimal design of prefeasibility laboratory campaign for foam evaluation, and the generalization of the findings to multifacies geological formations are also discussed. As such, the use of physical foam mobility reduction scaling law is highly recommended for foam process evaluation.
Key words: Multiphase flow / Porous media / Foam / Scaling laws / Heterogeneity / Models / Simulation / Soil remediation
Foreword: This paper includes a large Appendix A.2 which reports some selected insitu foam flow dynamic behaviour. These figures are discussed in the main text but have been put in appendix to facilitate the reading. Its consultation remains optional, but the authors believe the reader will find there a good complement to the production data overall response of the porous medium under foam flooding.
© The Author(s), published by EDP Sciences, 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Understanding and predicting foam flow behaviour in heterogeneous porous media, matricial or naturally fractured, is a generic and extremely complex problem which is often encountered in geosciences applications involving the displacement of a nonaqueous phase liquid (NAPL) by water and/or gas [1–6], given that the behaviour of foams in porous media is far from being elucidated and remains a subject of active research. Whether the targeted application is soil remediation, carbon capture, utilization and storage or hydrocarbon production, foam suffers from two major fundamental limitations: surfactant adsorption on rock and the antifoaming character of the NAPL to displace, as is the case of most pollutants and oil.
The result of these two strong constraints, if one considers foam as a NAPL frontal displacing fluid, is a fundamental contradiction that is very difficult to resolve in practice between the need to limit the amount of surfactants injected for obvious economic reasons [2, 4] and the dual penalty foam suffers from, namely not only chemicals adsorption on rock but also the antifoaming character of the NAPL to displace that yields additional losses and requires even more chemicals to compensate for. In this regard, the use of foams for CO_{2} storage in aquifer does not present the antifoaming NAPL issue but is more questionable technically unless foam improves the effective storage capacity by smoothing out viscous fingering that may develop between CO_{2} and brine and/or channeling due to heterogeneity through conformance control [7], provided foams are stable over the longterm.
In this context, one can only hope that new foaming formulations that are cheaper, more resistant to NAPL, less adsorbent, more stable over the long term, and more respectful of the environment will be developed. Applicationwise, the story could end there, but hopefully, in addition to their large apparent viscosity which could have been very effective in frontal displacement, unfortunately impacted by the antifoaming nature of fluids such as NAPL and oil, foams have attractive selectivity properties in porous media which may change the game and are the subject of this paper.
In the perspective of a largescale deployment of this type of injection, the investments necessary for the implementation of such a process on geological formations of hectometric or kilometric size and the associated risks require predictive models making it possible to determine optimal injection configurations that are robust to the inherent subsurface uncertainties beforehand. One must admit it is illusory, at this geological scale under consideration, to aspire to model the physics of foams in porous media at the pore scale, even assuming it would be completely known at this scale. Indeed, if the numerical calculation allowed it one day, on the condition of a fully known physics, the problem of the time and lengthscales that can be sampled for subsurface characterization hence the setup of the initial boundary value problem, and the resulting unavoidable uncertainties would remain. Even if we could integrate the porescale motion equations in a general form, it would be completely impossible to substitute in the general solution the (uncertain) initial conditions for the velocities, pressures, saturations and compositions of all the phases at reservoir scale. As an example, considering a typical pore size of about 50 μm, a porous medium of 20% porosity, 200 × 200 m extension and 50 m thickness counts very roughly about a billion pores.
It goes without saying that high performance computing and rocksolid physical laws are useless without an accurate characterization of the porous medium under consideration, which certainly does not extend to the scale of the cubic centimeter and even less of the pore over an entire geological formation. Thus, if a reasonable and convenient modelling tradeoff is to build semiempirical augmented generalized multiphase Darcytype models to describe coarsely foam flow as we shall see in the following, a certain void remains between the pore scale and that already homogenized of the fictitious continuous medium considered at the Darcy’s scale.
Given the atrocious complexity of the problem, it is likely that no answer will be provided anytime soon regarding the upscaling of foam flow within a multiphase framework from the pore to a representative geological scale. In industrial practice, either one resorts to casebycase evaluations that involve massive laboratory studies sampling the heterogeneity variability of the porous medium under consideration, which are necessarily extremely expensive (therefore extremely rare and most often proprietary), or goes with the wet finger hunch with very few measurements and some literature survey (win or loose strategy when the risks allow it). One of the major obstacles lies in the absolutely phenomenal amount of work (and costs) that it represents to engage in feasibility studies without having the certainty of choosing foams as a displacement fluid at the end.
This work transposes a knowledge coming from reservoir engineering to environmental engineering, with however differences between those two application fields, mostly (i) in term of permeability that is usually much higher for aquifers, and (ii) in term of air that can be present in the transition zone of nearsubsurface undersaturated soils. The thermodynamic conditions considered are representative of deep subsurface porous media, not of shallow aquifers, but have little impact on the results of this work that is mainly phenomenological in nature. This work represents a first conceptual attempt and does not claim to cover all possible cases.
We show the contribution of a sensible physical coarsegrain approach based on scaling laws that is resolutely pragmatic and which has the advantage of being neither simplistic nor involving massive workflows which are known to be difficult to implement in practice. From a workflow perspective, as we shall see, the authors are not suggesting any additional complexity but rather a smart and optimal guidance for foam evaluation in relation to the porous medium characterization on which precisely this type of approach opens. The main result, which lies in a quantitative permeabilitydependent effective foam mobility, is quite directly and easily exploitable in the framework of augmented Darcytype models addressing largescale heterogeneous porous media, and indeed comply with the “if injected with a large gas fraction, foam reduces mobility more in highpermeability layers and diverts flow to lowpermeability layers” saying, of which we illustrate quantitatively and as concretely as can be demonstrated by a representative generic minimal model what it is about and the implications.
To sum up, it is an understatement to say that foam flow modelling in porous media still has a long way to go. Until further progress, let us see what works well with geosciences foam flow modelling despite the dizzying change of scale and the hocuspocus. The paper is organized as follows: first the generic problem to be solved is presented in Section 2 where the notations and driving equations are introduced. In particular it is shown that the emergence of the scaling of foam maximum mobility reduction as the square root of permeability is not limited to a specific Darcytype model (implicittexture or populationbalance). Section 3 is devoted to the setup of a minimal model heterogeneous porous medium that is used in Sections 4 (detailed insitu flow behaviour) and 5 (production data analysis) to study the impact on the flow performance of taking into account such a scaling law by varying the porous medium degree of heterogeneity. The results of a sensitivity analysis are summarized in Section 6. The paper concludes with a discussion in Section 7. A first Appendix A provides elements for estimating the parameters of the considered foam model from measurements while a second Appendix A.2 reports the time evolution of simulated saturation and composition fields related to foam for a few selected configurations which cannot fit in the main text without harming its reading but are however essential to the understanding of foam flow and of the arguments advanced in this work. The simulations presented in this work were carried out with a research version of IFP Energies nouvelles Puma reservoir simulator [8] in which new automatic scaling functionalities have been developed for the only purpose of demonstrating the impact of foam selectivity. Let us first recall the different types of generic foam models considered in geosciences and the scaling law that we propose to test.
2 Darcyscale model types, driving equations and scaling relationships
2.1 Darcyscale model types and driving equations
We start with general considerations about threephase flow in porous media in the presence of foam. We distinguish three phases: an aqueous phase denoted w, a NAPL phase o and a gas phase g. This flow is modified by the presence of foam. Modelling foam requires the presence of surfactant, which is transported by the water phase and requires to solve an additional mass balance equation. Surfactant is either mobile or adsorbed on the rock.
It has been shown that the transport of liquid is not affected by the presence of foam [9–11]. On the opposite, the gas velocity is significantly reduced by the presence of foam. Thus, to describe the water and nonaqueous phases, we consider a “blackNAPL” model (in reference to the “blackoil” model developed in hydrocarbon context [12, 13]) where the gas phase involves a modified velocity denoted . The mass conservation equations read:(1)where Φ is the rock porosity and the fluxes J_{w}, , and are defined as(2)
For each phase denoted φ = w, o, g, U_{φ} is the Darcy velocity, S_{φ} the saturation, ρ_{φ} the density and q_{φ} the source/sink term per unit volume of porous medium. The gas phase contains a single volatile component denoted v, whereas the NAPL phase contains a heavy component denoted h and previous volatile component v, with C_{h} and C_{v} = 1 − C_{h} their respective mass fraction in the NAPL phase. The equilibrium constant K_{v} = 1/X_{v}, with X_{v} the molar fraction of volatile component in the NAPL phase, is defined from molar masses of v and h, and mass fractions C_{v} and C_{h}. K_{v} is a function of pressure. stands for the flowing surfactant mass fraction in the water phase and for the adsorbed surfactant mass fraction on the rock with ρ_{r} the rock mass density. Mobile and adsorbed surfactant mass fractions are related with an adsorption law such as the Langmuir isotherm.
Under creeping (lowvelocity) flow conditions, the phase velocities in permeable porous media are governed by the generalized multiphase Darcy’s law [14]:(3)where U_{φ} is the Darcy velocity, K the porous medium permeability tensor, λ_{φ} = k_{rφ}/μ_{φ} the mobility, μ_{φ} the viscosity, p_{φ} the pressure and V_{φ} = U_{φ}/(ΦS_{φ}) the interstitial velocity of phase φ [14].
The flow of gas in the presence of foam is modelled differently whether a PopulationBalance (PB) or a SemiEmpirical (SE)/ImplicitTexture (IT) modelling approach is used. These models are based on Darcytype laws which are extended to obtain a modified gas velocity. Most SE/IT foam flow models apply a gas mobility reduction factor when surfactant is present in the water phase. Specifically, the gas mobility of IT model is scaled by a multiparameter interpolation function FM assigned to the relative permeability to gas, whereas the gas viscosity is assumed unchanged whether foam is present or not:(4)where is the modified relative permeability to gas for the IT model in the presence of foam and k_{rg} is the conventional relative permeability to gas. Of course, this mobility reduction can also be interpreted as an increase in the viscosity of the gas in the form of foam, as we shall see hereafter for PB models. FM is a multiparameter interpolation functional form that includes the contributions of physical parameters impacting the gas mobility reduction. FM is formulated as follows:(5)where M_{ref} is the reference (maximum) gas mobility reduction under optimal conditions of the rockfluidadditive system under consideration, and F_{i} are functions of four physical parameters that are surfactant concentration, water and NAPL saturations, gas velocity, or equivalently, capillary number. Following [15–17] these functions may be written(6)where is the capillary number and σ_{wg} the interfacial tension between water and gas. Positive numbers e_{s}, Θ, e_{o} and e_{c} drive the shape of F_{1}, F_{2}, F_{3} and F_{4} functions, respectively, while , , and are thresholds. Specifically:

M_{ref} can be estimated from the ratio of steadystate measurements of differential pressure across coreplugs, Δp_{f}/Δp, with and without foam, as measured in the laboratory, given that Δp corresponds to a gas injection in the IT model according to (4).

F_{1} function accounts for lamellas stability that is obtained for surfactant concentrations larger than some threshold [18, 19].

F_{2} function governs the sharpness of the transition from the low to the highquality regime, through the dimensionless constant Θ, as water saturation decreases in the vicinity of [15–17].

F_{3} function accounts for the impact of the antifoaming character of NAPL on foam.

F_{4} function accounts for foam shearthinning (or shearthickening) behaviour through the capillary number Ca.
Examples of such functions can be found in [15–17, 20] and cited references therein. Estimation of M_{ref} from measurements is detailed in Appendix A. Chemicals natural degradation with time can be accounted for through a halflife time associated to an exponential decay.
The threephase flow relative permeability to NAPL is classically modelled as a combination of the underlying water/NAPL and gas/NAPL twophase flow subsystems following Stone’s models [21, 22]. Eventually, the system (1)–(6) with closure relationships ∑S_{φ} = 1, p_{cow}(S_{w}) = p_{o} − p_{w} and p_{cog}(S_{g}) = p_{g} − p_{w}, where p_{cow} (p_{cog}) is the capillary pressure between NAPL and water (gas) phases, is solved for the initial boundary value problem in the unknowns p_{φ}(x, t), S_{φ}(x, t) and .
Populationbalance models share the same generalized multiphase Darcy framework but consider foam texture n_{f}, defined as the number of lamellas per unit of foam volume, as an additional unknown to be solved, and a foamed gas viscosity constitutive relationship which depends on n_{f} and adjustable parameters, based on Hirasaki et al [23] extension of Bretherton [24] work to porous media. While it has been shown that PB models are equivalent to IT ones at local steady state [16, 17], they obviously have the advantage, at least in principle (we come back to this below), of being able to better model foam transient generation and collapse which is driven by source terms in the texture transport equation (7) given hereafter. Reviews of such PB models, that are not addressed in this work, can be found in [25, 26]. Here are the main useful features for the work presented here.
Foam mobility is strongly related to its texture which is a key variable in foam modelling in porous media: as foam texture increases, the resistance to gas flow in porous media increases. Populationbalance models were designed in order to relate explicitly the gas mobility reduction to the foam texture. Thus, for PB models, the dynamics of foam texture needs to be modelled and a lamella populationbalance equation is considered. The impact of foam texture on gas mobility is modelled through an effective gas viscosity and a gas trapping for more elaborated models [11, 27, 28]. The populationbalance model involves a lamellas balance that includes lamellas advection at the modified gas velocity and source/sink terms taking into account lamellas creation and destruction [29, 30]. The lamellas population balance equation reads(7)where n_{f} is the foam texture, that is the number of flowing lamellas (or foam bubbles) per unit volume of gas and q_{f} is the external lamellas source/sink term (number of lamellas per unit of time and per unit volume of porous medium). r_{g} and r_{c} are the rates of lamellas creation and coalescence. Different formulas have been proposed in the PB models literature to model these rates [28, 31, 32].
As mentioned earlier, the gas velocity is significantly reduced when foam develops. The viscositytype PB models extend the generalized Darcy’s equation for the gas phase as follows:(8)where is the effective gas viscosity when flowing as a succession of bubbles. This relationship is similar to (4) but its interpretation and implementation differ. The rheology of PB models of foam flow is classically based on Bretherton’s flow model of a single bubble within a capillary [24], that was later extended by Hirasaki and coworkers to a train of bubbles [23] thus leading to the following expression for the effective gas viscosity:(9)where V_{g} is the interstitial gas velocity in the presence of foam. c_{f} is a constant depending on the surfactant concentration and the permeability of the porous medium [33]. The velocity power law involved in the expression of expresses the shearthinning effect of flow on foam bubbles. The exponent was determined by Bretherton for the motion of a single bubble in a capillary tube with smooth walls [24]. However, for so complex pore geometries as natural porous media, the shearthinning effect cannot be modelled with the same power law [15].
In a nutshell, PB models solve foam texture, then foam effective viscosity which depends on the latter by means of a constitutive relationship, and finally phase saturations S_{φ} and pressures p_{φ} related to the generalized multiphase Darcy framework. While PB models are certainly more consistent from a physical point of view because of the consideration of foam texture that is solved provided source terms r_{g} and r_{c} accounting for foam insitu generation and collapse in (7) are known, the reality of practice is somewhat different. Conceptually, it is not so much the relation with foam texture steadystate behaviour that is the main plus of PB models, which is well accommodated by IT models although it is not obvious [17], as the local and transient character of foam creation and destruction, driven by r_{g} and r_{c}. In this respect, PB models are good research tools. However, for the time being, when considering natural porous media it remains tricky to take advantage of this benefit of PB models because foam local and timedependent generation and collapse, modelled by source terms r_{g} and r_{c}, are intimately related to the porescale structure of the porous medium (involved in leavebehind, snapoff, and lamella division mechanisms), which again is difficult to determine accurately at large geological scales. Thus, rather paradoxically, in practice PB models may not provide more predictivity but more degrees of freedom and require solving one more (texture) equation compared to IT models. Given that PB models are equivalent to IT models at local equilibrium steady state, setting r_{g} = r_{c} in (7) [16, 17], this assumption not being intended to be general but remaining the simplest to make in practice for lack of anything better, it is legitimate to wonder about the advantage of considering PB models to simulate foambased displacements in natural porous media which are intrinsically uncertain. Notwithstanding, it is certainly not the model but the determination and the taking into account of the structure of the porous medium at small scales nested in large ones which is an issue and leads to rather conservative and basic choices in practice.
To finish with, IT models account for foam apparent viscosity that may be written(10)and depends on foam quality f_{g} = U_{g}/U, where the total velocity U = ∑U_{φ} is the sum of the phase velocities. When plotted as a function of f_{g}, has the shape of an asymmetrical bell curve that is maximum for the optimal quality in the neighborhood of which the sharpness of the transition from the low to the highquality regime is driven by the F_{2} function (6). This is a key feature for assessing foam in porous media (see [15, 16] and cited references therein for calibration examples).
In any case, given that each model has specific drawbacks and little predictive capability due to the large number of parameters that need to be adjusted beforehand, the following scaling law considerations, while developed in the context of IT models, also apply to PB models under the assumption of local equilibrium steady state.
2.2 Scaling relationships
Now we focus on the effect of permeability on foam properties. Many observations have been made by several authors, especially [34–37] in the context of oil recovery, but also in other contexts that will not be reviewed herein. Regarding soil remediation however, it is worth citing a laboratory study on a twolayer stack of porous media that evidenced the selective nature of foam [4]. In the field of enhanced oil recovery, a laboratory and phenomenological simulation study has been performed by considering simplified heterogeneous porous media, with most emphasis on chemistry however [38]. Eventually, a recent work based on laboratory measurements analysis and physical interpretation within the framework of implicittexture modelling has brought new trends for foam main drivers with respect to permeability, assuming porous media to be homothetic [15]. As the present paper is a companion of this previous published work, the reader is referred to it for further details that cannot be given here. Here are the main findings useful for this work.
For a given foam, the mobility reduction M_{ref} represents the ratio between the continuous gas mobility k_{rg}/μ_{g} and the foamedgas mobility that is defined as in IT models, or as in PB models, at optimal foaming conditions. Scaling M_{ref} for different porous media characterized by the same k_{rg} function can then be reduced to the scaling of . Analogy between flow in a capillary tube (originally formulated by Bretherton [24]) and flow in a porous medium has led most authors to use the empirical rheological law given by (9) for foamedgas.
Kovscek and Bertin [33] derived an expression of scaled with respect to rock permeability and capillary pressure by establishing an equivalence between rheological law given by (9) and Hirasaki and Lawson’s formula for the apparent viscosity of a train of bubbles of predetermined volume flowing in a capillary tube, that is [23]:(11)where R is the tube radius and R_{c} the curvature radius of the Plateau borders [18] separating gas bubbles and n_{f} the lineic foam texture, i.e. the number of lamellas per flowlength unit. This formula assumes that touching bubbles are flowing through the capillary and that the surface tension gradient effects on effective gas viscosity are negligible. Actually, that expression of the effective viscosity can be considered for scaling purpose, because it is related to the geometrical characteristics of capillary tubes, or to the properties of the equivalent porous medium represented as a bundle of such capillary tubes.
For a given foam flowing in a porous medium of permeability K and porosity Φ, capillary tube radius R is equivalent to the quantity and R_{c} is related to the disjoining pressure of the films of the foam under consideration [39]. R_{c} is therefore assumed invariant for the same foam displaced through different porous media. Furthermore, we assume that R_{c} is negligible compared to pore radii (such an assumption may not be valid however in verylowpermeability media). Therefore, scales as n_{f}R^{2} under given velocity conditions. For homothetic porous media of different permeabilities with similar porosities, R scales as and the lineic foam texture n_{f}, that is proportional to 1/R [15], scales as , hence scales as . The above analysis indicates that M_{ref} is expected to increase as the square root of permeability for homothetic porous media of different permeabilities:(12)
That scaling approximation neglects small porosity variations compared to permeability variations in many natural porous media, and is the main trend that emerges among the scaling trends reported in [15] that are reproduced in Table 1.
It is worth mentioning that aquifer permeability are commonly much higher than oil reservoir ones. Aquifer permeability can be very high, ranging from a few Darcies to hundreds of Darcies. When K is high, the scaling relationship (12) is questionable because it is based roughly on the assumption that the characteristic length of the gas bubble is of same order than pore length (one bubble per pore). This assumption is not true when permeability is high and there are several bubbles per pore.
Probing the impact of that M_{ref} scaling at the reservoir scale is the purpose of this work, with implications on the assignment of M_{ref} values to rocktypes in reservoir models.
3 Model case study setup
3.1 General considerations related to the scaling with respect to permeability
Scaling law (12) highlights the role that foam can play as a NAPL displacing fluid in heterogeneous porous media. If, as an example, we represent schematically the porous medium under consideration as a stack of layers^{1}, highpermeability layers may lead to a bypass of lowpermeability ones when a fluid (say water or gas) is injected to flush NAPL. Thus, remaining NAPL in lowpermeability layers remains trapped and may constitute a primary target for foam, not in a frontal displacement perspective but rather as foam blocking selectively the flow in more permeable layers thus generating fluid diversion in lowpermeability layers. It is precisely the diverted fluid (water and/or gas), and not necessarily foam itself, that will displace the trapped NAPL in lowpermeability zones. It is thus more a problem of conformance control than of displacement. Before moving on to the setup of the toy porous medium used to evidence the impact of foam mobility reduction scaling with respect to permeability, described in Section 3.3, let us return to the scaling law (12) to account for permeability anisotropy and raise some practical questions.
Any feasibility study implies measurements, here foam injections practiced on cores, assuming that an adequate chemistry has been developed for the conditions of pressure, temperature and brine salinity of the porous medium under consideration. As such, the question of cores selection is raised: on which core samples would it be necessary to characterize foam beforehand? Ultimately, how many measurements are needed? Obviously, scaling law (12) may guide a smart core sampling in relation to the porous medium heterogeneity characterization, rocktype per rocktype.
For the sake of demonstration, let us assume the layercake reservoir overall average horizontal permeability is mD where h_{i} is the ith layer height, is the ith layer horizontal permeability and H = ∑ h_{i} is the height of the reservoir interval. Let us further assume the interpretation of laboratory foam displacements on coreplugs yields, among other things, a reference mobility reduction of M_{ref} = 500 for that permeability 〈K_{x}〉 = 100 mD, and most of the potential issues related to foam have been solved, e.g. it does not drain/collapse too quickly and surfactant does not adsorb too much on the rock (this is not the question asked here but these are strong assumptions).
In the minimal layercake example discussed in this paper, whose reservoir structure and heterogeneity are detailed below in Section 3.3, we propose to quantify the impact of taking into account a dependence of the foam mobility reduction on the permeability, following scaling law (12), on NAPL removal. Such a scaling law practice is common for capillary pressure and (less often) for chemicals adsorption, which both^{2} scale as .
Let us then consider a layercake reservoir of fixed average horizontal permeability, say 〈K_{x}〉 = 100 mD to which is assigned a mobility reduction of M_{ref} = 500, with contrasts in permeability from one layer to another. This contrast is varied from one realization to another in order to probe the foam flow behaviour on independent configurations whose heterogeneity degree varies gradually, from very large to moderate deviations to the average permeability (preserving the mean), and illustrate the implications of M_{ref} scaling on the assignment of M_{ref} values to rocktypes in reservoir models. In each case, M_{ref} is associated with K_{ref} and the mobility reduction is computed for each cell of the reservoir mesh following (12) as:(13)
According to this relationship, fixing K_{ref} and varying the assigned M_{ref} is equivalent to fixing M_{ref} and varying the associated K_{ref}. Thus, varying the M_{ref} value associated with K_{ref} = 100 mD, say either M_{ref} = 500 or 1581, is equivalent to assigning M_{ref} = 500 to either K_{ref} = 100 or 10 mD, respectively. These two examples will be considered hereafter.
Let us now consider how to manage the directional character of the flow in relation to the permeability tensor K. In the simplest layercake geometry where the layers are all flat and parallel (they can however be tilted), the permeability tensor is diagonal in the Cartesian coordinate system (x, y, z) such that the x and y coordinates span the geological layers, thus K =d iag(K_{x}, K_{y}, K_{z}) and (13) reads for the scaled mobility reduction:(14)
This relation writes identically if the layers are not flat but parabolic (as in the case of an anticline structure) or of any shape in relation to a deposit structure, by substituting to the Cartesian coordinate system (x, y, z) the curvilinear one that spans the layers. Doing so, the mobility function FM given by (5) also acquires a directional character upon M_{ref} scaling, and it will now be denoted as follows:(15)
As the vertical permeability K_{z} usually differs from the horizontal one K_{x} (assuming K_{x} = K_{y} for the considered simple layercake geometry), which is driven by the anisotropy ratio K_{z}/K_{x}, so do the reference permeabilities K_{ref,z} and K_{ref,x}, and the mobility reductions M_{ref,z} and M_{ref,x}, FM_{z} and FM_{x} in the corresponding directions.
To finish with, we do not consider chemicals adsorption nor the impact of NAPL on foam in this work since they both primarily impact economics and chemicals slug size, which is not the topic addressed here. Of course, they are of critical interest for technical and economical evaluation (the reader is referred to [20] for such an evaluation for naturally fractured oilwet reservoirs, e.g.). It is also assumed that the reservoir can withstand the imposed injection flow rate without suffering from any injectivity issues due to pressure buildup. This is in practice one of the key points to be addressed for any viscous fluid injection in lowpermeability formations, but since we assume the existence of very permeable layers as detailed in Section 3.3, this assumption remains reasonable. It is worth noting that the injectivity of shearthinning (or shearthickening) fluids whose viscosity varies a lot for large velocity gradients in the nearwellbore, as it can be the case of preformed foam or polymer solution, requires a specific treatment of injectivity and productivity indices [40] whose traditional formulation presupposes a constant viscosity over the well gridblock [41], which underestimates the injectivity of shearthinning fluids.
Commingled foam injection is considered first as a baseline prior to any selective injection: gas and surfactant solution are injected over the full 50 m height perforated reservoir interval (see Sect. 3.3), not in a specific layer. The permeability tensor with components K_{x}(x, z) and K_{z}(x, z) is considered either isotropic or strongly anisotropic, i.e. K_{z}/K_{x }= 1 or , respectively, where K_{z}(x, z) is deduced from K_{x}(x, z) which is geostatistically generated by imposing a constant mean 〈K_{x}〉 = 100 mD, as explained in Section 3.3. Two types of injection are considered: either a coinjection of surfactant solution and gas (with a quality of f_{g} = 0.8) which corresponds to a preformed foam injection, and surfactant solution alternating gas injection for which foam is formed insitu, as reported in Figure 1, with an imposed constant total injection flow rate of 54.8 m^{3}/day for which one pore volume (PV) of the reservoir is injected within one year. Indeed, since foams are not a proven largescale recovery process, the mode of fluids injection is also being investigated, with the literature reporting roughly as many coinjection tests as alternate injection. However, preformed foam injection may be preferred for fractured reservoirs with poor insitu foam generation mechanisms in the fracture network. In both cases, 0.25 PV of a 5000 ppm concentrated surfactant solution is injected. All the simulations are performed with IFP Energies nouvelles reservoir simulator Puma, which accounts for the abovementioned foam features, including the M_{ref} scaling law relationships (14) and (15) (see [15] and cited references therein).
Fig. 1 Timeline of surfactant solution/gas coinjection (blue) and surfactant solution alternating gas injection (red) scenarios, where W, G and W* denote water, gas and surfactant solution, respectively. Timescale is indicated in injected pore volumes (PV, volume of fluid injected over the pore volume of the porous medium). Both cases have a prefoam water or gas injection history of 1 PV. In each case the volume of surfactant solution injected of 0.25 PV is identical: in the coinjection case, it is distributed over a total volume of gas and surfactant solution of 1.25 PV according to the imposed quality of f_{g} = 0.8 that corresponds to optimal foaming conditions. There is no such quality control for the surfactant solution alternating gas injection, where foam is formed by fingering of gas (injected after surfactant solution) through surfactant solution, or at least mixing of the two, due to the different mobility of the fluids. 
3.2 Rockfluid properties
The reservoir under consideration has an average temperature of T_{res} = 100 °C and an average pressure of p_{res} = 100 bar. Rockfluid properties are representative of Berea sandstone with a foam made up of a surfactant aqueous solution and a methane/CO_{2} mixture (62% CO_{2}, 38% CH_{4}) with a quality f_{g} = 0.8. The reader is referred to the previous work [42] for further details. These conditions are far from nearsubsurface aquifer ones where the temperature is close to 15 °C and the pressure is about a few bar. In the same vein, the gas that will be used to generate foam for soil remediation will be nitrogen or simply air. However, these points have no bearing on the work presented, which remains phenomenological in nature.
Threephase flow is modelled as a combination of the underlying water/NAPL and gas/NAPL twophase subsystems following Stone I model [21, 22]. Water is considered as the wetting phase while gas is the nonwetting phase, NAPL being the common and intermediate phase to both subsystems. Each twophase subsystem is set up symmetrically with respect to its wetting phase following BrooksCorey relationships [43, 44], assuming powerlaw relative permeabilities and capillary pressures. Relative permeabilities to water and NAPL, k_{rw} and k_{row}, and capillary pressure p_{cow} = p_{o} − p_{w}, depend on the normalized mobile water saturation S_{W} = (S_{w} − S_{wi})/(1 − S_{orw} − S_{wi}), and read:(16)where S_{wi} is the irreducible water saturation, S_{orw} is the residual NAPL saturation with respect to water, κ_{w} and κ_{ow} are the maximum relative permeabilities to water and NAPL, and π_{ow} is the entry capillary pressure. The gas/NAPL subsystem, for which NAPL is the wetting phase, is set up symmetrically as a function of the normalized mobile gas saturation S_{G} = (S_{g} − S_{gr})/(1 − S_{org} − S_{wi} − S_{gr}):(17)with p_{cog} = p_{g} − p_{o} and where S_{gr} is the irreducible gas saturation and S_{org} is the residual NAPL saturation with respect to gas. Figure 2 shows flow curves (16) and (17) with parameters given in Table 2.
Fig. 2 Water/NAPL and gas/NAPL twophase subsystems powerlaw relative permeabilities and capillary pressures (16) and (17) with parameters given in Table 2. 
To finish with, the specific choice of Stone’s model is arbitrary and only justified by its widespread use in reservoir engineering. This “choice of the majority” puts a restriction on the work carried out and would merit further specific investigations related to foam flow, disentangling properly transport and adsorption, that are out of the scope of the present work.
3.3 Model layered porous medium setup
The reservoir under consideration is a twodimensional crosssectional xz layered porous medium of 200 m length, 50 m height and 10 m width^{3}. Porosity is set constant to Φ = 0.2. The average horizontal permeability in the xcoordinate is set constant to 〈K_{x}〉 = 100 mD while the permeability anisotropy ratio K_{z}/K_{x}, from which the vertical permeability K_{z} is deduced, is set to either (strongly anisotropic) or 1 (isotropic). Several permeability map realizations K_{x}(x, z) of increasing degree of heterogeneity are considered. They range from quasihomogeneous to highly heterogeneous by varying the horizontal permeability standard deviation , 50, 100, 150, 200, 250 mD, as reported in Figures 3 and 4, using a fast Fourier transform moving average algorithm [45, 46].
Fig. 3 Comparison of several twodimensional horizontal permeability field K_{x}(x, z) realizations (approximately lognormal distributed) with constant mean 〈K_{x}〉 = 100 mD and standard deviation , 50, 100, 150, 200, 250 mD, or equivalently, DykstraParsons coefficient V_{DP} = 0.22, 0.37, 0.57, 0.66, 0.72, 0.75. Typical layers thickness is driven by and close to the correlation length L_{z} = 10 m. 
Fig. 4 Comparison of several twodimensional horizontal permeability field K_{x}(x, z) realizations (approximately lognormal distributed) with constant mean 〈K_{x}〉 = 100 mD and standard deviation , 50, 100, 150, 200, 250 mD, or equivalently, DykstraParsons coefficient V_{DP} =0.22, 0.37, 0.57, 0.66, 0.72, 0.75. Typical layers thickness is driven by and close to the correlation length L_{z} = 5 m. 
In each case, K_{x}(x, z) is approximately lognormal distributed with a constant mean 〈K_{x}〉 = 100 mD. Layers typical thickness over which K_{x}(x, z) is approximately constant is controlled by a Gaussian autocorrelation function whose correlation length L_{z} in the zcoordinate is set to 10 or 5 m, while the correlation length in the xcoordinate is arbitrarily large (L_{x} = 1000 m) in order to render a layercake with excellent lateral continuity. Under these conditions, the characteristic layers thickness is close to the correlation length L_{z}, and we will abuse the language by sometimes referring to L_{z} as the typical layers thickness.
The degree of heterogeneity in permeability may be quantified by the DykstraParsons coefficient [47] defined as V_{DP }= (e^{ μ } − e^{ μ } − ^{ σ })/e^{ μ } = 1 − e^{−σ}) for a lognormal random variable K of parameters μ = 〈ln K〉 and σ^{2} = 〈(lnK − μ)^{2}〉. V_{DP} writes in terms of the permeability mean 〈K〉 and variance as^{4}:(18)where the subscript x has been dropped for the sake of readability. Thus, V_{DP} → 0 for a homogeneous permeability distribution such that σ_{K}/〈K〉 → 0 whereas V_{DP} → 1 for an extremely heterogeneous distribution with σ_{K}/〈K〉 → ∞. Abovementioned permeability standard deviations of , 50, 100, 150, 200, 250 mD are equivalent to the DykstraParsons coefficients V_{DP} = 0.22, 0.37, 0.57, 0.66, 0.72, 0.75, as indicated in Figures 3 and 4. The mesh is discretized in x and zcoordinate with 100 and 50 cells, respectively; cells are therefore of constant size of 2 × 1 m thus the mesh is really very refined. Injection and production wells are located at the boundaries x = 0 and 200 m respectively, and are both perforated over all the layers for all z ∈ [0, H]. Noflow conditions are assigned to the top and the bottom of the layercake z = 0 and H.
To finish with, it is worth mentioning the question of flow direction with respect to strata in a stratified system (which is already a limitation, since, in the case of NAPL problems, many shallow aquifers are not necessarily simple layered systems) is rather tricky. This work is restricted to the minimal complexity of a stratified bedding in order to quantify the coupling of foam viscous flow with heterogeneity. In practice, many configurations can arise, ranging from a very thick homogeneous shallow aquifertype reservoir interval to heterogeneous puzzletype geometries, with well architectures of varying complexity, that are not addressed in the present work.
4 Numerical simulations
We now compare foam flow performance on NAPL removal, depending on whether foam is preformed under coinjection of surfactant solution and gas, or forms insitu under surfactant solution alternating gas injection, according to the injection scenarios given in Figure 1. The foam maximum mobility reduction M_{ref} is scaled in relation to the porous medium permeability distribution following (14) and (15), or not. Both injection scenarios (coinjection or alternate injection) have a prefoam water or gas injection history of one injected pore volume in order to diagnose flow performance under water or gas flooding, to identify the reservoir unswept zones and to anticipate the contribution that foam would have. In each case the 0.25 PV volume of surfactant solution injected is identical: in the coinjection case, it is distributed over a total volume of gas and surfactant solution of 1.25 PV according to the imposed quality of f_{g} = 0.8 that corresponds to optimal foaming conditions. There is no such quality control for the surfactant solution alternating gas injection, where foam is formed by fingering of gas (injected after surfactant solution) through surfactant solution, or at least mixing of the two, due to the different mobility of the fluids. The twodimensional xz layercake porous medium under consideration is initially saturated with NAPL up to the irreducible water saturation^{5} S_{wi} given in Table 2. We proceed as follows:

A mobility reduction of either M_{ref,x} = 500 or 1581 is assigned to the reference horizontal permeability K_{ref,x} = 100 mD that is equal to the reservoir average horizontal permeability 〈K_{x}〉. As explained in Section 3.1, this is equivalent to assigning M_{ref,x} = 500 to either K_{ref,x} = 100 or 10 mD, respectively. This will give two configurations for further comparisons.

The permeability tensor with components K_{x}(x, z) and K_{z}(x, z), is considered either isotropic or strongly anisotropic setting either K_{z}/K_{x }= 1 or , where K_{z}(x, z) is deduced from K_{x}(x, z) which is geostatistically generated by imposing a constant mean 〈K_{x}〉 = 100 mD, as explained in Section 3.3.

Several twodimensional horizontal permeability field K_{x}(x, z) realizations (approximately lognormal distributed) are considered by imposing the mean 〈K_{x}〉 = 100 mD and varying the standard deviation , 50, 100, 150, 200, 250 mD, or equivalently, the DykstraParsons coefficient V_{DP} = 0.22, 0.37, 0.57, 0.66, 0.72, 0.75, as shown in Figures 3 and 4.

Typical layers thickness driven by and close to the correlation length L_{z} of either 10 or 5 m are considered.
To start with, insitu foam flows without and with M_{ref} scaling are compared in the two following Sections 4.2 and 4.3 for the largest DykstraParsons coefficient V_{DP} = 0.75 for foam flooding obtained under coinjection of surfactant solution and gas or surfactant solution alternating gas injection. The initial insitu behaviour under water or gas flooding is first analyzed in Section 4.1. Figures that compare the insitu dynamic behaviour of saturations S_{w}, S_{g}, S_{o}, foaming agent concentration C_{f}, mobility reduction function in the xcoordinate (FM_{x})^{−1}, as well as the horizontally and vertically averaged corresponding profiles, have been moved in Appendix A.2 for the sake of readability.
Production data such as NAPL recovery, produced chemicals and gas including chemicals breakthrough time, are compared in Sections 5. Eventually, flow performance with and without M_{ref} scaling are compared in Section 6 and characterized by sweep trends in relation with the DykstraParsons coefficient.
4.1 Flow performance baseline without foam
To start with, let us analyze the flow without chemical additives, i.e. without foam, considering the two injection scenarios given in Figure 1, coinjection of water and gas or water alternating gas injection, without surfactant in the water phase. This baseline will serve as a comparison point for simulated foam flows with and without M_{ref} scaling, which will also be compared with each other.
Very classically, we find the following wellknown result: all other things being equal, the higher the reservoir degree of heterogeneity that is transcribed by V_{DP}, the lower the vertical sweep of the reservoir by water and/or gas flooding. This can be seen very clearly on NAPL production curves shown in (top) Figures 5a and 8a: as V_{DP} increases from 0.22 to 0.75 (see Figs. 3 and 4), the final NAPL recovery, expressed as the ratio of the cumulative volume of NAPL recovered over the initial volume of NAPL in place, drops from 70% to 30%, approximately. In the most heterogeneous case V_{DP} = 0.75, only the most permeable fraction of the reservoir is swept as shown in Figure A1. Once the injected displacing fluid has broken through to the producing well, a large fraction of the reservoir is hydraulically bypassed and will mostly remain so, as evidenced by the very slow NAPL recovery latetime evolution reported in Figures 5a and 8a which reflects the latetime microscopic displacement efficiency, due to the propagation in the most permeable layers of the very loose tail of the displacing saturation front that relaxes very slowly and is very inefficient in NAPL mobilization. Put more simply: once the displacement front has broken through, not much happens, or at least very slowly.
Fig. 5 Coinjection case with M_{ref,x} = 500, K_{ref,x }= 100 mD, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Comparison of (a) NAPL recovery factor (RF, produced NAPL volume over initial NAPL volume in place, OOIP) and (b) produced surfactant concentration (C_{f}) and cumulative mass (M_{f}) time evolution as a function of V_{DP} for chemicalsfree injection (“nofoam”), foam with constant M_{ref,x} (“ref”, dashed lines) and scaled M_{ref,x} (“mod”, solid lines). Relative variations are either in relation to a chemicalsfree twin injection (“/ nofoam”) or a constant M_{ref,x} case (“/ foam”). 
Fig. 6 Coinjection case with K_{ref,x} = 100 mD, (a) M_{ref,x} = 500, (b) M_{ref,x} = 1581 (or equivalently M_{ref,x} = 500, (a) K_{ref,x} = 100 mD, (b) K_{ref,x} = 10 mD), K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Comparison of produced cumulative gas volume (V_{g}) time evolution as a function of V_{DP} for chemicalsfree injection (“nofoam”), foam with constant M_{ref,x} (“ref”, dashed lines) and scaled M_{ref,x} (“mod”, solid lines). Relative variations are either in relation to a chemicalsfree twin injection (“/ nofoam”) or a constant M_{ref,x} case (“/ foam”). 
Fig. 7 Coinjection case with M_{ref,x} = 1581, K_{ref} = 100 mD (or equivalently M_{ref,x} = 500, K_{ref,x} = 10 mD), K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Comparison of (a) NAPL recovery factor (RF, produced NAPL volume over initial NAPL volume in place, OOIP) and (b) produced surfactant concentration (C_{f}) and cumulative mass (M_{f}) time evolution as a function of V_{DP} for chemicalsfree injection (“nofoam”), foam with constant M_{ref,x} (“ref”, dashed lines) and scaled M_{ref,x} (“mod”, solid lines). Relative variations are either in relation to a chemicalsfree twin injection (“/ nofoam”) or a constant M_{ref,x} case (“/ foam”). 
Fig. 8 Alternate injection case with M_{ref,x} = 500, K_{ref,x} = 100 mD, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Comparison of (a) NAPL recovery factor (RF, produced NAPL volume over initial NAPL volume in place, OOIP) and (b) produced surfactant concentration (C_{f}) and cumulative mass (M_{f}) time evolution as a function of V_{DP} for chemicalsfree injection (“nofoam”), foam with constant M_{ref,x} (“ref”, dashed lines) and scaled M_{ref,x} (“mod”, solid lines). Relative variations are either in relation to a chemicalsfree twin injection (“/ nofoam”) or a constant M_{ref,x} case (“/ foam”). 
Fig. 9 Alternate injection case with M_{ref,x} = 1581, K_{ref,x}= 100 mD (or equivalently M_{ref,x} = 500, K_{ref,x} = 10 mD), K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Comparison of NAPL recovery factor (RF, produced NAPL volume over initial NAPL volume in place, OOIP) time evolution as a function of V_{DP} for chemicalsfree injection (“nofoam”), foam with constant M_{ref,x} (“ref”, dashed lines) and scaled M_{ref,x} (“mod”, solid lines). Relative variations are either in relation to a chemicalsfree twin injection (“/ nofoam”) or a constant M_{ref,x} case (“/ foam”). 
Fig. 10 Comparison of the maximum NAPL recovery factor (RF, produced NAPL volume over initial NAPL volume in place, OOIP) relative variation under foam flooding with scaled M_{ref,x} relatively to constant M_{ref}, as a function of the injection sequence (coinjection or alternate injection), V_{DP}, K_{ref,x} = 100, 10 mD assigned to M_{ref,x} = 500 (or equivalently M_{ref,x} = 500, 1581 assigned to K_{ref,x} = 100 mD), K_{z}/K_{x} and L_{z}. 
Fig. 11 Comparison of the (a) final produced surfactant cumulative mass (M_{f}) and (b) chemicals breakthrough time (BT) relative variation under foam flooding with scaled M_{ref,x} relatively to constant M_{ref,x}, as a function of the injection sequence (coinjection or alternate injection), V_{DP}, K_{ref,x} = 100, 10 mD assigned to M_{ref,x} = 500 (or equivalently M_{ref,x} = 500, 1581 assigned to K_{ref,x} = 100 mD), K_{z}/K_{x} and L_{z}. 
Fig. 12 Comparison of the produced cumulative gas volume (V_{g}) relative variation under foam flooding with scaled M_{ref,x} relatively constant M_{ref,x}, as a function of the injection sequence (coinjection or alternate injection), V_{DP}, K_{ref,x} = 100, 10 mD assigned to M_{ref,x} = 500 (or equivalently M_{ref,x} = 500, 1581 assigned to K_{ref,x} = 100 mD), K_{z}/K_{x} and L_{z}. 
Fig. 13 Comparison of the differential behaviour of production results under foam flooding with scaled M_{ref,x} relatively to constant M_{ref,x} as a function of the injection sequence (coinjection or alternate injection), V_{DP} = 0.22, 0.37, 0.57, K_{ref,x} = 100, 10 mD assigned to M_{ref,x} = 500 (or equivalently M_{ref,x} = 500, 1581 assigned to K_{ref,x} = 100 mD), K_{z}/K_{x} and L_{z} (subscripts “ani” and “iso” stand for K_{z}/K_{x} = 0.1 and 1 while “L_{z}” and “L_{z}/2” correspond to L_{z} = 10 and 5 m, respectively). 
Fig. 14 Comparison of the differential behaviour of production results under foam flooding with scaled M_{ref,x} relatively to constant M_{ref,x} as a function of the injection sequence (coinjection or alternate injection), V_{DP} = 0.66, 0.72, 0.75, K_{ref,x} =100, 10 mD assigned to M_{ref,x} = 500 (or equivalently M_{ref,x} = 500, 1581 assigned to K_{ref,x} = 100 mD), K_{z}/K_{x} and L_{z} (subscripts “ani” and “iso” stand for K_{z}/K_{x} = 0.1 and 1 while “L_{z}” and “L_{z}/2” correspond to L_{z} = 10 and 5 m, respectively). 
We will not go into further detail on this wellknown type of behaviour which has been the subject of much work [48–50] and it should be retained, since the sweep efficiency generally writes as the product of the macroscopic (here vertical) and microscopic sweep efficiencies, how crucial it is to improve the volumetric sweep (vertical and/or areal) when addressing very heterogeneous porous media, because otherwise whatever the improvement in the microscopic displacement efficiency (by increasing the displacing fluid viscosity and/or modifying the interfacial rockfluids capillary forces), it may concern only a tiny volumetric fraction of the porous medium under consideration. For example, Figures 5a and 8a clearly illustrate, for V_{DP} = 0.22 (very moderate heterogeneity), the impact of the displacing fluid viscosity, depending on whether it is water or gas, on the transient earlytime NAPL recovery: water indeed presents a much larger displacement saturation front than gas (assumed here to be immiscible with NAPL), mostly due to the difference in viscosity, thus the transient earlytime NAPL recovery is more efficient and much faster for water than for gas flooding. It will also take much longer for gas than for water flooding to reach a pseudo steady state such that the tail of the displacing saturation front has completely relaxed, as saturation levels of NAPL to be displaced by the tail at the back of the displacing saturation front are larger for gas flooding, with the practical difference that in general gas injects more easily than water, again due to its viscosity (however, this is not the case for the simulations presented which were carried out at an imposed flow rate).
It is worth noting that in moderate to low V_{DP} configurations, the vertical sweep is much better because more and more layers are accessible to the flow and are swept, which is enhanced by the alternating injection of a dense fluid (water) and a less dense fluid (gas), known as WAG for water alternating gas, in order to compensate for gravity segregation and poor vertical sweep that can follow independently of heterogeneity. Coinjection of water and gas, which is equivalent to the limiting case of infinitesimal alternating water and gas slugs, is obviously much more efficient in NAPL recovery kinetics as can be seen by comparing Figures 5a and 8a, because of fluids gravity segregation which is permanently compensated, unlike the water alternating gas injection, the NAPL recovery of which is done in two stages, as shown in Figure 8a, which correspond to the propagation of gas then water. In the very heterogeneous cases, the mitigation of gravity segregation by coinjecting or alternating water and gas is much less efficient as heterogeneity rules. This trend, which however is not a general rule and depends on the height of the most permeable layers and interwell spacing, is very clear in Figure 8a as V_{DP} increases.
This is a quite classic and predictable behaviour that foam is precisely expected to remedy. We will indeed show that the more heterogeneous the reservoir is (the higher V_{DP} is), the smaller the vertical sweep of the reservoir by water or gas is (only the most permeable layers are swept), but the more efficient foam is and the better the incremental vertical sweep and NAPL recovery are compared to a twin injection type without foam (without surfactant). Let us see why and how.
4.2 Comparison of insitu foam flow with and without M_{ref} scaling: coinjection
We now focus on foam injection in the form of a 1.25 PV coinjection of surfactant solution and gas that follows a 1 PV water injection, and continues with a 1.75 PV water injection again, as indicated in Figure 1. Due to the imposed quality f_{g} = 0.8, the coinjection of surfactant solution and gas consists of a 0.25 PV surfactant solution slug. The behaviour of the prefoam water injection has been commented in the previous section. Let us see now what the foam brings.
Figures A1, A2 and A3 compare the insitu dynamic bevaviour of foam flow (i) after 1 PV of water injection prior to surfactant injection, (ii) after 0.25 PV of surfactant injection (end of the surfactant injection), and (iii) at the end of the full injection sequence, without and with M_{ref} scaling, for V_{DP} = 0.75, M_{ref} = 500, K_{ref,x} = 100 mD, K_{z}/K_{x} = 0.1 and L_{z} = 10 m.
After 1 PV of water injection, the most permeable central layer has been well swept, as can be seen from water and NAPL saturation profiles. This layer clearly dominates the flow so that the rest of the reservoir is little swept, as evidenced in Figure A1.
After 0.25 PV of surfactant solution injection, foam has propagated in the most permeable layers where there is preferential flow as shown by the mobility reduction function component (FM_{x})^{−1} in Figure A2 that is much larger in highpermeability layers when scaling M_{ref} compared to a constant M_{ref}, and, on the contrary, that is smaller in the less permeable layers. With a constant M_{ref}, the surfactant concentration profile mostly invades the zones that were previously swept by water, that is mainly the most permeable layers and a small fraction of the less permeable peripheral layers. An improvement in the sweeping of the latter can be noted, which did not occur after the 1 PV water injection and which would certainly not have occured if the water injection had been continued. In addition, the dominant drain is better swept microscopically compared to water flooding due to the sharp increase in the displacing foamed gas saturation front, as can be seen from the NAPL saturation profile, and which is the first expected effect according to a classic BuckleyLeverett calculation [14, 49–51].
Alternatively, it can be noted qualitatively from the saturation profiles that the insitu quality of the foam is very different from the quality imposed at the injection well of f_{g} = 0.8. While this is due to the coupling of the differential behaviour of the phases mobilities and reservoir heterogeneity (and also capillary trapping), it may raise a question about the conclusions that can be drawn from laboratory evaluations by corescale displacements, where the quality is often imposed and constant, all the more so at steady state. At least, the inevitable variation of f_{g} in the reservoir implies paying special attention to the behaviour of the foam apparent viscosity as a function of f_{g}, and it would obviously be desirable that the latter varies little in the neighborhood of the optimal foam quality . In any case this is monitorable by simulation provided that a reliable and accurate characterization of the rockfluids system is available, and that onsite pilot tests are suitably instrumented for monitoring purposes.
When scaling M_{ref}, (FM_{x})^{−1} increases in the most permeable layer while it decreases in the lowpermeable ones, compared to a constant M_{ref}, as expected from (14) and (15). This leads to a more significant invasion of the peripheral layers as can be seen from the surfactant concentration and water, gas and NAPL saturation profiles in Figure A2.
At the end of the 0.25 PV surfactant injection, the profiles shown in Figure A3 are very different when comparing a constant M_{ref} to a scaled one. One can observe an almost complete desaturation of the central drain in both cases, but above all a much better sweep of the peripheral layers due to a redistribution of the flow from the central drain when M_{ref} is scaled. The differences in surfactant concentration and saturation profiles are very significant in the peripheral layers and are not primarily due to differences in (FM_{x})^{−1} values in the peripheral layers, which are roughly speaking close there whether M_{ref} is scaled or not.
Eventually, at the end of the total injection sequence shown in Figure A3, surfactant has been evacuated from the central drain (which again behaves like a bypass) and remains further trapped in the peripheral layers when M_{ref} is scaled, which is consistent with the previously observed flow redistribution.
A similar behaviour is observed when considering M_{ref,x} =1581 and K_{ref,x} = 100 mD (or equivalently M_{ref,x }= 500 and K_{ref,x} = 10 mD), as reported in Figures A4 and A5. The differences are more pronounced and can be explained by significantly larger values of M_{x} and (FM_{x})^{−1}.
Foam flow for K_{ref} = 100 mD and K_{z}/K_{x} = 1, which is shown in Figure A6, also behaves similarly with more interlayer transfers due to the isotropic nature of the permeability which globally enhances the vertical sweep.
A last comparison is performed in Figure A7 where the impact of the typical layers thickness, L_{z} = 10 or 5 m, is shown for a scaled M_{ref} with M_{ref,x} = 500, K_{ref,x} = 100 mD and K_{z}/K_{x} = 0.1. One can note a slightly better sweep for thinner layers such that L_{z} = 5 m which promote interlayer vertical transfers.
Although we do not report insitu profiles for every V_{DP} values, it will be shown in Section 5 where production results are compared that the differences which have been observed for V_{DP} = 0.75 gradually fade as V_{DP} decreases. We will show further in Section 6 from production data that a V_{DP} threshold exists from which this sweep improvement manifests itself significantly when scaling M_{ref}.
4.3 Comparison of insitu foam flow with and without M_{ref} scaling: alternate injection
Let us now move on to foam flow where foam is formed insitu by alternating surfactant solution and gas slugs, as indicated in Figure 1. Unlike the coinjection case described in the previous section, there is no quality control for the surfactant solution alternating gas injection as foam is formed by fingering of gas (injected after surfactant solution) through surfactant solution, or at least mixing of the two, due to their different mobility. The injected surfactant solution volume of 0.25 PV is the same as for the coinjection case. The behaviour of the prefoam gas injection has been already commented in Section 4.1.
Figures A9, A10 and A11 compare the insitu dynamic behaviour of foam flow (i) after 0.25 PV of surfactant solution injection (end of the surfactant solution injection), (ii) after 0.9 PV of gas injection (post surfactant solution injection), and (iii) after 1.8 PV of gas injection (post surfactant solution injection), without and with M_{ref} scaling, for V_{DP} = 0.75, M_{ref,x} = 500, K_{ref,x} = 100 mD, K_{z}/K_{x} = 0.1 and L_{z} = 10 m.
There is a similar behaviour to coinjection with a very significant sweep improvement when scaling M_{ref}, but with an additional period of time for the foam to form and lead to the same effects, which is expected for alternate injection and obviously depends on its design in terms of slug size. This can be seen very clearly by comparing Figures A2 and A9: while after 0.25 PV of surfactant solution injected the fluids diversion due to foam is clear for coinjection, it is hardly observed for the alternate injection. It is indeed necessary to wait until the following gas slug gives rise to the formation of foam. Regarding chemical additives, the main difference with coinjection is the less surfactant trapping in lowpermeability layers. It should be noted that the remaining surfactant that can be observed in the central dominant layer in Figure A11 is not significant at all since it corresponds to very a very low water saturation tail that results from the poor displacement efficiency of gas injection. We will show further on by inspecting production results the differences that exist between alternate injection and coinjection, as there are. Let us quantify all this now in terms of the fluids produced by varying the degree of heterogeneity transcribed by V_{DP}.
5 Comparison of foam flow production results with and without M_{ref} scaling
We now explore the behaviour of the production results time evolution as a function of the DykstraParsons coefficient over the range V_{DP} = 0.22, …,0.75. Again, we compare the relative impact of the scaling of M_{ref} with respect to simulations where M_{ref} is constant. The analyzed production results are (i) NAPL recovery, expressed as the ratio of the cumulative volume of NAPL recovered over the initial volume of NAPL in place, (ii) produced surfactant concentration and cumulative mass, and (iii) produced cumulative gas volume.
The relative variations in these quantities are quantified in relation to (i) the twin chemicalsfree injection type when scaling M_{ref} or not (coinjection with chemicals is compared to coinjection without chemicals, and so on for alternate injection), and (ii) to a foam injection with a constant M_{ref} when considering a scaled M_{ref} (coinjection with M_{ref} scaling is compared to coinjection without M_{ref} scaling, and so on for alternate injection).
Although all the configurations reported in Table 3 were analyzed, it is sufficient to comment only the results obtained for K_{ref,x} = 100 mD and either M_{ref,x} = 500 or 1581 (or equivalently M_{ref,x} = 500 and either K_{ref,x} = 100 or 10 mD) and K_{z}/K_{x} = 0.1 within a coinjection (Sect. 5.1) or alternate injection (Sect. 5.2), for the sake of conciseness and readability. All the results of Table 3 configurations are presented in a more synthetic manner in Section 6.
5.1 Comparison of foam flow production results with and without M_{ref} scaling: coinjection
To begin with, we note the extreme sensitivity of the NAPL recovery of a chemicalsfree injection (in absence of foam) to the DykstraParsons coefficient, as reported in Figure 5a: as V_{DP} increases from 0.22 to 0.75, NAPL recovery approximately drops from 70% to 30%. The more heterogeneous the reservoir is (the higher V_{DP} is), the smaller the vertical sweep of the reservoir by water or gas is (only the most permeable layers are swept). This behaviour is well known and has been studied in detail by numerous authors for water flooding [48–50].
Under foam flooding, results differences whether M_{ref} is scaled or fixed are all the higher as the degree of heterogeneity, quantified byV_{DP}, is higher. Thus the more heterogeneous the reservoir is, the smaller the vertical sweep of the reservoir by water or gas is, but the more efficient the foam is in its diversion role and the better the incremental vertical sweep and NAPL recovery are compared to a twin injection type without foam. The differences in NAPL recovery can reach about 20 and 10% during the transient regime, when scaling M_{ref} with M_{ref,x} = 500, K_{ref,x} = 100 mD, for V_{DP} = 0.72 and 0.75, respectively and relatively to a foam injection with a constant M_{ref}, which is very significant. This is due to a better vertical sweep as explained in Section 4.2.
Also, we note a very strong impact on the surfactant breakthrough time and on the mass of surfactant recovered for V_{DP} ≥ 0.57 when M_{ref} is scaled with M_{ref,x} = 500 and K_{ref,x} = 100 mD, as reported in Figure 5b. The higher V_{DP} is, the earlier the chemicals breakthrough and the larger the amount of recovered chemicals are for a constant M_{ref}. This trend remains when scaling M_{ref} but there is a significant delay in chemicals breakthrough and a reduction in recovered chemicals of about 20–40% for large V_{DP} values (V_{DP} ≥ 0.57) over constant M_{ref} cases. Breakthrough time relative differences are quantified in Section 6. Again, this is due to a better vertical sweep that leads to additional trapping of surfactant in lowpermeability layers, as explained in Section 4.2.
To finish with, the produced volume of gas is sensitive to V_{DP} but much less than NAPL and surfactant recovery, as reported in Figure 6a. It is lowered by 5% for V_{DP} = 0.75 when scaling M_{ref} compared to a constant M_{ref}. Such less contrasted behaviour is due to the high mobility of gas in the absence of foam generally speaking.
The behaviours observed for M_{ref,x} = 1581 and K_{ref,x} = 100 mD (or equivalently M_{ref,x} = 500 and K_{ref,x} = 10 mD) are similar and more pronounced when M_{ref} is scaled compared to a constant M_{ref}, as reported in Figures 6b, 7a and 7b, and as already observed from the insitu profiles in Section 4.2. Let us give some orders of magnitude: when scaling M_{ref}, transient NAPL recovery is improved by about 35% and 20% for V_{DP} = 0.75 and 0.72, respectively, while chemicals production is lowered by about 50–60% for V_{DP} ≥ 0.57. Gas production does not exhibit significant differences compared to the M_{ref,x} = 500 mD and K_{ref,x} = 100 mD case.
5.2 Comparison of foam flow production results with and without M_{ref} scaling: Alternate injection
Let us come back the surfactant solution alternating gas injection. The same production results time evolution as previously is shown in Figures 8a and 8b for M_{ref,x} = 500 and K_{ref,x} = 100 mD, and in Figure 9 for M_{ref,x} = 1581 and K_{ref,x} = 100 mD (or equivalently M_{ref,x} = 500 and K_{ref,x} = 10 mD).
A critical dependence of NAPL recovery on V_{DP} for a chemicalfree injection is again noted: as V_{DP} increases from 0.22 to 0.75, NAPL recovery drops from 70% to 20%, approximately. By scaling M_{ref} it is improved by about 25%, 10% and 7% relatively to a constant M_{ref} for V_{DP} = 0.75, 0.72 and 0.66 when M_{ref,x} = 500 and K_{ref} = 100 mD, and by 40%, 23% and 10% relatively to a constant M_{ref} for the same V_{DP} values when M_{ref,x} = 1581 and K_{ref} = 100 mD. Again, the more heterogeneous the reservoir is (the higher V_{DP} is), the smaller the vertical sweep of the reservoir by water or gas is, but the more efficient the foam is in its diversion role and the better the incremental vertical sweep and NAPL recovery are compared to a twin injection type without surfactant.
6 Sensitivity study on foam flow performance differential behaviour with and without M_{ref} scaling
We now give a more global and synthetic view of the differentiated behaviour of foam flow performance by covering all the configurations reported in Table 3, 96 in number. To sum up:

Foam flow is considered through either coinjection of surfactant solution and gas or surfactant solution alternating gas injection. In both cases 0.25 PV of a 5000 ppm concentrated surfactant solution is injected as indicated in Figure 1. All foam flow simulations are performed at a constant total injection flow rate such that 1 PV of the reservoir is injected within one year.

A mobility reduction of M_{ref,x} = 500 or 1581 is assigned to the reference horizontal permeability K_{ref,x} = 100 mD that fits the average horizontal permeability 〈K_{x}〉. The case M_{ref,x} = 1581 and K_{ref,x} = 100 mD is equivalent to M_{ref,x} = 500 and K_{ref,x} = 10 mD that approximately corresponds to the lowpermeability layers whenever the permeability map is heterogeneous enough.

The considered permeability map is such that K_{z}/K_{x} = 0.1 or 1.

Permeability maps of typical layers thickness of 10 m or 5 m are considered.

The DykstraParsons coefficient that accounts for permeability heterogeneity mD is set to V_{DP} = 0.22, …, 0.75.
We now compare all the corresponding production indicators, NAPL recovery, chemicals and gas production, chemicals breakthrough time to identify trends. Specifically, special attention is paid to the thresholds in the DykstraParsons coefficient, such that when scaling M_{ref} significant differences are obtained due to a better vertical sweep, which is a consequence of flow diversion from most permeable layers by the play of the foam selective character, as previously reported in Sections 4 and 5.
Regarding NAPL recovery, we consider below the maximum variation in NAPL recovery during the transient regime under foam flooding while the cumulative mass of recovered chemicals and gas are taken at the end of Figure 1 injection sequences. All the reported production indicators variations when scaling M_{ref} are relative to the same injection type, as in the previous sections, all things being equal except a constant M_{ref}.
Figure 10a compares the maximum NAPL recovery relative variation when scaling M_{ref} relatively to a constant M_{ref}, as a function of the foam injection type (coinjection or alternate injection), V_{DP}, K_{ref,x,} K_{z}/K_{x} and L_{z}. Regardless of the configurations in the injection type, reference permeability, permeability anisotropy ratio, typical layers thickness, relative variation of NAPL recovery appears to approximately follow a power law of the DykstraParsons coefficient, above some V_{DP} threshold, as shown in Figure 10b in loglog scale.
This onset in the DykstraParsons coefficient for activating additional interlayer transfers is about V_{DP} ∼ 0.60.7 for all cases, roughly speaking. NAPL recovery improvement can be as high as about 35% for the most heterogeneous case with V_{DP} = 0.75. More precisely, this relative increase is about 20–35% in a coinjection mode for K_{ref,x} = 100 or 10 mD, and of about 25–40% in an alternate injection mode for K_{ref,x} = 100 or 10 mD, respectively.
The power law exponents reported in Figure 10b using a loglog scale approximately range from 8 to 13 and are close to 10 in terms of order of magnitude. The differential impact on NAPL recovery is very significant but less pronounced for thin layers with L_{z} = 5 m, as well as whenever permeability is isotropic with K_{z}/K_{x }= 1, which basically favors vertical flow.
The same comparison is reported in Figures 11a and 11b for the produced cumulative mass of chemicals and chemicals breakthrough time, respectively. No specific dependence on V_{DP} is observed as for NAPL recovery. However, a V_{DP} threshold (much less clearcut than for NAPL recovery) about V_{DP} ~ 0.6 emerges in all cases, such that for V_{DP} ≥ 0.6 there is a significant (i) reduction of produced chemicals for coinjection up to 40–60% for K_{ref,x} = 100 or 10 mD, and (ii) increase or little variation of produced chemicals for alternate injection, up to 5–12% for K_{ref,x} = 100 or 10 mD.
Chemicals breakthrough time, shown in Figure 11b exhibits very significant delays for coinjection only and for V_{DP} ≥ 0.5–0.6. It can be delayed up to 25–40% for K_{ref,x} = 100 or 10 mD. No delay is noticeable for the alternate injection. Again, the impact of scaling M_{ref} is still significant but less pronounced for thin layers with L_{z} = 5 m, as well as whenever permeability is isotropic with K_{z}/K_{x }= 1, which basically favors interlayer transfers.
The differential behaviour of gas production is reported in Figure 12. There are very significant differences compared to the constant M_{ref} cases as well, that mainly occur for V_{DP} ≥ 0.7 for coinjection and lead to a reduction of gas production or an increase of it or little variation otherwise. There are little differences for the alternate injection, and the impact is again less pronounced for thin layers with L_{z} = 5 m, as well as whenever permeability is isotropic. Generally speaking, however, the relative variations of gas production are much less marked than for the previous indicators.
Figures 13 and 14 provide an overview of these results by distinguishing coinjection from alternate injection for each DykstraParsons coefficient in order to quickly capture the key variations over the bundle of configurations of Table 3. A quick look will capture very significant variations for the most heterogeneous configurations for V_{DP} > 0.6 and mainly for NAPL recovery and surfactant amounts produced.
7 Discussion
This work, as a followup of previous work [15, 52], tested a new scaling law that drives the maximum mobility reduction yielded by foam of Darcyscale implicittexture foam models as a function of the permeability of the reservoir rock. This scaling law for the effect of permeability on foam properties was inferred from an analogy between foam flow in porous media and foam flow in capillary tubes and was found consistent with the modelling of available experimental data.
We have shown step by step, by analyzing the insitu saturation and composition profiles time evolution and the resulting production data, varying the degree of heterogeneity of a model stackoflayers porous medium, how such a result, newly implemented in IFP Energies nouvelles reservoir simulator Puma, impacts the hydrodynamic behaviour of foam and causes, in the most heterogeneous cases, a flow diversion from high to lowpermeability layers that increases the vertical sweep.
A sensitivity study was carried out to account for foam injection type (coinjection of surfactant solution and gas or surfactant solution alternating gas injection), layers thickness, permeability anisotropy and, of course, the degree of heterogeneity of the reservoir transcribed by the DykstraParsons coefficient.
The impact of this scaling law has been found moderate to very strong compared to simulations that do not account for it, according to the degree of reservoir heterogeneity. In a nutshell, the higher the heterogeneity is, the lower the sweep is but the higher the sweep improvement by foam is. Ignoring mobility reduction dependence on permeability in the foam process assessment of heterogeneous formations leads to an underestimation of mobility reduction benefits to improve flow conformance. As such, the question of cores selection is raised regarding foam evaluation in the laboratory and model calibration on selected rock types. Reservoir characterization remains more than ever a must.
NAPL relative incremental recovery, when scaling M_{ref}, in relation to a constant M_{ref}, has been found to follow a powerlaw behaviour as a function of the DykstraParsons coefficient. The V_{DP} thresholds above which these interlayer transfers are effectively activated have been identified. Chemicals production is very significantly impacted as well.
Even though the porous medium considered in this work remains a twodimensional toy model – it is however a deliberate choice in order to remain generic and moreover quite representative of the vertical sweep efficiency issue –, and given that the scaling law tested probably requires further verification and is not intended to apply to any rock type but rock type per rock type (there can very well be several regionalized K_{ref}, M_{ref} couples when considering multifacies geological formations), it is clear, given the variability observed as a function of the scaling performed, that no Darcean simulation of foam flow in heterogeneous porous media can be taken seriously without substantial management of heterogeneity in relation to foam. This could be a decisive contribution to the predictivity of Darcytype simulations for design purposes carried out prior to onsite tests, especially since investments are limited in the field of soil remediation. Good characterization of the considered porous medium and parameterization of the flow model, however imperfect and rustic it may be, are definitively unavoidable to achieve reliable predictions. In this respect, scaling laws, even approximate ones, should help guide smart and optimal sampling.
Surfactant adsorption and the impact of NAPL on foam, which were not taken into account in this work, should not significantly impact the observed trends. Indeed, the adsorption will have to be compensated by larger amounts of surfactant, and one can intuit that the antifoaming impact of NAPL, initially redhibitory when considering foam as a frontal displacement fluid, may be secondary due to the placement of the foam in most permeable areas already desaturated from NAPL. A major issue not addressed in this work remains injectivity. Indeed, there is no guarantee in practice that foams developing relatively high apparent viscosities can inject and propagate in porous formations due to pressure constraints at injection and in the reservoir [40].
Another concern is about possible crossflow between close layers when gas is blocked in a high permeability layer. In this case gas flows in the low permeability layer but will be in contact with the high permeability layer and may be diverted towards the high permeability zone. This happens when the layers are in capillary contact and may be enhanced by buoyancy effect. The considered model porous medium, with a relatively low permeability, is not typical of subsurface aquifers commonly encountered with NAPL problems, and is likely to minimize gravity effects compared to classical subsurface aquifers.
To finish with, it is worth mentioning that the foam flow transient regime studied in this work with a steadystate mobility reduction implicittexture model carries an ambiguity which remains to be clarified. The name, perhaps complicated on purpose, is indeed somewhat confusing, as this socalled steadystate model regarding foam actually accounts for the transient regime through the relative permeabilities interplay (the relative permeability to gas being modified in the presence of foam), like any standard Darcyscale multiphase flow model. It seems to us that the ambiguity, let us say rather lack of validation, can only be resolved by means of wellcontrolled transient experiments, by minimizing the model adjustable parameters of which in practice they are numerous. A flawless characterization of threephase flow relative permeabilities in the absence of foam is indeed required, possibly including hysteresis, as well as a thorough understanding and modelling of adsorption processes, devoid of any fudge factors. A certain research effort remains to be made in this area, apart from modelling.
Naturally fractured reservoirs with a poor recovery prognosis due to a preferential wettability to NAPL of the matrix rock are the most heterogeneous and present a similar problem to the one addressed in this work. Foam combines a reduction of fluid flow in the fracture network and the enhancement of fluids displacement in the matrix blocks thanks to the filtration of foam surfactant solution. Main expected effect is a forced twophase displacement of the matrix NAPL by foam surfactant filtrate. A possible (but not guaranteed) additional benefit is the reduction of oil trapping through a modification of interfacial rockfluids properties [20]. Needless to say, there is still room for improvement in dualmedium modelling that is necessary for performance purposes.
Notwithstanding, the authors hope that this work will contribute to a more reliable assessment of reservoirscale foam processes, whether for soil remediation, carbon capture, utilization, and storage, or enhanced oil recovery, and that it will be followed up in the community to advance the knowledge of foams in heterogeneous porous media. The scaling trends reported in Table 1 [15] certainly remain to be further investigated as well as the subsequent examination of their exploitation in practice. Considerable progress also remains to be made to accurately predict the distances over which foam holds during its propagation in porous media, in relation to the transient mechanisms of destruction and generation that depend on the porescale structure of the porous medium. In this respect, it is perhaps not Darcytype models that should be challenged first (populationbalance models were made precisely for this purpose) but rather the characterization of the porous medium at the pore scale and its transcription to the Darcy scale. This may be an impossible task in terms of predictivity at the geological scale of interest and considering reduced pilot scales of the order of a hundred or fifty meters seems appropriate. One may wonder how much faster this science would progress if more proprietary data of multiscale experiments were made public.
This layercake configuration is the most frequent and the simplest for sedimentary deposits [54]. In the context of this work, it is the minimum geological complexity needed in terms of permeability heterogeneity to highlight the selective behavior of foam in natural porous media.
This is known since the work of Leverett [55] for capillary pressure, being invariant for a given rocktype (θ is the contact angle). Regarding adsorption, a basic calculation indicates rock adsorption capacity scales as .
References
 Hirasaki G.J., Miller C.A., Szafranski R., Lawson J.B., Akiya N. (1997) Surfactant/foam process for aquifer remediation (Paper SPE 37 257), in: 1997 SPE International Symposium on Oilfield Chemistry, Houston, TX, USA, 18–21 February. [Google Scholar]
 Mamun C.K., Rong J.G., Kam S.I., Liljestrand H.M., Rossen W.R. (2002) Extending foam technology form improved oil recovery to environmantal remediation (Paper SPE 77557MS), in: SPE Annual Technical Conference and Exhibition, San Antonio, TX, USA, 29 September–2 October. [Google Scholar]
 Atteia O., Del Campo Estrada E., Bertin H. (2013) Soil flushing: a review of the origin of efficiency variability, Rev. Environ. Sci. Biotechnol. 12, 379–389. [CrossRef] [Google Scholar]
 Bertin H., Del Campo Estrada E., Atteia O. (2017) Foam placement for soil remediation, Environ. Chem. 14, 338–343. [CrossRef] [Google Scholar]
 Karthick A., Roy B., Chattopadhyay P. (2019) A review on the application of chemical surfacant and surfacant foam for remediation of petroleum oil contaminated soil, J. Envoron. Manage. 243, 187–205. [CrossRef] [Google Scholar]
 Flaifel H., Izadi M., Park S., Gupta I., Lee G., Kam S.I. (2020) Shallow subsurface environmental remediation by using tracersurfacantfoam processes: historymatching and performance prediction, Trans. Porous Media 134, 565–592. [CrossRef] [Google Scholar]
 Bouquet S., Douarche F., Roggero F., Leray S. (2020) Characterization of viscous fingering and channeling for the assessment of polymerbased heavy oil displacements, Trans. Porous Media 131, 873–906. [CrossRef] [Google Scholar]
 Puma is comarketed by BeicipFranlab and KAPPA Engineering, Available at https://www.beicip.com/reservoirsimulation and https://www.kappaeng.com/software/rubispuma/overview. [Google Scholar]
 Bernard G., Jacobs W.L. (1965) Effect of foam on trapped gas saturation and on permeability of porous media to water, Soc. Pet. Eng. J. 5, 4, 295–300. [CrossRef] [Google Scholar]
 Lawson J.B., Reisberg J. (1980) Alternate slugs of gas and dilute surfactant for mobility control during chemical flooding (Paper SPE8839MS), in: SPE/DOE Enhanced Oil Recovery Symposium, Tulsa, Oklahoma, April. [Google Scholar]
 Friedmann F., Chen W.H., Gauglitz P.A. (1991) Experimental and simulation study of hightemperature foam displacement in porous media, SPE Res. Eng. 6, 1, 37–45. [CrossRef] [Google Scholar]
 Peaceman D.W. (1977) Fundamentals of numerical reservoir simulation, volume 6 of Developments in Petroleum Science, Elsevier Science, Amsterdam. [Google Scholar]
 Trangenstein J.A., Bell J.B. (1989) Mathematical structure of the blackoil model for petroleum reservoir simulation, SIAM J. Appl. Math. 49, 3, 749–783. [Google Scholar]
 Marle C.M. (1981) Multiphase flow in porous media, 3rd edn., Gulf Publishing Company. [Google Scholar]
 Gassara O., Douarche F., Braconnier B., Bourbiaux B. (2020) Calibrating and scaling semiempirical foam flow models for the assessment of foambased EOR processes (in heterogeneous reservoirs), Trans. Porous Media 131, 1, 193–221. [CrossRef] [MathSciNet] [Google Scholar]
 Gassara O., Douarche F., Braconnier B., Bourbiaux B. (2017) Calibrating and interpreting implicittexture models of foam flow through porous media of different permeabilities, J. Pet. Sci. Eng. 159, 588–602. [Google Scholar]
 Gassara O., Douarche F., Braconnier B., Bourbiaux B. (2017) Equivalence between semiempirical and populationbalance foam models, Trans. Porous Media 120, 3, 473–493. [CrossRef] [MathSciNet] [Google Scholar]
 Weaire D., Hutzler S. (1999) The physics of foams, Oxford University Press. [Google Scholar]
 de Gennes P.G., BrochardWyart F., Quéré D. (2004) Capillarity and wetting phenomena: drops, bubbles, pearls, waves, SpringerVerlag. [CrossRef] [Google Scholar]
 Bouquet S., Douarche F., Roggero F., Bourbiaux B. (2020) Foam processes in naturally fractured carbonate oilwet reservoirs: technical and economic analysis and optimization, J. Pet. Sci. Eng 190, 107111. [CrossRef] [Google Scholar]
 Aziz K., Settari A. (1985) Petroleum reservoir simulation, Elsevier. [Google Scholar]
 Fayers F.J. (1989) Extension of Stone’s method 1 and conditions for real characteristics in threephase flow, SPE Res. Eng. 4, 4, 437–445. [CrossRef] [Google Scholar]
 Hirasaki G.J., Lawson J.B. (1985) Mechanisms of foam flow in porous media: apparent viscosity in smooth capillaries, SPE J. 25, 2, 176–190. [Google Scholar]
 Bretherton F.P. (1961) The motion of long bubbles in tubes, J. Fluid Mech. 10, 2, 166–188. [CrossRef] [MathSciNet] [Google Scholar]
 Ma K., Ren G., Mateen K., Morel D., Cordelier P. (2015) Modeling techniques for foam flow in porous media, SPE J. 20, 3, 453–470. [CrossRef] [Google Scholar]
 Ma K., Mateen K., Ren G., Luo H., Bourdarot G., Morel D. (2018) Mechanistic modeling of foam flow through porous media in the presence of oil: review of foamoil interactions and an improved bubble populationbalance model (Paper SPE 191 564), in: 2018 SPE Annual Technical Conference and Exhibition, Dallas, TX, USA, 24–26 September. [Google Scholar]
 Chen Q., Kovscek A.R., Gerritsen M. (2010) Modeling foam displacement with the localequilibrium approximation: theory and experimental verification, SPE J. 15, 1, 171–183. [CrossRef] [Google Scholar]
 Kovscek A.R., Radke C.J. (1994) Fundamentals of foam transport in porous media, in: M.J. Comstock, L.L. Schramm (eds), Foams: fundamentals and applications in the petroleum industry, American Chemical Society Advances in Chemistry. [Google Scholar]
 Patzek T.W. (1988) Description of foam flow in porous media by the population balance method, in: H.S. Duane (ed.), Surfactantbased mobility control, American Chemical Society Symposium Series, ACS, pp. 326–341, Available at https://pubs.acs.org/doi/pdf/10.1021/bk19880373.ch016. [Google Scholar]
 Falls A.H., Hirasaki G.J., Patzek T.W., Gauglitz D.A., Miller D.D., Ratulowski T. (1988) Development of a mechanistic foam simulator: the population balance and generation by snapoff, SPE Reserv. Eng. 3, 3, 884–892. [CrossRef] [Google Scholar]
 Kam S.I., Nguyen Q.P., Li Q., Rossen W.R. (2007) Dynamic simulations with an improved model for foam generation, SPE J. 12, 1, 35–48. [CrossRef] [Google Scholar]
 Kam S.I. (2008) Improved mechanistic foam simulation with foam catastrophe theory, Colloids Surf. A Physicochem. Eng. Asp. 318, 1–3, 62–77. [CrossRef] [Google Scholar]
 Kovscek A.R., Bertin H.J. (2003) Foam mobility in heterogeneous porous media, Trans. Porous Media 52, 1, 37–49. [CrossRef] [Google Scholar]
 Farajzadeh R., Eftekhari A. (2017) Effect of foam on liquid phase mobility in porous media, Sci. Rep. 43870, 3011–3018. [Google Scholar]
 Farajzadeh R., Lotfollahi M., Eftekhari A.A., Rossen W.R., Hirasaki G.J.H. (2015) Effect of permeability on implicittexture foam model parameters and the limiting capillary pressure, Energy Fuels 29, 5, 3011–3018. [Google Scholar]
 Khatib Z.I., Hirasaki G.J., Falls A.H. (1988) Effects of capillary pressure on coalescence and phase mobilities in foams flowing through porous media, SPE Res. Eng. 3, 3, 919–926. [CrossRef] [Google Scholar]
 Kapetas L., VincentBonnieu S., Farajzadeh R., Eftekhari A.A., MohdShafian S.R., Kamarul Bahrim R.Z., Rossen W.R. (2015) Effect of permeability on foammodel parameters – an integrated approach from coreflood experiments through to foam diversion calculations, in 18th European Improved Oil Recovery Symposium, Dresden, Germany, 14–16 April. [Google Scholar]
 Jian G., Zhang L., Da C., Puerto M., Johnston K.P., Biswal S.L., Hirasaki G.J. (2019) Evaluating the transport behavior of CO_{2} foam in the presence of crude oil under hightemperature and highsalinity conditions for carbonate reservoirs, Energy Fuels 33, 6038–6047. [CrossRef] [Google Scholar]
 Bergeron V., Radke C.J. (1992) Equilibrium measurements of oscillatory disjoining pressures in aqueous foam films, Langmuir 8, 12, 3020–3026. [CrossRef] [Google Scholar]
 Soulat A., Douarche F., Flauraud E. (2020) A modified well index to account for shearthinning behavior in foam EOR simulation, J. Pet. Sci. Eng. 191, 107146. [CrossRef] [Google Scholar]
 Peaceman D. (1978) Interpretation of wellblock pressures in numerical reservoir simulation, SPE J. 18, 3, 182–194. [Google Scholar]
 Batôt G., Delaplace P., Bourbiaux B., Pedroni L.G., Nabzar L., Douarche F., Chabert M. (2016) WAG management with foams: influence of injected gas properties and surfactant adsorption (Paper SPE 183 326), in: Abu Dhabi International Petroleum Exhibition & Conference, Abu Dhabi, UAE, 7–10 November. [Google Scholar]
 Brooks R.H., Corey A.T. (1966) Properties of porous media affecting fluid flow, J. Irrig. Drain. Div. 92, 2, 61–88. [CrossRef] [Google Scholar]
 Standing M.B. (1974) Notes on relative permeability relationships, Technical report, Division of Petroleum Engineering and Applied Geophysics, Norwegian Institute of Technology, University of Trondheim. Available at https://blasingame.engr.tamu.edu/z_zCourse_Archive/P620_18C/P620_zReference/PDF_Txt_Std_Kr_Notes_(1973).pdf. [Google Scholar]
 de Figueiredo L., Grana D., Le Ravalec M. (2020) Revisited formulation and applications of FFT moving average, Math. Geosci. 52, 6, 810–816. [Google Scholar]
 Le Ravalec M., Nœtinger B., Hu L.Y. (2000) The fft moving average (FFTMA) generator: an efficient numerical method for generating and conditioning Gaussian simulations, Mathe. Geol. 32, 6, 701–723. [CrossRef] [Google Scholar]
 Dykstra H., Parsons R.L. (1950) The prediction of oil recovery by waterflood, in: Secondary recovery of oil in the United States, American Petroleum Institute, Washington, DC, pp. 160–174. [Google Scholar]
 Craig F.F. (1993) The reservoir engineering aspects of waterflooding, volume 3 of SPE Monograph Series, Society of Petroleum Engineers. [Google Scholar]
 Willhite G.P. (1986) Waterflooding, volume 3 of SPE Textbook Series, Society of Petroleum Engineers. [Google Scholar]
 Lake L.W. (1989) Enhanced oil recovery, Prentice Hall. [Google Scholar]
 Green D.W., Willhite G.P. (1998) Enhanced oil recovery, volume 6 of SPE Textbook Series, Society of Petroleum Engineers. [Google Scholar]
 Douarche F., Braconnier B., Bourbiaux B. (2020) Scaling foam flow models in heterogeneous reservoirs for a better improvement of sweep efficiency (Paper ThB04), in: 17th European Conference of the Mathematics of Oil Recovery (ECMOR), Edinburgh, Scotland, 14–17 September. [Google Scholar]
 Pedroni L., Nabzar L. (2016) New insights on foam rheology in porous media, in: Rio Oil and Gas Expo and Conference Proceedings, Rio de Janeiro/Brazil, 24–27 October, (IBP. ISSN 2525 7560). [Google Scholar]
 Weber K.J., van Geuns L.C. (1990) Framework for constructing clastic reservoir simulation models, J. Pet. Technol. 42, 1248–1297. [CrossRef] [Google Scholar]
 Leverett M.C. (1941) Capillary behavior in porous solids (Paper SPE941152G), Trans. AIME 142, 152–169. [CrossRef] [Google Scholar]
Appendix
A.1 Estimation of M_{ref} from measurements
We have indicated in Section 2.1 that M_{ref} can be computed from the ratio of watergas coinjection steadystate measurements of differential pressure across coreplugs, Δp_{f}/Δp, with and without foam, as measured in the laboratory, given that Δp corresponds to a gas injection in the IT model according to (4). In fact it is not that simple, even assuming the reference injection leading to a stabilized steadystate Δp is a gas injection, in which case it would be tempting to consider that the ratio Δp_{f}/Δp corresponds to M_{ref}, where Δp_{f} denotes the steadystate pressure drop obtained under a coinjection of surfactant solution and gas. Here is why and how to estimate M_{ref} from such measurements. The point we wish to stress here is that the model M_{ref} is generally larger than the measured M_{ref}, which will be denoted M_{lab} hereafter. Moreover, in practice, it is common in the laboratory for reasons of convenience to compare the steadystate Δp_{f} obtained under foam injection (most often coinjection of surfactant solution and gas at a prescribed quality close to the optimal one), not to the Δp obtained under the equivalent coinjection without surfactant, but to a water or gas injection. We illustrate that with a few examples and orders of magnitude.
Let us calculate the maximum mobility reduction M_{ref} knowing the measured mobility reduction M_{lab} = Δp_{f}/Δp where Δp_{f} and Δp are respectively the steadystate pressure differences measured across coreplugs flooded by a surfactant solution/gas coinjection yielding foam and by the corresponding coinjection without surfactant, performed at the same quality f_{g} (i.e. with the same volumes of injected water and gas). Darcy’s law, for foam flow and for the corresponding coinjection without surfactant, at the same imposed total flow rate and quality, writes for the gas phase(A1)where U is the total velocity. Indeed, not only the relative permeability to gas is modified in the presence of foam (), but the saturation of the displacing gas front is also (). This is precisely what a BuckleyLeverett analytical calculation would predict [14]. From relations (A1) we get(A2)
Moreover, in the presence of foam, we have(A3)according to (4). Under optimal reference conditions, FM given by (5) writes, using (A3) (A4) i.e., using (A2) (A5)
The determination of M_{ref} is thus reduced to that of the saturations and , depending on S_{g}.
Determination of under foam injection. Darcy’s law, for foam flow and for the corresponding coinjection without surfactant, writes for the water phase(A6) i.e. by forming the ratio(A7)where denotes the reciprocal of k_{rw}.
Determination of S_{g} under coinjection. The gas saturation S_{g} verifies the relation(A8)
A similar calculation gives for a reference water injection(A9)and for a gas injection(A10)
Orders of magnitude. Powerlaw relative permeabilities (16) and (17) rewrite for a water and gas twophase flow in the absence of NAPL, setting S_{orw} = S_{org}(A11)with S_{G} = (S_{g} − S_{gr})/(1 − S_{wi} − S_{gr}). Having measured M_{lab}, to determine M_{ref} we compute S_{g} according to (A8), then following (A7) and finally M_{ref} from (A5). Some examples are given in Table A1. The M_{ref} values of 500 and 1581 assigned to K_{ref} = 100 mD in this work correspond to M_{lab} = 71 and 196 if the reference injection is a coinjection of water ans gas, and to M_{lab} = 272 and 746 if the reference injection is water flooding, respectively. These are quite reasonable values that are not abnormally large at all. They are even relatively low compared to what is generally measured [15, 35, 53].
Examples of M_{lab} to M_{ref} conversion for a reference injection which is eiher a coinjection of water and gas or a water injection (the gas considered is CO_{2} with T = 100 °C and p = 100 bar).
A.2 Comparison of insitu profiles with and without M_{ref} scaling
Fig. A1 Coinjection case after 1 PV of water injection prior to surfactant solution/gas coinjection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x} . Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 
Fig. A2 Coinjection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 
Fig. A3 Coinjection case final state, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 
Fig. A4 Coinjection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 1581, K_{ref,x} = 100 mD (or equivalently M_{ref,x} = 500, K_{ref,x} = 10 mD), V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 
Fig. A5 Coinjection case final state, with M_{ref,x} = 1581, K_{ref,x} = 100 mD (or equivalently M_{ref,x} = 500, K_{ref,x} = 10 mD), V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 
Fig. A6 Coinjection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 
Fig. A7 Coinjection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75 and L_{z} = 10 m. Left: scaled M_{ref,x} with K_{z}/K_{x} = 1. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x} with K_{z}/K_{x} = 0.1. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling (“ani” and “iso” stand for K_{z}/K_{x} = 0.1 and 1, respectively). 
Fig. A8 Coinjection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75 and K_{z}/K_{x} = 0.1. Left: scaled M_{ref,x} with L_{z} = 10 m. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x} with L_{z} = 5 m. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling (“L_{z}” and “L_{z}/2” stand for L_{z} = 10 and 5 m, respectively). 
Fig. A9 Alternate injection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, and K_{z}/K_{x} =0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 
Fig. A10 Alternate injection case post surfactant solution injection after approximately 0.9 PV of gas injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 
Fig. A11 Alternate injection case post surfactant solution injection after approximately 1.8 PV of gas injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 
All Tables
Examples of M_{lab} to M_{ref} conversion for a reference injection which is eiher a coinjection of water and gas or a water injection (the gas considered is CO_{2} with T = 100 °C and p = 100 bar).
All Figures
Fig. 1 Timeline of surfactant solution/gas coinjection (blue) and surfactant solution alternating gas injection (red) scenarios, where W, G and W* denote water, gas and surfactant solution, respectively. Timescale is indicated in injected pore volumes (PV, volume of fluid injected over the pore volume of the porous medium). Both cases have a prefoam water or gas injection history of 1 PV. In each case the volume of surfactant solution injected of 0.25 PV is identical: in the coinjection case, it is distributed over a total volume of gas and surfactant solution of 1.25 PV according to the imposed quality of f_{g} = 0.8 that corresponds to optimal foaming conditions. There is no such quality control for the surfactant solution alternating gas injection, where foam is formed by fingering of gas (injected after surfactant solution) through surfactant solution, or at least mixing of the two, due to the different mobility of the fluids. 

In the text 
Fig. 2 Water/NAPL and gas/NAPL twophase subsystems powerlaw relative permeabilities and capillary pressures (16) and (17) with parameters given in Table 2. 

In the text 
Fig. 3 Comparison of several twodimensional horizontal permeability field K_{x}(x, z) realizations (approximately lognormal distributed) with constant mean 〈K_{x}〉 = 100 mD and standard deviation , 50, 100, 150, 200, 250 mD, or equivalently, DykstraParsons coefficient V_{DP} = 0.22, 0.37, 0.57, 0.66, 0.72, 0.75. Typical layers thickness is driven by and close to the correlation length L_{z} = 10 m. 

In the text 
Fig. 4 Comparison of several twodimensional horizontal permeability field K_{x}(x, z) realizations (approximately lognormal distributed) with constant mean 〈K_{x}〉 = 100 mD and standard deviation , 50, 100, 150, 200, 250 mD, or equivalently, DykstraParsons coefficient V_{DP} =0.22, 0.37, 0.57, 0.66, 0.72, 0.75. Typical layers thickness is driven by and close to the correlation length L_{z} = 5 m. 

In the text 
Fig. 5 Coinjection case with M_{ref,x} = 500, K_{ref,x }= 100 mD, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Comparison of (a) NAPL recovery factor (RF, produced NAPL volume over initial NAPL volume in place, OOIP) and (b) produced surfactant concentration (C_{f}) and cumulative mass (M_{f}) time evolution as a function of V_{DP} for chemicalsfree injection (“nofoam”), foam with constant M_{ref,x} (“ref”, dashed lines) and scaled M_{ref,x} (“mod”, solid lines). Relative variations are either in relation to a chemicalsfree twin injection (“/ nofoam”) or a constant M_{ref,x} case (“/ foam”). 

In the text 
Fig. 6 Coinjection case with K_{ref,x} = 100 mD, (a) M_{ref,x} = 500, (b) M_{ref,x} = 1581 (or equivalently M_{ref,x} = 500, (a) K_{ref,x} = 100 mD, (b) K_{ref,x} = 10 mD), K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Comparison of produced cumulative gas volume (V_{g}) time evolution as a function of V_{DP} for chemicalsfree injection (“nofoam”), foam with constant M_{ref,x} (“ref”, dashed lines) and scaled M_{ref,x} (“mod”, solid lines). Relative variations are either in relation to a chemicalsfree twin injection (“/ nofoam”) or a constant M_{ref,x} case (“/ foam”). 

In the text 
Fig. 7 Coinjection case with M_{ref,x} = 1581, K_{ref} = 100 mD (or equivalently M_{ref,x} = 500, K_{ref,x} = 10 mD), K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Comparison of (a) NAPL recovery factor (RF, produced NAPL volume over initial NAPL volume in place, OOIP) and (b) produced surfactant concentration (C_{f}) and cumulative mass (M_{f}) time evolution as a function of V_{DP} for chemicalsfree injection (“nofoam”), foam with constant M_{ref,x} (“ref”, dashed lines) and scaled M_{ref,x} (“mod”, solid lines). Relative variations are either in relation to a chemicalsfree twin injection (“/ nofoam”) or a constant M_{ref,x} case (“/ foam”). 

In the text 
Fig. 8 Alternate injection case with M_{ref,x} = 500, K_{ref,x} = 100 mD, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Comparison of (a) NAPL recovery factor (RF, produced NAPL volume over initial NAPL volume in place, OOIP) and (b) produced surfactant concentration (C_{f}) and cumulative mass (M_{f}) time evolution as a function of V_{DP} for chemicalsfree injection (“nofoam”), foam with constant M_{ref,x} (“ref”, dashed lines) and scaled M_{ref,x} (“mod”, solid lines). Relative variations are either in relation to a chemicalsfree twin injection (“/ nofoam”) or a constant M_{ref,x} case (“/ foam”). 

In the text 
Fig. 9 Alternate injection case with M_{ref,x} = 1581, K_{ref,x}= 100 mD (or equivalently M_{ref,x} = 500, K_{ref,x} = 10 mD), K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Comparison of NAPL recovery factor (RF, produced NAPL volume over initial NAPL volume in place, OOIP) time evolution as a function of V_{DP} for chemicalsfree injection (“nofoam”), foam with constant M_{ref,x} (“ref”, dashed lines) and scaled M_{ref,x} (“mod”, solid lines). Relative variations are either in relation to a chemicalsfree twin injection (“/ nofoam”) or a constant M_{ref,x} case (“/ foam”). 

In the text 
Fig. 10 Comparison of the maximum NAPL recovery factor (RF, produced NAPL volume over initial NAPL volume in place, OOIP) relative variation under foam flooding with scaled M_{ref,x} relatively to constant M_{ref}, as a function of the injection sequence (coinjection or alternate injection), V_{DP}, K_{ref,x} = 100, 10 mD assigned to M_{ref,x} = 500 (or equivalently M_{ref,x} = 500, 1581 assigned to K_{ref,x} = 100 mD), K_{z}/K_{x} and L_{z}. 

In the text 
Fig. 11 Comparison of the (a) final produced surfactant cumulative mass (M_{f}) and (b) chemicals breakthrough time (BT) relative variation under foam flooding with scaled M_{ref,x} relatively to constant M_{ref,x}, as a function of the injection sequence (coinjection or alternate injection), V_{DP}, K_{ref,x} = 100, 10 mD assigned to M_{ref,x} = 500 (or equivalently M_{ref,x} = 500, 1581 assigned to K_{ref,x} = 100 mD), K_{z}/K_{x} and L_{z}. 

In the text 
Fig. 12 Comparison of the produced cumulative gas volume (V_{g}) relative variation under foam flooding with scaled M_{ref,x} relatively constant M_{ref,x}, as a function of the injection sequence (coinjection or alternate injection), V_{DP}, K_{ref,x} = 100, 10 mD assigned to M_{ref,x} = 500 (or equivalently M_{ref,x} = 500, 1581 assigned to K_{ref,x} = 100 mD), K_{z}/K_{x} and L_{z}. 

In the text 
Fig. 13 Comparison of the differential behaviour of production results under foam flooding with scaled M_{ref,x} relatively to constant M_{ref,x} as a function of the injection sequence (coinjection or alternate injection), V_{DP} = 0.22, 0.37, 0.57, K_{ref,x} = 100, 10 mD assigned to M_{ref,x} = 500 (or equivalently M_{ref,x} = 500, 1581 assigned to K_{ref,x} = 100 mD), K_{z}/K_{x} and L_{z} (subscripts “ani” and “iso” stand for K_{z}/K_{x} = 0.1 and 1 while “L_{z}” and “L_{z}/2” correspond to L_{z} = 10 and 5 m, respectively). 

In the text 
Fig. 14 Comparison of the differential behaviour of production results under foam flooding with scaled M_{ref,x} relatively to constant M_{ref,x} as a function of the injection sequence (coinjection or alternate injection), V_{DP} = 0.66, 0.72, 0.75, K_{ref,x} =100, 10 mD assigned to M_{ref,x} = 500 (or equivalently M_{ref,x} = 500, 1581 assigned to K_{ref,x} = 100 mD), K_{z}/K_{x} and L_{z} (subscripts “ani” and “iso” stand for K_{z}/K_{x} = 0.1 and 1 while “L_{z}” and “L_{z}/2” correspond to L_{z} = 10 and 5 m, respectively). 

In the text 
Fig. A1 Coinjection case after 1 PV of water injection prior to surfactant solution/gas coinjection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x} . Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 

In the text 
Fig. A2 Coinjection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 

In the text 
Fig. A3 Coinjection case final state, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 

In the text 
Fig. A4 Coinjection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 1581, K_{ref,x} = 100 mD (or equivalently M_{ref,x} = 500, K_{ref,x} = 10 mD), V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 

In the text 
Fig. A5 Coinjection case final state, with M_{ref,x} = 1581, K_{ref,x} = 100 mD (or equivalently M_{ref,x} = 500, K_{ref,x} = 10 mD), V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 

In the text 
Fig. A6 Coinjection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 

In the text 
Fig. A7 Coinjection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75 and L_{z} = 10 m. Left: scaled M_{ref,x} with K_{z}/K_{x} = 1. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x} with K_{z}/K_{x} = 0.1. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling (“ani” and “iso” stand for K_{z}/K_{x} = 0.1 and 1, respectively). 

In the text 
Fig. A8 Coinjection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75 and K_{z}/K_{x} = 0.1. Left: scaled M_{ref,x} with L_{z} = 10 m. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x} with L_{z} = 5 m. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling (“L_{z}” and “L_{z}/2” stand for L_{z} = 10 and 5 m, respectively). 

In the text 
Fig. A9 Alternate injection case after 0.25 PV of surfactant solution injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, and K_{z}/K_{x} =0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 

In the text 
Fig. A10 Alternate injection case post surfactant solution injection after approximately 0.9 PV of gas injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 

In the text 
Fig. A11 Alternate injection case post surfactant solution injection after approximately 1.8 PV of gas injection, with M_{ref,x} = 500, K_{ref,x} = 100 mD, V_{DP} = 0.75, K_{z}/K_{x} = 0.1 and L_{z} = 10 m. Left: constant M_{ref,x}. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled M_{ref,x}. Right: comparison of horizontally averaged profiles without (“ref”) and with (“mod”) M_{ref,x} scaling. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.