REGULAR ARTICLE

OPEN ACCESS

# Foam placement for soil remediation: scaling foam flow models in heterogeneous porous media for a better improvement of sweep efficiency

Frederic Douarche<fnref n="1" />,*, Benjamin Braconnier<fnref n="2" />, and Bernard Bourbiaux<fnref n="1" />

Received: 2 September 2022 / Accepted: 10 November 2023

**Abstract.** If injected with a large gas fraction, foam reduces mobility more in high-permeability layers and diverts flow to low-permeability 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 Darcy-scale numerical simulation. After briefly reviewing opportunities and challenges related to the use of foams in porous media and its Darcy-scale implicit-texture and population-balance 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 Darcy-type implicit-texture 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 population-balance 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 low-permeability layers. The threshold in permeability heterogeneity for which such a foam-driven 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 rock-typing strategy in relation to the porous medium characterization may guide a smart and optimal design of pre-feasibility laboratory campaign for foam evaluation, and the generalization of the findings to multi-facies geological formations are also discussed. As such, the use of physical foam mobility reduction scaling law is highly recommended for foam process evaluation.

**Keywords:** Multiphase flow, Porous media, Foam, Scaling laws, Heterogeneity, Models, Simulation, Soil remediation.

**Foreword:** This paper includes a large **Appendix B** which reports some selected in-situ 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.

<footnote n="1">¹IFP Energies nouvelles, Earth Sciences and Environmental Technologies Division, 1 et 4 avenue de Bois Préau, 92852 Rueil-Malmaison Cedex, France</footnote>

<footnote n="2">²IFP Energies nouvelles, Digital Sciences and Technologies Division, 1 et 4 avenue de Bois Préau, 92852 Rueil-Malmaison Cedex, France</footnote>

<footnote>*Corresponding author: frederic.douarche@ifpen.fr</footnote>
# 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 non-aqueous 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 anti-foaming 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 anti-foaming 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₂ storage in aquifer does not present the anti-foaming NAPL issue but is more questionable technically unless foam improves the effective storage capacity by smoothing out viscous fingering that may develop between CO₂ and brine and/or channeling due to heterogeneity through conformance control [7], provided foams are stable over the long-term.

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. Application-wise, 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 anti-foaming 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 large-scale 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 length-scales that can be sampled for subsurface characterization hence the set-up of the initial boundary value problem, and the resulting unavoidable uncertainties would remain. Even if we could integrate the pore-scale 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 rock-solid 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 trade-off is to build semi-empirical augmented generalized multiphase Darcy-type 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 case-by-case 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 near-subsurface 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 coarse-grain 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 permeability-dependent effective foam mobility, is quite directly and easily exploitable in the framework of augmented Darcy-type models addressing large-scale heterogeneous porous media, and indeed comply with the “if injected with a large gas fraction, foam reduces mobility more in high-permeability layers and diverts flow to low-permeability
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 hocus-pocus. 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 Darcy-type model (implicit-texture or population-balance). Section 3 is devoted to the set-up of a minimal model heterogeneous porous medium that is used in Sections 4 (detailed in-situ 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 B 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 Darcy-scale model types, driving equations and scaling relationships

## 2.1 Darcy-scale model types and driving equations

We start with general considerations about three-phase 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 non-aqueous phases, we consider a “black-NAPL” model (in reference to the “black-oil” model developed in hydrocarbon context [12, 13]) where the gas phase involves a modified velocity denoted $U_g^f$. The mass conservation equations read:

$$ 
\begin{aligned}
\partial_t(\Phi\rho_w S_w) + \nabla \cdot \mathbf{J}_w &= q_w, \\
\partial_t(\Phi\rho_w S_w C_w^s + (1-\Phi)\rho_r C_r^s) + \nabla \cdot \mathbf{J}_w^s &= q_w C_w^s, \\
\partial_t(\Phi C_h \rho_o S_o) + \nabla \cdot \mathbf{J}_o &= C_h q_o, \\
\partial_t(\Phi\rho_g S_g + \Phi C_v \rho_o S_o) + \nabla \cdot \mathbf{J}_g &= q_g + C_v q_o,
\end{aligned} 
\quad (1) $$

where $\Phi$ is the rock porosity and the fluxes $\mathbf{J}_w, \mathbf{J}_w^s, \mathbf{J}_o$ and $\mathbf{J}_g$ are defined as

$$ 
\begin{aligned}
\mathbf{J}_w &= \rho_w \mathbf{U}_w & \mathbf{J}_o &= \rho_o C_h \mathbf{U}_o \\
\mathbf{J}_w^s &= \rho_w \mathbf{U}_w C_w^s & \mathbf{J}_g &= \rho_o C_v \mathbf{U}_o + \rho_g \mathbf{U}_g^f.
\end{aligned} 
\quad (2) $$

For each phase denoted $\varphi = w, o, g, \mathbf{U}_\varphi$ is the Darcy velocity, $S_\varphi$ the saturation, $\rho_\varphi$ the density and $q_\varphi$ 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. $C_w^s$ stands for the flowing surfactant mass fraction in the water phase and $C_r^s$ for the adsorbed surfactant mass fraction on the rock with $\rho_r$ the rock mass density. Mobile and adsorbed surfactant mass fractions are related with an adsorption law such as the Langmuir isotherm.

Under creeping (low-velocity) flow conditions, the phase velocities in permeable porous media are governed by the generalized multiphase Darcy’s law [14]:

$$ \mathbf{U}_\varphi = -\lambda_\varphi \mathbf{K} \cdot (\nabla p_\varphi - \rho_\varphi \mathbf{g}) = \Phi S_\varphi \mathbf{V}_\varphi, \quad (3) $$

where $\mathbf{U}_\varphi$ is the Darcy velocity, $\mathbf{K}$ the porous medium permeability tensor, $\lambda_\varphi = k_{r\varphi}/\mu_\varphi$ the mobility, $\mu_\varphi$ the viscosity, $p_\varphi$ the pressure and $\mathbf{V}_\varphi = \mathbf{U}_\varphi/(\Phi S_\varphi)$ the interstitial velocity of phase $\varphi$ [14].

The flow of gas in the presence of foam is modelled differently whether a Population-Balance (PB) or a Semi-Empirical (SE)/Implicit-Texture (IT) modelling approach is used. These models are based on Darcy-type 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 $\lambda_g^f$ of IT model is scaled by a multi-parameter interpolation function FM assigned to the relative permeability to gas, whereas the gas viscosity is assumed unchanged whether foam is present or not:

$$ \mathbf{U}_g^f = -\lambda_g^f \mathbf{K} \cdot (\nabla p_g - \rho_g \mathbf{g}) \quad \lambda_g^f = \frac{k_{rg}^f}{\mu_g} = \frac{FM k_{rg}}{\mu_g}, \quad (4) $$

where $k_{rg}^f$ 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 multi-parameter interpolation functional form that includes the contributions of
physical parameters impacting the gas mobility reduction. FM is formulated as follows:

$$FM = \frac{1}{1 + (M_{\text{ref}} - 1) \prod_{i=1}^{4} F_i}, \qquad (5)$$

where $M_{\text{ref}}$ is the reference (maximum) gas mobility reduction under optimal conditions of the rock-fluid-additive 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

$$\begin{aligned}
F_1(C_w^s) &= \left( \frac{\min(C_w^s, C_{w,ref}^s)}{C_{w,ref}^s} \right)^{e_s}, \\
F_2(S_w) &= \frac{1}{2} + \frac{1}{\pi} \arctan [\Theta(S_w - S_w^*)], \\
F_3(S_o) &= \left( \frac{\max(0, S_o^* - S_o)}{S_o^*} \right)^{e_o}, \\
F_4(Ca) &= \left( \frac{Ca_{\text{ref}}}{\max(Ca, Ca_{\text{ref}})} \right)^{e_c},
\end{aligned} \qquad (6)$$

where $Ca = \mu_g U_g^f / \Phi \sigma_{wg}$ is the capillary number and $\sigma_{wg}$ the interfacial tension between water and gas. Positive numbers $e_s$, $\Theta$, $e_o$ and $e_c$ drive the shape of $F_1$, $F_2$, $F_3$ and $F_4$ functions, respectively, while $C_{w,ref}^s$, $S_w^*$, $S_o^*$ and $Ca_{\text{ref}}$ are thresholds. Specifically:

- $M_{\text{ref}}$ can be estimated from the ratio of steady-state measurements of differential pressure across core-plugs, $\Delta p_f / \Delta p$, with and without foam, as measured in the laboratory, given that $\Delta 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 $C_{w,ref}^s$ [18, 19].

- $F_2$ function governs the sharpness of the transition from the low- to the high-quality regime, through the dimensionless constant $\Theta$, as water saturation decreases in the vicinity of $S_w^*$ [15–17].

- $F_3$ function accounts for the impact of the anti-foaming character of NAPL on foam.

- $F_4$ function accounts for foam shear-thinning (or shear-thickening) behaviour through the capillary number Ca.

Examples of such functions can be found in [15–17, 20] and cited references therein. Estimation of $M_{\text{ref}}$ from measurements is detailed in Appendix A. Chemicals natural degradation with time can be accounted for through a half-life time associated to an exponential decay.

The three-phase flow relative permeability to NAPL is classically modelled as a combination of the underlying water/NAPL and gas/NAPL two-phase flow subsystems following Stone's models [21, 22]. Eventually, the system (1)–(6) with closure relationships $\sum S_\phi = 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_\phi(\mathbf{x}, t)$, $S_\phi(\mathbf{x}, t)$ and $C_w^s(\mathbf{x}, t)$.

Population-balance 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. Population-balance 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 population-balance 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 population-balance 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

$$\partial_t (\Phi S_g n_f) + \nabla \cdot (\mathbf{U}_g n_f) = \Phi S_g (r_g - r_c) + q_f, \quad (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 viscosity-type PB models extend the generalized Darcy's equation for the gas phase as follows:

$$\mathbf{U}_g^f = -\lambda_g^f \mathbf{K} \cdot (\nabla p_g - \rho_g \mathbf{g}) \qquad \lambda_g^f = \frac{k_{rg}}{\mu_g^f}, \qquad (8)$$

where $\mu_g^f$ 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 co-workers to a train of bubbles [23] thus leading to the following expression for the effective gas viscosity:

$$\mu_g^f = \mu_g + \frac{c_f n_f}{|\mathbf{V}_g|^{1/3}}, \qquad (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 $\mu_g^f$ expresses the shear-thinning effect of flow on foam bubbles. The exponent $\frac{1}{3}$ 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 shear-thinning 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_\phi$ and pressures $p_\phi$ 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 in-situ 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 steady-state 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 time-dependent generation and collapse, modelled by source terms $r_g$ and $r_c$, are intimately related to the pore-scale structure of the porous medium (involved in leave-behind, snap-off, 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 foam-based 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

$$ \mu_{\text{app}}^{f} = \frac{1}{\lambda_{w} + \lambda_{g}^{f}} = \frac{1}{\frac{k_{rw}}{\mu_{w}} + \frac{FM k_{rg}}{\mu_{g}}}, \quad (10) $$

and depends on foam quality $f_g = |U_g|/|U|$, where the total velocity $U = \sum U_\phi$ is the sum of the phase velocities. When plotted as a function of $f_g$, $\mu_{\text{app}}^f$ has the shape of an asymmetrical bell curve that is maximum for the optimal quality $f_g^*$ in the neighborhood of which the sharpness of the transition from the low- to the high-quality 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 two-layer 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 implicit-texture 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_{\text{ref}}$ represents the ratio between the continuous gas mobility $k_{rg}/\mu_g$ and the foamed-gas mobility that is defined as $k_{rg}^f/\mu_g$ in IT models, or as $k_{rg}/\mu_g^f$ in PB models, at optimal foaming conditions. Scaling $M_{\text{ref}}$ for different porous media characterized by the same $k_{rg}$ function can then be reduced to the scaling of $\mu_g^f$. 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 foamed-gas.

Kovscek and Bertin [33] derived an expression of $\mu_g^f$ 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 pre-determined volume flowing in a capillary tube, that is [23]:

$$ \mu_g^f = 0.85 \mu_w \frac{n_f}{R_c} (R^2 + R_c^2) \left( \frac{\sigma_{wg}}{3\mu_w V_g} \right)^{\frac{1}{3}}, \quad (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 flow-length 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.
<table><caption>**Table 1.** Scaling trends summary for homothetic porous media [15].</caption><thead><tr><th>Theoretical model for films stability</th><th>Static picture</th><th>Dynamic picture</th><th>Observation</th></tr></thead><tbody><tr><td>Expected evolution of p<sup>*</sup><sub>c</sub> with K</td><td>Constant p<sup>*</sup><sub>c</sub></td><td>p<sup>*</sup><sub>c</sub> ~ 1/√K</td><td>p<sup>*</sup><sub>c</sub> ~ 1/K<sup>α</sup> with α ≤ 1/2</td></tr><tr><td>Expected evolution of S<sup>*</sup><sub>w</sub> with K</td><td>S<sup>*</sup><sub>w</sub> increase with K decrease</td><td>Constant S<sup>*</sup><sub>w</sub></td><td>Moderate S<sup>*</sup><sub>w</sub> increase with K decreaseas S<sup>*</sup><sub>w</sub> ~ K<sup>λ</sup>(x-1/2)</td></tr><tr><td>Expected evolution of Θ with K</td><td>p<sub>c</sub>-dependent</td><td>Constant Θ</td><td>Quasi-constant</td></tr><tr><td>Expected evolution of M<sub>ref</sub> with K</td><td></td><td>M<sub>ref</sub> ~ √K</td><td></td></tr></tbody></table>

For a given foam flowing in a porous medium of permeability K and porosity Φ, capillary tube radius R is equivalent to the quantity $\sqrt{8K/\Phi}$ 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 very-low-permeability media). Therefore, $\mu_g^f$ scales as $n_f R^2$ under given velocity conditions. For homothetic porous media of different permeabilities with similar porosities, R scales as $\sqrt{K}$ and the lineic foam texture $n_f$ that is proportional to $1/R$ [15], scales as $1/\sqrt{K}$, hence $\mu_g^f$ scales as $\sqrt{K}$. The above analysis indicates that $M_{ref}$ is expected to increase as the square root of permeability for homothetic porous media of different permeabilities:

$$M_{ref} \sim \sqrt{K}. \qquad (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 rock-types in reservoir models.

# 3 Model case study set-up

## 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<fnref n="1" />, high-permeability layers may lead to a by-pass of low-permeability ones when a fluid (say water or gas) is injected to flush NAPL. Thus, remaining NAPL in low-permeability 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 low-permeability layers. It is precisely the diverted fluid (water and/or gas), and not necessarily foam itself, that will displace the trapped NAPL in low-permeability zones. It is thus more a problem of conformance control than of displacement. Before moving on to the set-up 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, rock-type per rock-type.

For the sake of demonstration, let us assume the layer-cake reservoir overall average horizontal permeability is $\langle K_x \rangle = \sum h_i K_x^i / H = 100$ mD where $h_i$ is the ith layer height, $K_x^i$ is the ith layer horizontal permeability and $H = \sum h_i$ is the height of the reservoir interval. Let us further assume the interpretation of laboratory foam displacements on core-plugs yields, among other things, a reference mobility reduction of $M_{ref} = 500$ for that permeability $\langle K_x \rangle = 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

<footnote n="1">¹ This layer-cake 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.</footnote>
much on the rock (this is not the question asked here but these are strong assumptions).

In the minimal layer-cake 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<fnref n="2" /> scale as 1/√K.

Let us then consider a layer-cake reservoir of fixed average horizontal permeability, say $\langle K_x \rangle = 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 rock-types 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:

$$ \mathcal{M} = M_{ref} \sqrt{K/K_{ref}}. \quad (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 layer-cake 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:

$$ \mathcal{M}_x = M_{ref,x} \sqrt{K_x/K_{ref,x}} \quad \alpha = x, y, z. \quad (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:

$$ FM_x = \frac{1}{1 + (\mathcal{M}_x - 1) \prod_{i=1}^{4} F_i} \quad \alpha = x, y, z. \quad (15) $$

As the vertical permeability $K_z$ usually differs from the horizontal one $K_x$ (assuming $K_x = K_y$ for the considered simple layer-cake 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 oil-wet 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 build-up. This is in practice one of the key points to be addressed for any viscous fluid injection in low-permeability 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 shear-thinning (or shear-thickening) fluids whose viscosity varies a lot for large velocity gradients in the near-wellbore, as it can be the case of pre-formed 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 shear-thinning 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 $1/10$, respectively, where $K_z(x, z)$ is deduced from $K_x(x, z)$ which is geostatistically generated by imposing a constant mean $\langle K_x \rangle = 100$ mD, as explained in Section 3.3. Two types of injection are considered: either a co-injection of surfactant solution and gas (with a quality of $f_g = 0.8$) which corresponds to a pre-formed foam injection, and surfactant solution alternating gas injection for which foam is formed in-situ, as reported in Figure 1, with an imposed constant total injection flow rate of 54.8 m³/day for which one pore volume (PV) of the reservoir is injected within one year. Indeed, since foams are not a proven large-scale recovery process, the mode of fluids injection is also being investigated, with the literature reporting roughly as many co-injection tests as alternate injection. However, pre-formed foam injection may be preferred for fractured reservoirs with poor in-situ 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 above-mentioned foam features, including the $M_{ref}$ scaling law relationships (14) and (15) (see [15] and cited references therein).

## 3.2 Rock-fluid properties

The reservoir under consideration has an average temperature of $T_{res} = 100$ °C and an average pressure of

<footnote n="2">² This is known since the work of Leverett [55] for capillary pressure, $J = \sqrt{K/\Phi p_c/\sigma \cos \theta}$ being invariant for a given rock-type ($\theta$ is the contact angle). Regarding adsorption, a basic calculation indicates rock adsorption capacity scales as $\Phi^{3/2}/\sqrt{K}$.</footnote>
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Timeline of surfactant solution/gas co-injection (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 pre-foam 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 co-injection 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.</figcaption></figure>

$p_{res}$ = 100 bar. Rock-fluid properties are representative of Berea sandstone with a foam made up of a surfactant aqueous solution and a methane/CO₂ mixture (62% CO₂, 38% CH₄) with a quality $f_g = 0.8$. The reader is referred to the previous work [42] for further details. These conditions are far from near-subsurface 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.

Three-phase flow is modelled as a combination of the underlying water/NAPL and gas/NAPL two-phase subsystems following Stone I model [21, 22]. Water is considered as the wetting phase while gas is the non-wetting phase, NAPL being the common and intermediate phase to both subsystems. Each two-phase subsystem is set up symmetrically with respect to its wetting phase following Brooks-Corey relationships [43, 44], assuming power-law 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:

$$ \left\{ \begin{array}{l} k_{rw}(S_W) = \kappa_w (S_W)^{n_w}, \\ k_{row}(S_W) = \kappa_{ow} (1 - S_W)^{n_{ow}}, \\ p_{cow}(S_W) = \frac{\pi_{ow}}{(S_W)^{\frac{1}{k_{ow}}}}, \end{array} \right. \qquad (16) $$

where $S_{wi}$ is the irreducible water saturation, $S_{orw}$ is the residual NAPL saturation with respect to water, $\kappa_w$ and $\kappa_{ow}$ are the maximum relative permeabilities to water and NAPL, and $\pi_{ow}$ is the entry capillary pressure. The gas/NAPL subsystem, for which NAPL is the wetting

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Water/NAPL and gas/NAPL two-phase subsystems power-law relative permeabilities and capillary pressures (16) and (17) with parameters given in Table 2.</figcaption></figure>

phase, is set up symmetrically as a function of the normalized mobile gas saturation $S_G = (S_g - S_{gr})/(1 - S_{orw} - S_{wi} - S_{gr})$:

$$ \left\{ \begin{array}{l} k_{rg}(S_G) = \kappa_g (S_G)^{n_g}, \\ k_{rog}(S_G) = \kappa_{og} (1 - S_G)^{n_{og}}, \\ p_{cog}(S_G) = \frac{\pi_{og}}{(1 - S_G)^{\frac{1}{k_{og}}}}, \end{array} \right. \qquad (17) $$

with $p_{cog} = p_g - p_o$ and where $S_{gr}$ is the irreducible gas saturation and $S_{orw}$ is the residual NAPL saturation with respect to gas. Figure 2 shows flow curves (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 propery transport and adsorption, that are out of the scope of the present work.

### 3.3 Model layered porous medium set-up

The reservoir under consideration is a two-dimensional cross-sectional x-z layered porous medium of 200 m length, 50 m height and 10 m width<fnref n="3" />. Porosity is set constant to $\Phi = 0.2$. The average horizontal permeability in the x-coordinate is set constant to $\langle K_x \rangle = 100$ mD while the permeability anisotropy ratio $K_z/K_x$, from which the vertical permeability $K_z$ is deduced, is set to either $\frac{1}{10}$ (strongly anisotropic) or 1 (isotropic). Several permeability map

<footnote n="3">³ Since we focus on vertical sweep efficiency, width is quite arbitrary here and could be of 100 m as well, subject to a recalculation of the injection flow rate in proportion to the porous medium volume.</footnote>
<table><caption>**Table 2.** Parameters of water/NAPL and gas/NAPL two-phase subsystems power-law relative permeabilities and capillary pressures (16) and (17).</caption><thead><tr><th colspan="2">Water-NAPL</th><th colspan="2">Gas-NAPL</th></tr></thead><tbody><tr><td>S<sub>orw</sub></td><td>0.3</td><td>S<sub>org</sub></td><td>0.1</td></tr><tr><td>S<sub>wi</sub></td><td>0.15</td><td>S<sub>gr</sub></td><td>0</td></tr><tr><td>κ<sub>ow</sub></td><td>0.9</td><td>κ<sub>og</sub></td><td>0.9</td></tr><tr><td>κ<sub>w</sub></td><td>0.2</td><td>κ<sub>g</sub></td><td>0.2</td></tr><tr><td>n<sub>ow</sub></td><td>2</td><td>n<sub>og</sub></td><td>4</td></tr><tr><td>n<sub>w</sub></td><td>4</td><td>n<sub>g</sub></td><td>2</td></tr><tr><td>π<sub>ow</sub> (bar)</td><td>0.01</td><td>π<sub>og</sub> (bar)</td><td>0.04</td></tr><tr><td>λ<sub>ow</sub></td><td>2</td><td>λ<sub>og</sub></td><td>2</td></tr></tbody></table>

realizations $K_x(x, z)$ of increasing degree of heterogeneity are considered. They range from quasi-homogeneous to highly heterogeneous by varying the horizontal permeability standard deviation $\sigma_{K_x} = 25, 50, 100, 150, 200, 250$ mD, as reported in Figures 3 and 4, using a fast Fourier transform moving average algorithm [45, 46].

In each case, $K_x(x, z)$ is approximately log-normal distributed with a constant mean $\langle K_x \rangle = 100$ mD. Layers typical thickness over which $K_x(x, z)$ is approximately constant is controlled by a Gaussian auto-correlation function whose correlation length $L_z$ in the z-coordinate is set to 10 or 5 m, while the correlation length in the x-coordinate is arbitrarily large ($L_x = 1000$ m) in order to render a layer-cake 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 Dykstra-Parsons coefficient [47] defined as $V_{DP} = (e^{\mu} - e^{\mu - \sigma}) / e^{\mu} = 1 - e^{-\sigma})$ for a log-normal random variable $K$ of parameters $\mu = \langle \ln K \rangle$ and $\sigma^2 = \langle (\ln K - \mu)^2 \rangle$. $V_{DP}$ writes in terms of the permeability mean $\langle K \rangle$ and variance $\sigma_K^2 = \langle (K - \langle K \rangle)^2 \rangle$ as<fnref n="4" />:

$$V_{DP} = 1 - \exp\left(-\sqrt{\ln\left[1 + \left(\frac{\sigma_K}{\langle K \rangle}\right)^2\right]}\right), \quad (18)$$

where the subscript $x$ has been dropped for the sake of readability. Thus, $V_{DP} \to 0$ for a homogeneous permeability distribution such that $\sigma_K/\langle K \rangle \to 0$ whereas $V_{DP} \to 1$ for an extremely heterogeneous distribution with $\sigma_K/\langle K \rangle \to \infty$. Above-mentioned permeability standard deviations of $\sigma_{K_x} = 25, 50, 100, 150, 200, 250$ mD are equivalent to the Dykstra-Parsons 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

z-coordinate with 100 and 50 cells, respectively; cells are therefore of constant size of $2 \times 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 \in [0, H]$. No-flow conditions are assigned to the top and the bottom of the layer-cake $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 aquifer-type reservoir interval to heterogeneous puzzle-type 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 pre-formed under co-injection of surfactant solution and gas, or forms in-situ 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 (co-injection or alternate injection) have a pre-foam 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 co-injection 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 two-dimensional $x$-$z$ layer-cake porous medium under consideration is initially saturated with NAPL up to the irreducible water saturation<fnref n="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 $\langle K_x \rangle$. 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.

<footnote n="4">⁴ For a log-normal random variable $X$ with probability density function $f_X(x) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-(\ln x - \mu)^2 / 2\sigma^2}$ of parameters $\mu = \langle \ln X \rangle$ and $\sigma^2 = \langle (\ln X - \mu)^2 \rangle$, where $\langle \cdot \rangle$ is the mathematical expectation, one has $\langle X^n \rangle = e^{n\mu + \frac{1}{2}n^2\sigma^2}$ and $\sigma_X^2 = \langle (X - \langle X \rangle)^2 \rangle = e^{2\mu + \sigma^2}(e^{\sigma^2} - 1) = \langle X \rangle^2(e^{\sigma^2} - 1)$, hence (18).</footnote>

<footnote n="5">⁵ In practice, aquifers are not fully saturated with NAPL but at a much lower saturation. However, this simplified initial condition has been considered to better highlight channeling issues and sweep efficiency.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Comparison of several two-dimensional horizontal permeability field $K_x(x, z)$ realizations (approximately log-normal distributed) with constant mean $\langle K_x \rangle = 100$ mD and standard deviation $\sigma_{K_x} = 25, 50, 100, 150, 200, 250$ mD, or equivalently, Dykstra-Parsons 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.</figcaption></figure>

- 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 $1/10$, where $K_z(x, z)$ is deduced from $K_x(x, z)$ which is geostatistically generated by imposing a constant mean $\langle K_x \rangle = 100$ mD, as explained in Section 3.3.

- Several two-dimensional horizontal permeability field $K_x(x, z)$ realizations (approximately log-normal distributed) are considered by imposing the mean $\langle K_x \rangle = 100$ mD and varying the standard deviation $\sigma_{K_x} = 25, 50, 100, 150, 200, 250$ mD, or equivalently, the Dykstra-Parsons 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, *in-situ* foam flows without and with $M_{\text{ref}}$ scaling are compared in the two following Sections 4.2 and 4.3 for the largest Dykstra-Parsons coefficient $V_{DP} = 0.75$ for foam flooding obtained under co-injection of surfactant solution and gas or surfactant solution alternating gas

injection. The initial *in-situ* behaviour under water or gas flooding is first analyzed in Section 4.1. Figures that compare the *in-situ* dynamic behaviour of saturations $S_w, S_g, S_o$, foaming agent concentration $C_f$, mobility reduction function in the x-coordinate $(FM_x)^{-1}$, as well as the horizontally and vertically averaged corresponding profiles, have been moved in Appendix B 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_{\text{ref}}$ scaling are compared in Section 6 and characterized by sweep trends in relation with the Dykstra-Parsons 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, co-injection of water and gas or water alternating gas injection, without surfactant in the water phase. This baseline will serve as a comparison
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Comparison of several two-dimensional horizontal permeability field $K_x(x, z)$ realizations (approximately log-normal distributed) with constant mean $\langle K_x \rangle = 100$ mD and standard deviation $\sigma_{K_x} = 25, 50, 100, 150, 200, 250$ mD, or equivalently, Dykstra-Parsons 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.</figcaption></figure>

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 well-known 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 late-time evolution reported in Figures 5a and 8a which reflects the late-time 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.

We will not go into further detail on this well-known 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 rock-fluids 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 early-time NAPL recovery: water indeed presents a much larger displacement saturation front than gas
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Co-injection case with $M_{\text{ref},x} = 500$, $K_{\text{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_{\text{DP}}$ for chemicals-free injection ("nofoam"), foam with constant $M_{\text{ref},x}$ ("ref", dashed lines) and scaled $M_{\text{ref},x}$ ("mod", solid lines). Relative variations are either in relation to a chemicals-free twin injection ("/ nofoam") or a constant $M_{\text{ref},x}$ case ("/ foam").</figcaption></figure>
(assumed here to be immiscible with NAPL), mostly due to the difference in viscosity, thus the transient early-time 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. Co-injection 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 co-injecting 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 inter-well 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 *in-situ* foam flow with and without $M_{ref}$ scaling: co-injection

We now focus on foam injection in the form of a 1.25 PV co-injection 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 co-injection of surfactant solution and gas consists of a 0.25 PV surfactant solution slug. The behaviour of the pre-foam water injection has been commented in the previous section. Let us see now what the foam brings.

Figures A1, A2 and A3 compare the *in-situ* dynamic behaviour 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 high-permeability 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 Buckley-Leverett calculation [14, 49–51].

Alternatively, it can be noted qualitatively from the saturation profiles that the *in-situ* 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 core-scale 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 $f_g^*$. In any case this is monitorable by simulation provided that a reliable and accurate characterization of the rock-fluids system is available, and that on-site 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 low-permeable 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 inter-layer 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 inter-layer vertical transfers.

Although we do not report *in-situ* 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 *in-situ* foam flow with and without $M_{ref}$ scaling: alternate injection

Let us now move on to foam flow where foam is formed *in-situ* by alternating surfactant solution and gas slugs, as indicated in Figure 1. Unlike the co-injection 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 co-injection case. The behaviour of the pre-foam gas injection has been already commented in Section 4.1.

Figures A9, A10 and A11 compare the *in-situ* 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 co-injection 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 co-injection, 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 co-injection is the less surfactant trapping in low-permeability 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 co-injection, 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 Dykstra-Parsons coefficient over the range $V_{DP} = 0.22, \dots, 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 chemicals-free injection type when scaling $M_{ref}$ or not (co-injection with chemicals is compared to co-injection 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}$ (co-injection with $M_{ref}$ scaling is compared to co-injection 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 co-injection (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: co-injection

To begin with, we note the extreme sensitivity of the NAPL recovery of a chemicals-free injection (in absence of foam) to the Dykstra-Parsons 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_{\text{ref}}$ is scaled or fixed are all the higher as the degree of heterogeneity, quantified by $V_{\text{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_{\text{ref}}$ with $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 100$ mD, for $V_{\text{DP}} = 0.72$ and 0.75, respectively and relatively to a foam injection with a constant $M_{\text{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_{\text{DP}} \ge 0.57$ when $M_{\text{ref}}$ is scaled with $M_{\text{ref},x} = 500$ and $K_{\text{ref},x} = 100$ mD, as reported in Figure 5b. The higher $V_{\text{DP}}$ is, the earlier the chemicals breakthrough and the larger the amount of recovered chemicals are for a constant $M_{\text{ref}}$. This trend remains when scaling $M_{\text{ref}}$ but there is a significant delay in chemicals breakthrough and a reduction in recovered chemicals of about 20–40% for large $V_{\text{DP}}$ values ($V_{\text{DP}} \ge 0.57$) over constant $M_{\text{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 low-permeability layers, as explained in Section 4.2.

To finish with, the produced volume of gas is sensitive to $V_{\text{DP}}$ but much less than NAPL and surfactant recovery, as reported in Figure 6a. It is lowered by 5% for $V_{\text{DP}} = 0.75$ when scaling $M_{\text{ref}}$ compared to a constant $M_{\text{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_{\text{ref},x} = 1581$ and $K_{\text{ref},x} = 100$ mD (or equivalently $M_{\text{ref},x} = 500$ and $K_{\text{ref},x} = 10$ mD) are similar and more pronounced when $M_{\text{ref}}$ is scaled compared to a constant $M_{\text{ref}}$, as reported in Figures 6b, 7a and 7b, and as already observed from the *in-situ* profiles in Section 4.2. Let us give some orders of magnitude: when scaling $M_{\text{ref}}$, transient NAPL recovery is improved by about 35% and 20% for $V_{\text{DP}} = 0.75$ and 0.72, respectively, while chemicals production is lowered by about 50–60% for $V_{\text{DP}} \ge 0.57$. Gas production does not exhibit significant differences compared to the $M_{\text{ref},x} = 500$ mD and $K_{\text{ref},x} = 100$ mD case.

## 5.2 Comparison of foam flow production results with and without $M_{\text{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_{\text{ref},x} = 500$ and $K_{\text{ref},x} = 100$ mD, and in Figure 9 for $M_{\text{ref},x} = 1581$ and $K_{\text{ref},x} = 100$ mD (or equivalently $M_{\text{ref},x} = 500$ and $K_{\text{ref},x} = 10$ mD).

A critical dependence of NAPL recovery on $V_{\text{DP}}$ for a chemical-free injection is again noted: as $V_{\text{DP}}$ increases from 0.22 to 0.75, NAPL recovery drops from 70% to 20%, approximately. By scaling $M_{\text{ref}}$ it is improved by about 25%, 10% and 7% relatively to a constant $M_{\text{ref}}$ for $V_{\text{DP}} = 0.75$, 0.72 and 0.66 when $M_{\text{ref},x} = 500$ and $K_{\text{ref},x} = 100$ mD,

and by 40%, 23% and 10% relatively to a constant $M_{\text{ref}}$ for the same $V_{\text{DP}}$ values when $M_{\text{ref},x} = 1581$ and $K_{\text{ref}} = 100$ mD. Again, the more heterogeneous the reservoir is (the higher $V_{\text{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_{\text{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 co-injection 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_{\text{ref},x} = 500$ or 1581 is assigned to the reference horizontal permeability $K_{\text{ref},x} = 100$ mD that fits the average horizontal permeability $\langle K_x \rangle$. The case $M_{\text{ref},x} = 1581$ and $K_{\text{ref},x} = 100$ mD is equivalent to $M_{\text{ref},x} = 500$ and $K_{\text{ref},x} = 10$ mD that approximately corresponds to the low-permeability 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 Dykstra-Parsons coefficient that accounts for permeability heterogeneity $\sigma_{K_x} = 25, \dots, 250$ mD is set to $V_{\text{DP}} = 0.22, \dots, 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 Dykstra-Parsons coefficient, such that when scaling $M_{\text{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_{\text{ref}}$ are relative to the same injection type, as in the previous sections, all things being equal except a constant $M_{\text{ref}}$.
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Co-injection case with $K_{\text{ref},x} = 100$ mD, (a) $M_{\text{ref},x} = 500$, (b) $M_{\text{ref},x} = 1581$ (or equivalently $M_{\text{ref},x} = 500$, (a) $K_{\text{ref},x} = 100$ mD, (b) $K_{\text{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_{\text{DP}}$ for chemicals-free injection ("nofoam"), foam with constant $M_{\text{ref},x}$ ("ref", dashed lines) and scaled $M_{\text{ref},x}$ ("mod", solid lines). Relative variations are either in relation to a chemicals-free twin injection ("/ nofoam") or a constant $M_{\text{ref},x}$ case ("/ foam").</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Co-injection case with $M_{\text{ref},x} = 1581$, $K_{\text{ref}} = 100$ mD (or equivalently $M_{\text{ref},x} = 500$, $K_{\text{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_{\text{DP}}$ for chemicals-free injection ("nofoam"), foam with constant $M_{\text{ref},x}$ ("ref", dashed lines) and scaled $M_{\text{ref},x}$ ("mod", solid lines). Relative variations are either in relation to a chemicals-free twin injection ("/ nofoam") or a constant $M_{\text{ref},x}$ case ("/ foam").</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Alternate injection case with $M_{\text{ref},x} = 500$, $K_{\text{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_{\text{DP}}$ for chemicals-free injection ("nofoam"), foam with constant $M_{\text{ref},x}$ ("ref", dashed lines) and scaled $M_{\text{ref},x}$ ("mod", solid lines). Relative variations are either in relation to a chemicals-free twin injection ("/ nofoam") or a constant $M_{\text{ref},x}$ case ("/ foam").</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Alternate injection case with $M_{\text{ref},x} = 1581$, $K_{\text{ref},x} = 100$ mD (or equivalently $M_{\text{ref},x} = 500$, $K_{\text{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_{\text{DP}}$ for chemicals-free injection ("nofoam"), foam with constant $M_{\text{ref},x}$ ("ref", dashed lines) and scaled $M_{\text{ref},x}$ ("mod", solid lines). Relative variations are either in relation to a chemicals-free twin injection ("/ nofoam") or a constant $M_{\text{ref},x}$ case ("/ foam").</figcaption></figure>

<table><caption>**Table 3.** Summary of the configurations studied. Injection type refers to the way foam is formed, either *ex-situ* by co-injection of surfactant solution and gas or *in-situ* by surfactant solution alternating gas injection, as shown in Figure 1. $M_{\text{ref},x}$ is assigned to $K_{\text{ref},x} = \langle K_x \rangle = 100$ mD.</caption><thead><tr><th>Variable</th><th>NB of cases</th><th>Cases</th></tr></thead><tbody><tr><td>Injection type</td><td>2</td><td>Co-injection/alternate</td></tr><tr><td>L<sub>z</sub> (m)</td><td>2</td><td>10/5</td></tr><tr><td>K<sub>z</sub>/K<sub>x</sub></td><td>2</td><td>1/10</td></tr><tr><td>V<sub>DP</sub></td><td>6</td><td>0.22, ..., 0.75</td></tr><tr><td>M<sub>ref,x</sub></td><td>2</td><td>500/1581</td></tr></tbody></table>
Figure 10a compares the maximum NAPL recovery relative variation when scaling $M_{\text{ref}}$ relatively to a constant $M_{\text{ref}}$, as a function of the foam injection type (co-injection or alternate injection), $V_{\text{DP}}$, $K_{\text{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 Dykstra-Parsons coefficient, above some $V_{\text{DP}}$ threshold, as shown in Figure 10b in log-log scale.

This onset in the Dykstra-Parsons coefficient for activating additional inter-layer transfers is about $V_{\text{DP}} \sim 0.6-0.7$ for all cases, roughly speaking. NAPL recovery improvement can be as high as about 35% for the most heterogeneous case with $V_{\text{DP}} = 0.75$. More precisely, this relative increase is about 20–35% in a co-injection mode for $K_{\text{ref},x} = 100$ or 10 mD, and of about 25–40% in an alternate injection mode for $K_{\text{ref},x} = 100$ or 10 mD, respectively.

The power law exponents reported in Figure 10b using a log-log 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_{\text{DP}}$ is observed as for NAPL recovery. However, a $V_{\text{DP}}$ threshold (much less clear-cut than for NAPL recovery) about $V_{\text{DP}} \sim 0.6$ emerges in all cases, such that for $V_{\text{DP}} \ge 0.6$ there is a significant (i) reduction of produced chemicals for co-injection up to 40–60% for $K_{\text{ref},x} = 100$ or 10 mD, and (ii) increase or little variation of produced chemicals for alternate injection, up to 5–12% for $K_{\text{ref},x} = 100$ or 10 mD.

Chemicals breakthrough time, shown in Figure 11b exhibits very significant delays for co-injection only and for $V_{\text{DP}} \ge 0.5-0.6$. It can be delayed up to 25–40% for $K_{\text{ref},x} = 100$ or 10 mD. No delay is noticeable for the alternate injection. Again, the impact of scaling $M_{\text{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 inter-layer transfers.

The differential behaviour of gas production is reported in Figure 12. There are very significant differences compared to the constant $M_{\text{ref}}$ cases as well, that mainly occur for $V_{\text{DP}} \ge 0.7$ for co-injection 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 co-injection from alternate injection for each Dykstra-Parsons 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_{\text{DP}} > 0.6$ and mainly for NAPL recovery and surfactant amounts produced.

## 7 Discussion

This work, as a follow-up of previous work [15, 52], tested a new scaling law that drives the maximum mobility reduction yielded by foam of Darcy-scale implicit-texture 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 *in-situ* saturation and composition profiles time evolution and the resulting production data, varying the degree of heterogeneity of a model stack-of-layers 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 low-permeability layers that increases the vertical sweep.

A sensitivity study was carried out to account for foam injection type (co-injection 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 Dykstra-Parsons 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_{\text{ref}}$, in relation to a constant $M_{\text{ref}}$, has been found to follow a power-law behaviour as a function of the Dykstra-Parsons coefficient. The $V_{\text{DP}}$ thresholds above which these inter-layer 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 two-dimensional 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_{\text{ref}}$, $M_{\text{ref}}$ couples when considering multi-facies 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
<figure><img src="image_10.png" /><figcaption>**Figure 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 (co-injection 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$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 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 (co-injection 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$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 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 (co-injection 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$.</figcaption></figure>

management of heterogeneity in relation to foam. This could be a decisive contribution to the predictivity of Darcy-type simulations for design purposes carried out prior to on-site 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 anti-foaming 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
<figure><img src="image_13.png" /><figcaption>Figure 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 (co-injection 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).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 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 (co-injection 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).</figcaption></figure>
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 steady-state mobility reduction implicit-texture model carries an ambiguity which remains to be clarified. The name, perhaps complicated on purpose, is indeed somewhat confusing, as this so-called steady-state 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 Darcy-scale multiphase flow model. It seems to us that the ambiguity, let us say rather lack of validation, can only be resolved by means of well-controlled transient experiments, by minimizing the model adjustable parameters of which in practice they are numerous. A flawless characterization of three-phase 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 two-phase 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 rock-fluids properties [20]. Needless to say, there is still room for improvement in dual-medium modelling that is necessary for performance purposes.

Notwithstanding, the authors hope that this work will contribute to a more reliable assessment of reservoir-scale 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 pore-scale structure of the porous medium. In this respect, it is perhaps not Darcy-type models that should be challenged first (population-balance 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 multi-scale experiments were made public.

## References

1. 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.

2. Mamun C.K., Rong J.G., Kam S.I., Liljestrand H.M., Rossen W.R. (2002) Extending foam technology form improved oil recovery to environmental remediation (Paper SPE 77557-MS), in: *SPE Annual Technical Conference and Exhibition, San Antonio, TX, USA, 29 September–2 October*.

3. 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.

4. Bertin H., Del Campo Estrada E., Atteia O. (2017) Foam placement for soil remediation, *Environ. Chem.* **14**, 338–343.

5. 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. Environ. Manage.* **243**, 187–205.

6. Flaifel H., Izadi M., Park S., Gupta I., Lee G., Kam S.I. (2020) Shallow subsurface environmental remediation by using tracer-surfacant-foam processes: history-matching and performance prediction, *Trans. Porous Media* **134**, 565–592.

7. Bouquet S., Douarche F., Roggero F., Leray S. (2020) Characterization of viscous fingering and channeling for the assessment of polymer-based heavy oil displacements, *Trans. Porous Media* **131**, 873–906.

8. *Puma is co-marketed by Beicip-Franlab and KAPPA Engineering*, Available at https://www.beicip.com/reservoir-simulation and https://www.kappaeng.com/software/rubis-puma/overview.

9. 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**, 04, 295–300.

10. Lawson J.B., Reisberg J. (1980) Alternate slugs of gas and dilute surfactant for mobility control during chemical flooding (Paper SPE-8839-MS), in: *SPE/DOE Enhanced Oil Recovery Symposium*, Tulsa, Oklahoma, April.

11. Friedmann F., Chen W.H., Gauglitz P.A. (1991) Experimental and simulation study of high-temperature foam displacement in porous media, *SPE Res. Eng.* **6**, 01, 37–45.

12. Peaceman D.W. (1977) *Fundamentals of numerical reservoir simulation, volume 6 of Developments in Petroleum Science*, Elsevier Science, Amsterdam.

13. Trangenstein J.A., Bell J.B. (1989) Mathematical structure of the black-oil model for petroleum reservoir simulation, *SIAM J. Appl. Math.* **49**, 3, 749–783.

14. Marle C.M. (1981) *Multiphase flow in porous media*, 3rd edn., Gulf Publishing Company.

15. Gassara O., Douarche F., Braconnier B., Bourbiaux B. (2020) Calibrating and scaling semi-empirical foam flow models for the assessment of foam-based EOR processes (in heterogeneous reservoirs), *Trans. Porous Media* **131**, 1, 193–221.

16. Gassara O., Douarche F., Braconnier B., Bourbiaux B. (2017) Calibrating and interpreting implicit-texture models of foam flow through porous media of different permeabilities, *J. Pet. Sci. Eng.* **159**, 588–602.
17 Gassara O., Douarche F., Braconnier B., Bourbiaux B. (2017) Equivalence between semi-empirical and population-balance foam models, *Trans. Porous Media* **120**, 3, 473–493.

18 Weaire D., Hutzler S. (1999) *The physics of foams*, Oxford University Press.

19 de Gennes P.-G., Brochard-Wyart F., Quéré D. (2004) *Capillarity and wetting phenomena: drops, bubbles, pearls, waves*, Springer-Verlag.

20 Bouquet S., Douarche F., Roggero F., Bourbiaux B. (2020) Foam processes in naturally fractured carbonate oil-wet reservoirs: technical and economic analysis and optimization, *J. Pet. Sci. Eng* **190**, 107111.

21 Aziz K., Settari A. (1985) *Petroleum reservoir simulation*, Elsevier.

22 Fayers F.J. (1989) Extension of Stone's method 1 and conditions for real characteristics in three-phase flow, *SPE Res. Eng.* **4**, 4, 437–445.

23 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.

24 Bretherton F.P. (1961) The motion of long bubbles in tubes, *J. Fluid Mech.* **10**, 2, 166–188.

25 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.

26 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 foam-oil interactions and an improved bubble population-balance model (Paper SPE 191 564), in: *2018 SPE Annual Technical Conference and Exhibition*, Dallas, TX, USA, 24–26 September.

27 Chen Q., Kovscek A.R., Gerritsen M. (2010) Modeling foam displacement with the local-equilibrium approximation: theory and experimental verification, *SPE J.* **15**, 01, 171–183.

28 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.

29 Patzek T.W. (1988) Description of foam flow in porous media by the population balance method, in: H.S. Duane (ed.), *Surfactant-based mobility control, American Chemical Society Symposium Series*.

30 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 snap-off, *SPE Reserv. Eng.* **3**, 03, 884–892.

31 Kam S.I., Nguyen Q.P., Li Q., Rossen W.R. (2007) Dynamic simulations with an improved model for foam generation, *SPE J.* **12**, 01, 35–48.

32 Kam S.I. (2008) Improved mechanistic foam simulation with foam catastrophe theory, *Colloids Surf. A Physicochem. Eng. Asp.* **318**, 1–3, 62–77.

33 Kovscek A.R., Bertin H.J. (2003) Foam mobility in heterogeneous porous media, *Trans. Porous Media* **52**, 1, 37–49.

34 Farajzadeh R., Eftekhari A. (2017) Effect of foam on liquid phase mobility in porous media, *Sci. Rep.* **43870**, 3011–3018.

35 Farajzadeh R., Lotfollahi M., Eftekhari A.A., Rossen W.R., Hirasaki G.J.H. (2015) Effect of permeability on implicit-texture foam model parameters and the limiting capillary pressure, *Energy Fuels* **29**, 5, 3011–3018.

36 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.

37 Kapetas L., Vincent-Bonnieu S., Farajzadeh R., Eftekhari A. A., Mohd-Shafian S.R., Kamarul Bahrim R.Z., Rossen W.R. (2015) Effect of permeability on foam-model parameters – an integrated approach from coreflight experiments through to foam diversion calculations, in *18th European Improved Oil Recovery Symposium*, Dresden, Germany, 14–16 April.

38 Jian G., Zhang L., Da C., Puerto M., Johnston K.P., Biswal S.L., Hirasaki G.J. (2019) Evaluating the transport behavior of CO₂ foam in the presence of crude oil under high-temperature and high-salinity conditions for carbonate reservoirs, *Energy Fuels* **33**, 6038–6047.

39 Bergeron V., Radke C.J. (1992) Equilibrium measurements of oscillatory disjoining pressures in aqueous foam films, *Langmuir* **8**, 12, 3020–3026.

40 Soulat A., Douarche F., Flauraud E. (2020) A modified well index to account for shear-thinning behavior in foam EOR simulation, *J. Pet. Sci. Eng.* **191**, 107146.

41 Peaceman D. (1978) Interpretation of well-block pressures in numerical reservoir simulation, *SPE J.* **18**, 3, 182–194.

42 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.

43 Brooks R.H., Corey A.T. (1966) Properties of porous media affecting fluid flow, *J. Irrig. Drain. Div.* **92**, 2, 61–88.

44 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.

45 de Figueiredo L., Grana D., Le Ravalec M. (2020) Revisited formulation and applications of FFT moving average, *Math. Geosci.* **52**, 6, 810–816.

46 Le Ravalec M., Nøtinger B., Hu L.Y. (2000) The fft moving average (FFT-MA) generator: an efficient numerical method for generating and conditioning Gaussian simulations, *Mathe. Geol.* **32**, 6, 701–723.

47 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.

48 Craig F.F. (1993) *The reservoir engineering aspects of waterflooding, volume 3 of SPE Monograph Series*, Society of Petroleum Engineers.

49 Willhite G.P. (1986) *Waterflooding, volume 3 of SPE Textbook Series*, Society of Petroleum Engineers.

50 Lake L.W. (1989) *Enhanced oil recovery*, Prentice Hall.

51 Green D.W., Willhite G.P. (1998) *Enhanced oil recovery, volume 6 of SPE Textbook Series*, Society of Petroleum Engineers.

52 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.

53 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).

54 Weber K.J., van Geuns L.C. (1990) Framework for constructing clastic reservoir simulation models, *J. Pet. Technol.* **42**, 1248–1297.
55 Leverett M.C. (1941) Capillary behavior in porous solids (Paper SPE-941152-G), *Trans. AIME* **142**, 152–169.

# 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 water-gas co-injection steady-state measurements of differential pressure across core-plugs, $\Delta p_f/\Delta p$, with and without foam, as measured in the laboratory, given that $\Delta 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 steady-state $\Delta p$ is a gas injection, in which case it would be tempting to consider that the ratio $\Delta p_f/\Delta p$ corresponds to $M_{ref}$, where $\Delta p_f$ denotes the steady-state pressure drop obtained under a co-injection 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 steady-state $\Delta p_f$ obtained under foam injection (most often co-injection of surfactant solution and gas at a prescribed quality close to the optimal one), not to the $\Delta p$ obtained under the equivalent co-injection 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} = \Delta p_f/\Delta p$ where $\Delta p_f$ and $\Delta p$ are respectively the steady-state pressure differences measured across core-plugs flooded by a surfactant solution/gas co-injection yielding foam and by the corresponding co-injection 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 co-injection without surfactant, at the same imposed total flow rate and quality, writes for the gas phase

$$f_g U = \frac{K k_{rg}^f(S_g^f) \Delta p_f}{\mu_g L} \quad (\text{foam}), \quad (A1)$$

$$f_g U = \frac{K k_{rg}(S_g) \Delta p}{\mu_g L} \quad (\text{co-injection}),$$

where $U$ is the total velocity. Indeed, not only the relative permeability to gas is modified in the presence of foam ($k_{rg} \to k_{rg}^f$), but the saturation of the displacing gas front is also ($S_g \to S_g^f$). This is precisely what a Buckley-Leverett analytical calculation would predict [14]. From relations (A1) we get

$$M_{\text{lab}} = \frac{\Delta p_f}{\Delta p} = \frac{k_{rg}(S_g)}{k_{rg}^f(S_g^f)}. \quad (A2)$$

Moreover, in the presence of foam, we have

$$k_{rg}^f(S_g^f) = FM \cdot k_{rg}(S_g^f), \quad (A3)$$

$$k_{rw}^f(S_g^f) = k_{rw}(S_g^f),$$

according to (4). Under optimal reference conditions, $FM$ given by (5) writes, using (A3)

$$FM = \frac{1}{M_{\text{ref}}} = \frac{k_{rg}^f(S_g^f)}{k_{rg}(S_g^f)}, \quad (A4)$$

i.e., using (A2)

$$M_{\text{ref}} = \frac{k_{rg}(S_g^f)}{k_{rg}^f(S_g^f)} = M_{\text{lab}} \frac{k_{rg}(S_g^f)}{k_{rg}(S_g^f)}. \quad (A5)$$

The determination of $M_{ref}$ is thus reduced to that of the saturations $S_g^f$ and $S_g$, $S_g^f$ depending on $S_g$.

*Determination of $S_g^f$ under foam injection.* Darcy’s law, for foam flow and for the corresponding co-injection without surfactant, writes for the water phase

$$(1 - f_g) U = \frac{K k_{rw}(S_g^f) \Delta p_f}{\mu_w L} \quad (\text{foam}), \quad (A6)$$

$$(1 - f_g) U = \frac{K k_{rw}(S_g) \Delta p}{\mu_w L} \quad (\text{co-injection}),$$

i.e. by forming the ratio

$$S_g^f = k_{rw}^{-1} \left( \frac{k_{rw}(S_g)}{M_{\text{lab}}} \right), \quad (A7)$$

where $k_{rw}^{-1}$ denotes the reciprocal of $k_{rw}$.

*Determination of $S_g$ under co-injection.* The gas saturation $S_g$ verifies the relation

$$S_g = \left( \frac{k_{rg}}{k_{rw}} \right)^{-1} \left( \frac{f_g}{1 - f_g} \frac{\mu_g}{\mu_w} \right). \quad (A8)$$

A similar calculation gives for a reference water injection

$$\left\{ \begin{array}{l} M_{\text{lab}} = \frac{\mu_g}{\mu_w} \frac{f_g}{k_{rg}^f(S_g^f)} \\ M_{\text{ref}} = M_{\text{lab}} \frac{\mu_w}{\mu_g} \frac{k_{rg}(S_g^f)}{f_g} \\ S_g^f = k_{rw}^{-1} \left( \frac{1 - f_g}{M_{\text{lab}}} \right) \end{array} \right. \quad (A9)$$

and for a gas injection

$$\left\{ \begin{array}{l} M_{\text{lab}} = \frac{f_g}{k_{rg}^f(S_g^f)} \\ M_{\text{ref}} = M_{\text{lab}} \frac{k_{rg}(S_g^f)}{f_g} . \\ S_g^f = k_{rw}^{-1} \left( \frac{\mu_w}{\mu_g} \frac{1 - f_g}{M_{\text{lab}}} \right) \end{array} \right. \quad (A10)$$

*Orders of magnitude.* Power-law relative permeabilities (16) and (17) rewrite for a water and gas two-phase flow in the absence of NAPL, setting $S_{orw} = S_{org}$

$$k_{rg}(S_G) = \kappa_g(S_G)^{n_g} \quad (A11)$$

$$k_{rw}(S_G) = \kappa_w(1 - S_G)^{n_w}$$
<table><caption>**Table A1.** Examples of $M_{\text{lab}}$ to $M_{\text{ref}}$ conversion for a reference injection which is either a co-injection of water and gas or a water injection (the gas considered is CO$_2$ with $T = 100$ °C and $p = 100$ bar).</caption><thead><tr><th>Reference injection</th><th>Co-injection</th><th>Water</th></tr></thead><tbody><tr><td>$f_g$</td><td>0.8</td><td>—</td></tr><tr><td>$\mu_g$ (cP)</td><td>0.02</td><td>—</td></tr><tr><td>$\mu_w$ (cP)</td><td>0.28</td><td>—</td></tr><tr><td>$S_{gr}$</td><td>0</td><td>—</td></tr><tr><td>$S_{wi}$</td><td>0.15</td><td>—</td></tr><tr><td>$\kappa_g$</td><td>0.2</td><td>—</td></tr><tr><td>$\kappa_w$</td><td>0.2</td><td>—</td></tr><tr><td>$n_g$</td><td>2</td><td>—</td></tr><tr><td>$n_w$</td><td>4</td><td>—</td></tr><tr><td>$M_{\text{lab}}$</td><td>71</td><td>196</td></tr><tr><td>$M_{\text{ref}}$</td><td>500</td><td>1581</td></tr></tbody></table>

with $S_G = (S_g - S_{gr})/(1 - S_{wi} - S_{gr})$. Having measured $M_{\text{lab}}$, to determine $M_{\text{ref}}$ we compute $S_g$ according to (A8), then $S_g^f$ following (A7) and finally $M_{\text{ref}}$ from (A5). Some examples are given in Table A1. The $M_{\text{ref}}$ values of 500 and 1581 assigned to $K_{\text{ref}} = 100$ mD in this work correspond to $M_{\text{lab}} = 71$ and 196 if the reference injection is a co-injection of water and gas, and to $M_{\text{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].

## A.2 Comparison of in-situ profiles with and without $M_{\text{ref}}$ scaling
<figure><img src="image_15.png" /><figcaption>**Figure A1.** Co-injection case after 1 PV of water injection prior to surfactant solution/gas co-injection, with $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 100$ mD, $V_{\text{DP}} = 0.75$, $K_z/K_x = 0.1$ and $L_z = 10$ m. Left: constant $M_{\text{ref},x}$. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$. Right: comparison of horizontally averaged profiles without ("ref") and with ("mod") $M_{\text{ref},x}$ scaling.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure A2.** Co-injection case after 0.25 PV of surfactant solution injection, with $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 100$ mD, $V_{\text{DP}} = 0.75$, $K_z/K_x = 0.1$ and $L_z = 10$ m. Left: constant $M_{\text{ref},x}$. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$. Right: comparison of horizontally averaged profiles without ("ref") and with ("mod") $M_{\text{ref},x}$ scaling.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure A3.** Co-injection case final state, with $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 100$ mD, $V_{\text{DP}} = 0.75$, $K_z/K_x = 0.1$ and $L_z = 10$ m. Left: constant $M_{\text{ref},x}$. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$. Right: comparison of horizontally averaged profiles without ("ref") and with ("mod") $M_{\text{ref},x}$ scaling.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure A4.** Co-injection case after 0.25 PV of surfactant solution injection, with $M_{\text{ref},x} = 1581$, $K_{\text{ref},x} = 100$ mD (or equivalently $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 10$ mD), $V_{\text{DP}} = 0.75$, $K_z/K_x = 0.1$ and $L_z = 10$ m. Left: constant $M_{\text{ref},x}$. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$. Right: comparison of horizontally averaged profiles without ("ref") and with ("mod") $M_{\text{ref},x}$ scaling.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

**Figure A5.** Co-injection case final state, with $M_{\text{ref},x} = 1581$, $K_{\text{ref},x} = 100$ mD (or equivalently $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 10$ mD), $V_{\text{DP}} = 0.75$, $K_z/K_x = 0.1$ and $L_z = 10$ m. Left: constant $M_{\text{ref},x}$. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$. Right: comparison of horizontally averaged profiles without ("ref") and with ("mod") $M_{\text{ref},x}$ scaling.
<figure><img src="image_20.png" /><figcaption>**Figure A6.** Figure A6 Co-injection case after 0.25 PV of surfactant solution injection, with $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 100$ mD, $V_{\text{DP}} = 0.75$, $K_z/K_x = 1$ and $L_z = 10$ m. Left: constant $M_{\text{ref},x}$. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$. Right: comparison of horizontally averaged profiles without (‘ref’) and with (‘mod’) $M_{\text{ref},x}$ scaling.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure A7.** Co-injection case after 0.25 PV of surfactant solution injection, with $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 100$ mD, $V_{\text{DP}} = 0.75$ and $L_z = 10$ m. Left: scaled $M_{\text{ref},x}$ with $K_z/K_x = 1$. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$ with $K_z/K_x = 0.1$. Right: comparison of horizontally averaged profiles without ("ref") and with ("mod") $M_{\text{ref},x}$ scaling ("ani" and "iso" stand for $K_z/K_x = 0.1$ and 1, respectively).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure A8.** Co-injection case after 0.25 PV of surfactant solution injection, with $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 100$ mD, $V_{\text{DP}} = 0.75$ and $K_z/K_x = 0.1$. Left: scaled $M_{\text{ref},x}$ with $L_z = 10$ m. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$ with $L_z = 5$ m. Right: comparison of horizontally averaged profiles without ("ref") and with ("mod") $M_{\text{ref},x}$ scaling ("$L_z$" and "$L_z/2$" stand for $L_z = 10$ and 5 m, respectively).</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure A9.** Alternate injection case after 0.25 PV of surfactant solution injection, with $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 100$ mD, $V_{\text{DP}} = 0.75$, and $K_z/K_x = 0.1$ and $L_z = 10$ m. Left: constant $M_{\text{ref},x}$. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$. Right: comparison of horizontally averaged profiles without ("ref") and with ("mod") $M_{\text{ref},x}$ scaling.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure A10.** Alternate injection case post surfactant solution injection after approximately 0.9 PV of gas injection, with $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 100$ mD, $V_{\text{DP}} = 0.75$, $K_z/K_x = 0.1$ and $L_z = 10$ m. Left: constant $M_{\text{ref},x}$. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$. Right: comparison of horizontally averaged profiles without ("ref") and with ("mod") $M_{\text{ref},x}$ scaling.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure A11.** Alternate injection case post surfactant solution injection after approximately 1.8 PV of gas injection, with $M_{\text{ref},x} = 500$, $K_{\text{ref},x} = 100$ mD, $V_{\text{DP}} = 0.75$, $K_z/K_x = 0.1$ and $L_z = 10$ m. Left: constant $M_{\text{ref},x}$. Solid and dashed lines stand for horizontally and vertically averaged profiles, respectively. Center: scaled $M_{\text{ref},x}$. Right: comparison of horizontally averaged profiles without ("ref") and with ("mod") $M_{\text{ref},x}$ scaling.</figcaption></figure>