O Gassara
email: omar.gassara@ifpen.fr

F Douarche

B Braconnier
email: benjamin.braconnier@ifpen.fr

B Bourbiaux
email: bernard.bourbiaux@ifpen.fr

Equivalence between semi-empirical and population-balance foam models

Keywords: Multiphase flow, Porous media, Foam, Models, Reservoir simulation, Reservoir engineering, Enhanced oil recovery

Models for simulating foam-based displacements in enhanced oil recovery (EOR) processes fall into two categories: population-balance (PB) models that derive explicitly foam texture, or bubble size, evolution in porous media from pore-level mechanisms related to lamellas generation and coalescence, and semi-empirical (SE) models that account implicitly for foam texture effects through a gas mobility reduction factor that depends on fluid saturation, interstitial velocity, surfactant concentration, and other factors. This mobility reduction factor has to be calibrated from a large number of experiments on a case by case basis in order to match the physical effect of each considered parameter on foam behavior. This paper develops a method for identifying the SE models from the physics of foams as derived from PB models at local equilibrium (LE). The identification of both foam flow models leads to a method for calibrating SE models from the PB model translation of foam flow data. Application to a set of foam quality-scan experiments at fixed total flow rate shows that the SE and PB models at LE match equally well the measurements and generate almost the same results in both the so-called high-and low-quality regimes. We demonstrate that the two approaches are equivalent at local equilibrium and differ only in the way in which the complex dynamic mechanisms of lamellas are handled. This physical approach of foam flow could circumvent some difficulties in the direct calibration of SE models from foam mobility (or apparent viscosity) data.

Introduction

In view of the recent oil production challenges resulting from rising oil demand and declining new fields discoveries, enhanced oil recovery (EOR) methods are increasingly becoming important. They follow primary production by depletion and secondary recovery methods by water or gas injection. Gas injection processes suffer from poor volumetric sweep efficiency due to both channeling and gravity override. Foam has revealed great potential to overcome these detrimental effects by reducing gas mobility, and consequently, it proves itself as a promising EOR method [START_REF] Bond | Gas drive oil recovery[END_REF][START_REF] Holm | Foam injection test in the Siggins field[END_REF][START_REF] Jonas | Evaluation of a CO 2 foam field trial: Rangely Weber Sand Unit[END_REF][START_REF] Watts | Field test of CO 2 mobility control at Rock Creek[END_REF].

Foam consists in a succession of gas bubbles separated by thin liquid films called lamellas. Foam is usually characterized by its quality and texture. The foam quality is the gas volumetric fraction in the bulk foam, or the gas fractional rate when injected in a porous medium, and the foam texture is the number of lamellas per unit volume of foam. The gas phase flowing as a dispersed phase within a foam is much less mobile than the same gas flowing as a continuous phase. That mobility reduction can be seen as the result of an increased effective viscosity of gas, of a decreased relative permeability of gas [START_REF] Falls | Development of a mechanistic foam simulator: the population balance and generation by snap-off[END_REF], of a high trapping of gas [START_REF] Kovscek | Fundamentals of foam transport in porous media. in Foams: fundamentals and applications in the petroleum industry[END_REF], or also as a conjunction of several such phenomena. Increased effective viscosity of gas flowing as a foam can be explained thanks to a microscopic approach of foam flowing as a train of gas bubbles separated by lamellas [START_REF] Hirasaki | Mechanisms of foam flow in porous media: apparent viscosity in smooth capillaries[END_REF], which are responsible for an additional resistance to gas flow according to Bretherton's study of bubble flow in a capillary [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF]. Therefore, the effective gas viscosity in presence of foam depends strongly on foam texture: the finer the texture, the larger the number of lamellas, hence the higher resistance to the flow of the gas phase. However, foam flow properties within a porous medium are not controlled by the foam texture alone, because porous structure and pore thresholds have a determinant impact on the flow behavior of a given foam. Furthermore, foam texture can change within porous media as the result of multiple pore-level mechanisms resulting in foam bubbles creation and coalescence [START_REF] Ransohoff | Mechanisms of foam generation in glass bead packs[END_REF][START_REF] Khatib | Effects of capillary pressure on coalescence and phase mobilities in foams flowing through porous media[END_REF][START_REF] Kovscek | Fundamentals of foam transport in porous media. in Foams: fundamentals and applications in the petroleum industry[END_REF][START_REF] Chambers | Capillary phenomena in foam flow through porous media[END_REF].

Experimental evidence of lamellas coalescence in foam models led Khatib et al [START_REF] Khatib | Effects of capillary pressure on coalescence and phase mobilities in foams flowing through porous media[END_REF] to introduce the limiting capillary pressure P * c concept. Above that threshold, foam films become unstable, and foam becomes coarser and coarser as a result of bubbles coalescence (cf Fig. 1). The magnitude of P * c varies with surfactant concentration, electrolyte concentration, gas flow rate and porous medium permeability [START_REF] Khatib | Effects of capillary pressure on coalescence and phase mobilities in foams flowing through porous media[END_REF][START_REF] Farajzadeh | Effect of permeability on implicit-texture foam model parameters and the limiting capillary pressure[END_REF].

The behavior of foam flowing at a given velocity within a porous medium exhibits two distinct flow regimes depending on foam quality [START_REF] Osterloh | Effects of gas and liquid velocity on steady-state foam flow at high temperature[END_REF]: a low-quality regime such that the pressure gradient, or the apparent foam viscosity, increases with foam quality until an optimal foam quality f * g is obtained, and a high-quality regime such that the pressure gradient and the apparent foam viscosity decrease.

A variety of models has been proposed for modeling foam flow in porous media, which can be classified into two main groups. The first group includes population-balance (PB) models that account for pore-level mechanisms related to lamella creation and coalescence in order to track the dynamical behavior of foam texture in porous media along with its effect on gas mobility [START_REF] Falls | Development of a mechanistic foam simulator: the population balance and generation by snap-off[END_REF][START_REF] Friedmann | Experimental and simulation study of high-temperature foam displacement in porous media[END_REF][START_REF] Ettinger | Influence of texture on steady foam flow in Berea sandstone[END_REF][START_REF] Kovscek | Foam mobility in heterogeneous porous media[END_REF]. Alternatively, a local steady state version of PB models [START_REF] Ettinger | Influence of texture on steady foam flow in Berea sandstone[END_REF][START_REF] Kovscek | Foam mobility in heterogeneous porous media[END_REF][START_REF] Chen | Modeling foam displacement with the local-equilibrium approximation: theory and experimental verification[END_REF] comes out by equating the lamella creation and coalescence rates; solving this equation yields the foam texture. Semi-empirical (SE) models, that constitute the second group, assume that an equilibrium state is attained instantaneously such that the foam texture effects can be represented implicitly through a gas mobility reduction factor that depends on fluid saturation, interstitial velocity, surfactant concentration and other factors [START_REF] Martinsen | Foam assisted water alternating gas (FAWAG) process on Snorre[END_REF][START_REF] Cheng | Simulating foam processes at high and low foam qualities[END_REF]. SE models incorporate implicitly the limiting capillary pressure with a steep transition function which represents the coalescence effects on the mobility reduction factor. This mobility reduction factor has to be calibrated from a large number of experiments on a case by case basis in order to match the physical effect of each considered parameter on foam behavior.

In the laboratory, the foam flow behavior in a porous medium is generally characterized by its apparent viscosity µ app f . Considering foam as a single homogeneous phase in a porous medium whose boundaries are subjected to a pressure difference ∆ P, one can apply Darcy law and define the foam apparent viscosity as µ

app f = k |u| ∆ P L ,
where k is the average porous media permeability, L the porous medium length across which the pressure difference is applied, and |u| the foam velocity. Procedures to fit SE models to foam quality-scan experiments at fixed total flow rate have been successfully provided by several authors [START_REF] Cheng | Simulating foam processes at high and low foam qualities[END_REF][START_REF] Ma | Nonuniqueness, numerical artifacts, and parameter sensitivity in simulating steady-state and transient foam flow through porous media[END_REF][START_REF] Boeije | Fitting foam-simulation-model parameters to data: I. coinjection of gas and liquid[END_REF][START_REF] Farajzadeh | Effect of permeability on implicit-texture foam model parameters and the limiting capillary pressure[END_REF]. They are discussed later on, and compared to the present paper method.

This paper develops a method for identifying the SE models from the physics of foams as derived from PB models. The motivation of this work is to strengthen further the physical basis of SE models and constrain the model calibration by reproducing as closely as possible the foam texture effects on gas mobility, as suggested by PB models. The addressed items are covered in the following, namely:

• Summary of PB and SE models features and differences, • Equivalence between SE and PB models at local equilibrium, • Validation of the proposed SE models identification method against experimental data, • Comparison between the results of this work with others found in the literature.

Foam modeling approaches

In this section, we detail the two main modeling categories that describe foam in the framework of a reservoir simulator, namely the population-balance (PB) models and the semiempirical (SE) models. These models extend the usual three phase reservoir simulators based on the generalized Darcy law. Thus, for these models, foam is not modeled as an additional phase: rather, it is still described as a mixture of two water and gas phases, where the gas phase mobility is modified and the liquid phase mobility remains unchanged [START_REF] Bernard | Effect of foam on trapped gas saturation and on permeability of porous media to water[END_REF][START_REF] Lawson | Alternate slugs of gas and dilute surfactant for mobility control during chemical flooding[END_REF][START_REF] Friedmann | Experimental and simulation study of high-temperature foam displacement in porous media[END_REF]. In addition, both models require to model a surfactant which is needed for the foam or lamellas formation and stability. Before underlining PB and SE models specific features regarding the gas mobility, we present the material balance equations for the water, hydrocarbon and surfactant components that underlie both types of model.

Hydrocarbons, surfactant and water mass balance equations

We consider a model for a three-phase flow in a porous medium in the presence of foam. We distinguish three phases: an aqueous phase w, an oil phase o and a gas phase g. This flow is modified by the presence of foam. Modeling foam requires the presence of a surfactant, which is transported by the water phase, and which is described by an additional mass balance equation. The surfactant is either mobile or adsorbed on the rock. When present in the water phase, we assume that the surfactant mass does not affect the water mass.

It has been shown [START_REF] Bernard | Effect of foam on trapped gas saturation and on permeability of porous media to water[END_REF][START_REF] Lawson | Alternate slugs of gas and dilute surfactant for mobility control during chemical flooding[END_REF][START_REF] Friedmann | Experimental and simulation study of high-temperature foam displacement in porous media[END_REF] that the transport of liquid is not affected by the presence of foam. On the opposite, the gas velocity is significantly reduced by the presence of foam. Thus, to describe the water and hydrocarbon phases, we consider a black-oil model [START_REF] Peaceman | Fundamentals of numerical reservoir simulation[END_REF][START_REF] Trangenstein | Mathematical structure of the black-oil model for petroleum reservoir simulation[END_REF] where the gas phase involves a modified velocity which will be denoted u f g . The mass conservation equations read

           ∂ t (Φρ w S w ) + ∇ • (ρ w u w ) = q w , ∂ t (Φρ w S w C s w + (1 -Φ)ρ r C s r ) + ∇ • (ρ w u w C s w ) = q w C s w , ∂ t (ΦC h ρ o S o ) + ∇ • (ρ o C h u o ) = C h q o , ∂ t (Φρ g S g + ΦC v ρ o S o ) + ∇ • (ρ o C v u o + ρ g u f g ) = q g +C v q o , (1)
where Φ is the rock porosity. For each phase denoted ψ = w, o, g, S ψ is the saturation, ρ ψ the mass density and q ψ the source/sink term per unit volume of porous medium. The gas phase contains a single volatile component denoted v, whereas the oil 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 that phase). The equilibrium constant K v = 1 X v , with X v the molar fraction of volatile component in the oil phase, is defined from C v and C h . K v is a function of pressure. C s w stands for the flowing surfactant mass fraction in the water phase and C s r for the adsorbed surfactant mass fraction on the rock with ρ r the rock mass density. Mobile and adsorbed surfactant mass fractions are related with an adsorption law such as the Langmuir isotherm [START_REF] Lake | Enhanced oil recovery[END_REF][START_REF] Green | of SPE Textbook Series[END_REF].

Under creeping (i.e. low-velocity) flow conditions, the pure phase velocities in permeable porous media are governed by the generalized Darcy law:

u ψ = - kk rψ µ ψ ∇P ψ -ρ ψ g , ( 2
)
where k is the rock permeability, µ ψ the pure phase viscosity, P ψ the pressure of the phase ψ and g the gravity. k rψ is the relative permeability for the pure phase ψ, i.e. without lamellas. We suppose the relative permeabilities and capillary pressure functions are known, with given saturation end points. In our specific context of foam modeling, we will only consider the gas and water phases, thus S g ∈ [S gr , 1 -S wi ] where S gr is the residual gas saturation and S wi the irreducible water saturation.

In order to simplify the notations in the following, we introduce the phase mobility λ ψ = k rψ /µ ψ and the interstitial velocity v ψ = u ψ /(ΦS ψ ) (see [START_REF] Marle | Multiphase flow in porous media[END_REF] for more details). We also introduce the phase fractional flow f ψ = |u ψ |/|u| where the total velocity u = ∑ ψ∈{w,o,g} u ψ is the sum of the phase velocities.

The flow of gas in the presence of foam is discussed later on as it is modeled differently whether a population-balance or an semi-empirical modeling approach is used. These models are based on Darcy type laws which are extended to obtain a modified gas velocity so that

u f g = u PB g for PB models, u SE g for SE models, (3)
where u PB g stands for the gas velocity in the presence of foam according to PB models and u SE g is the gas velocity in the presence of foam for the SE models. These two gas velocities models are detailed in the two next sections.

Population-balance lamellas models

Foam mobility is strongly related to its texture which is a key variable in foam modeling 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 modeled and a lamella population-balance equation is considered. The impact of foam texture on gas mobility is modeled through an effective gas viscosity and a gas trapping for more elaborated models [START_REF] Friedmann | Experimental and simulation study of high-temperature foam displacement in porous media[END_REF][START_REF] Chen | Modeling foam displacement with the local-equilibrium approximation: theory and experimental verification[END_REF][START_REF] Kovscek | Fundamentals of foam transport in porous media. in Foams: fundamentals and applications in the petroleum industry[END_REF]. In this study, we focus on viscosity-type PB models and do not consider gas trapping.

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 [START_REF] Patzek | Description of foam flow in porous media by the population balance method[END_REF][START_REF] Falls | Development of a mechanistic foam simulator: the population balance and generation by snap-off[END_REF]. The lamellas population balance equation reads

∂ t (ΦS g n f ) + ∇ • (u PB g n f ) = ΦS g (r g -r c ) + q f , (4)
where n f 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). As mentioned earlier, the gas velocity is significantly reduced when foam develops. The viscosity-type PB models extend the generalized Darcy equation for the gas phase as follows:

u PB g = - kk rg µ PB g (∇P g -ρ g g) ,
where µ PB g is the effective gas viscosity when flowing as a succession of bubbles. The rheology of PB models of foam flow is classically based on Bretherton's flow model of a single bubble within a capillary [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF], that was later extended by Hirasaki and co-workers to a train of bubbles [START_REF] Hirasaki | Mechanisms of foam flow in porous media: apparent viscosity in smooth capillaries[END_REF] thus leading to the following expression for the effective gas viscosity:

µ PB g = µ g + c f n f |v PB g | 1/3 , (5)
where v PB g is the interstitial gas velocity in the presence of foam, deduced from u PB g . c f is a constant depending on the surfactant concentration and the permeability of the porous medium [START_REF] Kovscek | Foam mobility in heterogeneous porous media[END_REF]. r g and r c are the rates of lamellas creation and coalescence.

Different formulas have been proposed in the PB models of the literature to model these rates [START_REF] Kovscek | Fundamentals of foam transport in porous media. in Foams: fundamentals and applications in the petroleum industry[END_REF][START_REF] Kam | Dynamic simulations with an improved model for foam generation[END_REF][START_REF] Kam | Improved mechanistic foam simulation with foam catastrophe theory[END_REF]. For the purpose of our paper, we propose to focus on the Kam et al [START_REF] Kam | Dynamic simulations with an improved model for foam generation[END_REF] model for which the creation and coalescence rates read:

r g = c g S w ∇P g m , r c = c c n f (S w -S * -PB w ) n , (6)
where c g is the generation rate coefficient, c c the coalescence rate coefficient, n a coalescence exponent, m a model parameter, and where the foam texture n f cannot exceed a maximum value n max f . The formulations of these rates need to be commented as they involve porous medium control on the typical foam behavior that can be observed in core flow experiments. Regarding lamellas generation, Kam et al [START_REF] Kam | Dynamic simulations with an improved model for foam generation[END_REF] suggested that the rate of lamellas creation in porous media is proportional to the water saturation and to a power-law expression of the gas-phase pressure gradient. Specifically, foam generation is easier at higher water saturation because more liquid lenses can be generated within the pore network. Regarding lamellas coalescence, the rate r c is expressed as a function of the foam texture n f and the limiting water saturation S * -PB w , corresponding to the limiting capillary pressure P * c . The expression of r c predicts that lamellas destruction diverges toward infinity as S w approaches S * -PB w . From a more mathematical point of view, let us observe that the local steady state version of the population-balance model is obtained by assuming local steady state and instantaneous equilibrium between the rates of foam generation and coalescence. This equality between r g and r c yields the foam texture, as will be shown in Section 4.

As regards the foam texture in a porous medium, it is acknowledged that the size of foam bubbles (that is the inverse of n f ) is closely depending on the porous structure [START_REF] Falls | Development of a mechanistic foam simulator: the population balance and generation by snap-off[END_REF][START_REF] Friedmann | Experimental and simulation study of high-temperature foam displacement in porous media[END_REF]. In first approximation, it can be estimated as the volume of spherical bubbles with a radius r in the order of characteristic pore radius, the latter being a function of porous medium flow properties [START_REF] Kovscek | Foam mobility in heterogeneous porous media[END_REF]. Following Alvarez et al's working hypothesis [START_REF] Alvarez | Unified model for steady-state foam behavior at high and low foam qualities[END_REF], Kam et al [START_REF] Kam | Dynamic simulations with an improved model for foam generation[END_REF] assume that the size of individual foam bubbles is constant in the low-quality regime once the conditions of foam generation (velocity) are met, then increases in the high-quality regime because of coalescence. To end with, according to Kam et al model, the maximum foam texture n max f is obtained in the low-quality regime and can be estimated from a characteristic pore radius. However, that question of foam texture remains controversial because the mechanisms of lamellas generation are complex [START_REF] Kovscek | Fundamentals of foam transport in porous media. in Foams: fundamentals and applications in the petroleum industry[END_REF] and hardly tractable in the context of natural porous media, but such a discussion is beyond the scope of present work.

Overall, expressions (6) incorporate the qualitative physical trends of foam texture dependence on the porous medium as discussed above, but they still have to be calibrated from experiments, such as the ones from Alvarez et al [START_REF] Alvarez | Unified model for steady-state foam behavior at high and low foam qualities[END_REF] and Moradi-Araghi et al [START_REF] Moradi-Araghi | Laboratory evaluation of surfactants for CO 2 -foam applications at the South Cowden Unit[END_REF].

Semi-empirical foam models

Most semi-empirical (SE) foam flow simulators are based on the same modeling approach that consists in applying a gas mobility reduction factor when foam is present. A SE model implemented in IFPEN simulator PumaFlow [START_REF] Braconnier | Efficient scheme for chemical flooding[END_REF] is considered herein for the only purpose of demonstrating the identification methodology of that paper. Herein, a two-phase gaswater fluid system is considered. Surfactant is transported by the water phase and assumed not to adsorb on the rock. Hence, the mass conservation equations system (1) reads

     ∂ t (Φρ w S w ) + ∇ • (ρ w u w ) = q w , ∂ t (Φρ w S w C s w ) + ∇ • (ρ w u w C s w ) = q w C s w , ∂ t (Φρ g S g ) + ∇ • (ρ g u f g ) = q g . (7)
The main features of that SE model are described hereafter. The gas mobility λ SE g of SE model is scaled by a multi-parameter interpolation function FM assigned to the relative permeability while the gas viscosity is assumed unchanged whether foam is present or not:

u SE g = -kλ SE g (∇P g -ρ g g) with λ SE g = k SE rg µ g = FM • k rg µ g , (8)
where k SE rg is the modified gas relative permeability for the SE model in the presence of foam and k rg is the conventional relative permeability of gas. 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 = 1 1 + (M ref -1) ∏ 4 i=1 F i , (9)
where M 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 saturation, oil saturation, and gas velocity, or equivalently, the gas capillary number. In this paper, we study the effect of water saturation and gas velocity on foam performance. Thus, we focus on the so called dry-out function F 2 and shear-thinning function F 4 (F 1 and F 3 functions account for surfactant concentration and oil saturation effect on foam performance, and are not developed herein). In our model,

F 2 reads F 2 (S w ) = 1 2 + 1 π arctan Θ S w -S * -SE w , (10)
where the dimensionless constant Θ governs the sharpness of the transition from the lowquality regime to the high-quality regime when water saturation decreases below S * -SE w . A very high value of Θ (several thousands) leads to a steep transition and in this case foam coalescence occurs at the given water saturation S * -SE w , whereas a low value (less than 100 or 10) means a smooth transition in the vicinity of S * -SE w . Furthermore, FM includes a shear-thinning function F 4 expressing the non-Newtonian behavior of foam flow in porous media: at fixed foam parameters (surfactant concentration, foam texture and fluids saturations), foam viscosity decreases with gas velocity, and subsequently, foam mobility increases [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF][START_REF] Hirasaki | Mechanisms of foam flow in porous media: apparent viscosity in smooth capillaries[END_REF][START_REF] Falls | Development of a mechanistic foam simulator: the population balance and generation by snap-off[END_REF][START_REF] Friedmann | Experimental and simulation study of high-temperature foam displacement in porous media[END_REF]. In SE models, this dependence is expressed via the local gas capillary number N cg which represents the relative effect of viscous forces versus capillary ones. The shear-thinning function is defined as

F 4 (N cg ) = N ref cg N cg e c with N cg = µ g |u SE g | Φσ S g , (11)
where N ref cg is a reference capillary number for which M ref is measured, e c an exponent representing the extent of shear-thinning behavior and σ the interfacial tension between gas and water. Shear-thinning behavior can be ignored by setting e c = 0.

Identification of SE and PB models at local steady state

As only local steady state foam flow is modeled by SE model, the identification procedure developed hereafter will involve only the solutions of PB models at local equilibrium. As explained in the previous section, the gas-phase mobility λ PB g in most PB models is based on a modified gas viscosity whereas the foam mobility λ SE g in SE models is based on a modified relative permeability of gas. Hence, the identification of SE-model gas mobility to PB-model gas mobility is

λ SE g = λ PB g , that is k SE rg /µ g = k rg /µ PB g .
Reintroducing the expressions of k SE rg and µ PB g given by Eqs ( 8), ( 9) and ( 5) we obtain:

(M ref -1) 4 ∏ i=1 F i = c f n f µ g |v PB g | 1/3 . ( 12
)
Assuming no oil and optimal concentration of foaming additive, the ∏ 4 i=1 F i product only involves the dry-out and shear thinning functions F 2 and F 4 given by Eqs ( 10) and [START_REF] Chambers | Capillary phenomena in foam flow through porous media[END_REF]. In order to go further with the identification, the interstitial velocity is expressed as a function of the capillary number N cg defined in Eq. [START_REF] Chambers | Capillary phenomena in foam flow through porous media[END_REF], such that Eq. ( 12) can be rewritten as:

(M ref -1) F 2 N ref cg N cg e c = c f n f µ g µ g N cg σ 1/3 . ( 13
)
Rearranging shows that the right member emanating from PB model includes the N cgdependent F 4 function corresponding to the shear-thinning regime:

(M ref -1) F 2 N ref cg N cg e c = c f n f µ g µ g N ref cg σ 1/3 N ref cg N cg 1/3 . (14)
At that stage, a first identification item can be pointed out, that concerns the shearthinning impact on the mobility reduction function. Previous equation shows that dependence on velocity, here expressed through N cg , involves an exponent e c equal to 1 3 according to the motion law of bubbles established by Bretherton [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF]. We note that both our SE model and under consideration PB model account for the shear-thinning effect on foam performance for any gas velocity, because the latter is assumed higher than a minimum value corresponding to N ref cg and satisfying foam generation. To go further on with F 2 identification, we note that, apart from the velocity-dependent N cg term, on the one hand the remaining variable is n f , and on the other hand F 2 is a function that is scaled between 0 and 1. Therefore, F 2 can be identified as follows

F 2 = An f , ( 15
)
where A is a constant. Allowing F 2 to reach the unity value requires to set A to the inverse of the maximum foam texture n max f . It is worth underlining that the shear-thinning mechanism is ignored in F 2 identification, because that mechanism is held by N cg -dependent function.

The dry-out function F 2 models the decreasing impact of foam on mobility when the porous medium induces a coalescence of foam bubbles because of capillary pressure effects on foam films stability [START_REF] Farajzadeh | Effect of permeability on implicit-texture foam model parameters and the limiting capillary pressure[END_REF][START_REF] Boeije | Fitting foam-simulation-model parameters to data: I. coinjection of gas and liquid[END_REF][START_REF] Ma | Nonuniqueness, numerical artifacts, and parameter sensitivity in simulating steady-state and transient foam flow through porous media[END_REF]. In other words, the volumetric concentration of bubbles, i.e. the foam texture n f , decreases rapidly as gas saturation exceeds a certain threshold value that limits films stability. However, that coalescence of the foam occurs over a certain range of capillary pressure values (or of films disjoining pressures). Such a behavior is well represented by the F 2 function given by Eq. ( 10), F 2 (S w ) = + δ ] over which foam collapses and the Θ parameter drives the width of that saturation range.

The F 2 function accounts directly for the dependence of foam performance on the porous medium saturation whereas in PB models at local steady state, the same dependence is modeled indirectly via the foam texture that is itself modeled as a function of saturation. Considering the Kam's model taken as PB model type, the threshold saturation S * -PB w involved in that PB model is analog but not exactly equal to the water saturation threshold S * -SE w of the SE model, because S * -PB w is the lowest water saturation value below which foam no more exists, which corresponds to a value of the F 2 function equal to 0. For the sake of clarity, we will then keep on distinguishing S * -SE w and S * -PB w , whose difference is given by δ

= S * -SE w -S * -PB w = F -1 2 ( 1 2 ) -F -1 2 (0) where F -1
2 is the inverse of the F 2 (S w ) function. Obviously, regarding the F 2 function, the identification of SE model to PB model on the basis of their respective gas mobilities at local steady state as written before in Eq. ( 14) does not help to determine directly Θ and S * -SE w as an expression of PB model parameters. Thus, F 2 can only be identified globally as the dimensionless foam texture n f /n max f , knowing however that all the underlying physics of porous medium effects as described before is held in the foam texture, that is the local steady state solution of the PB model. That is, we can write

F 2 = n f n max f = n D , (16)
where n D is the dimensionless foam texture. Different values of the foam texture n i f and of the corresponding saturation S i w can be determined by solving PB model at local steady state for a set of foam quality values f i g , as explained in Section 4. That is, Θ and S * -SE w constants can be determined through an optimization of a set equations written as w )] = ε, that can also be written as 1 2 -1 π arctan (Θ δ ) = ε, which yields Θ = 1 δ tan π 1 2ε . The choice of tolerance ε is discussed later on in Section 5.3. To end up with the identification of SE model empirical parameters to PB model, we rewrite Eq. ( 14) as

n i f (S i w ) n max f = 1 2 + 1 π arctan Θ S i w -S * -SE w . (17)
(M ref -1) F 2 N ref cg N cg e c = c f n max f µ g µ g N ref cg σ 1/3 n D N ref cg N cg 1/3 . ( 18
)
Taking into account previous identifications, M ref expression versus PB model parameters is

M ref = 1 + c f n max f µ g µ g N ref cg σ 1/3 . ( 19
)
Identification summary is then:

M ref = 1 + c f n max f µ g µ g N ref cg σ 1/3 , F 2 = n f n max f = n D , F 4 = N ref cg N cg 1/3 . ( 20
)
To proceed with the identification of an SE model using the results of the above identification, we need to know the foam texture n f for various foam flow conditions. This is provided by the solution of PB models at local steady state which is detailed hereafter.

4 Analytical approach to relate gas saturation and foam texture at local steady state

In this section, we propose an analytical approach to calculate foam texture at local steady state. The texture is obtained by equating the generation and coalescence rates of the PB model given by Eq. ( 6). The basic idea consists in solving jointly the foam-gas fractional flow equation and the local equilibrium equation, with gas saturation and foam texture as unknowns. Indeed, the foam-gas fractional flow f g is a known parameter of the problem at local steady state. The problem set-up is worth being detailed due to the peculiarity of entailed resolution.

Herein, a two-phase gas-water fluid system is considered. Capillary pressure and gravity are neglected. Porosity is assumed constant. A one-dimensional incompressible flow is considered without any sink terms. The PB model equations ( 1) and ( 4) can be combined into the following fractional flow equation

∂ t (S g ) + u Φ ∂ x f g (S g , n f ) = 0, (21)
where, for the PB model, the gas fractional flow rate depends on gas saturation and foam texture:

f g (S g , n f ) = 1 1 + µ g + c f n f (v PB g ) 1/3 k rw µ w k rg . (22)
That equation can be inverted in order to express the foam texture as a function of saturation and gas fractional flow, that is, denoting n FF f this fractional-flow-derived foam texture:

n FF f = u ΦS g 1/3 1 -f g 1 + µ g µ w k rw k rg f 2/3 g c f µ w k rw k rg . ( 23
)
At local equilibrium, the generation and coalescence rates of the PB model are equal. According to Eq. ( 6) and considering the maximum texture limitation, this local equilibrium yields the following expression of the foam texture at local equilibrium denoted n LE f :

n LE f = c g c c S w S w -S * -PB w n |∇P| m if n f < n max f , n max f otherwise, ( 24
)
where ∇P is the pressure gradient which can be expressed from Darcy equation applied to the water phase, whose flow properties are not affected by the presence of foam:

|∇P| = µ w u(1 -f g ) kk rw . ( 25
)
As explained in Section 2.2, the maximum foam texture n max f is estimated as the characteristic pore size of the porous medium under consideration. This is an approximation for natural porous media characterized by a more or less complex pore size distribution. Therefore, it remains a variable subject to adjustment in PB models at local equilibrium [START_REF] Chen | Modeling foam displacement with the local-equilibrium approximation: theory and experimental verification[END_REF].

Following the expression of foam texture at local equilibrium given by Eq. ( 24), problem solving should distinguish two cases:

-For any f g within the low-quality regime, foam texture is constant and is equal to the maximum value (n f = n max f ). Thus, only gas saturation remains to be calculated from the gas fractional flow equation [START_REF] Lawson | Alternate slugs of gas and dilute surfactant for mobility control during chemical flooding[END_REF], which yields

S g = u Φ   1 -f g 1 + µ g µ w k rw k rg n max f f 2/3 g c f µ w k rw k rg   3 . (26)
-For any given f g within the high-quality regime (n f < n max f ), equating Eqs ( 24) and ( 23) yields the following gas saturation equation:

n LE f (S g ) c g c c (1 -S g ) 1 -S g -S * -PB w n µ w u(1 -f g ) kk rw m = n FF f (S g ) u ΦS g 1/3 1 -f g 1 + µ g µ w k rw k rg f 2/3 g c f µ w k rw k rg , (27)
where n LE f represents the foam texture obtained from the local equilibrium of the PB model and n FF f is the foam texture computed from the gas fractional flow equation in the presence of foam.

Equations [START_REF] Lake | Enhanced oil recovery[END_REF] and [START_REF] Green | of SPE Textbook Series[END_REF] have to be solved with numerical methods, as described later on in Section 5.

Finally, a series of measurements varying the foam quality at fixed total injection rate is needed in order to relate the gas saturation to the foam texture at local equilibrium for any quality of the foam under consideration. This relationship underlies the possibility to identify an SE model to a PB model.

Application to a foam data set: identification of our SE model to calibrated Kam et al PB model

We apply the previous identification method to calibrate our SE model to a given experimental data set published by Alvarez et al [START_REF] Alvarez | Unified model for steady-state foam behavior at high and low foam qualities[END_REF]. Specifically, the authors measured the steady state apparent foam viscosity for several foam qualities at constant total injection rate. All these experiments are two-phase flow displacements involving water, gas and surfactant.

The considered PB model for the application of the calibration method is previouslydescribed in the Section 2.2. We focus on the steady state foam flow that is representative of flow conditions achieved in a reservoir beyond the near-wellbore region where foam generation conditions are fulfilled. [START_REF] Alvarez | Unified model for steady-state foam behavior at high and low foam qualities[END_REF] The experiments of Alvarez et al [START_REF] Alvarez | Unified model for steady-state foam behavior at high and low foam qualities[END_REF] are conducted on a Berea sandstone core with a N 2based foam at ambient temperature and an outlet pressure of 41 bar. The core is 5 cm in diameter and 28 cm in length. The aqueous phase of foam is a 1 wt % solution of surfactant in a brine containing 3 wt % NaCl and 0.01 wt % CaCl 2 .

Case study: Alvarez et al experiments

A pre-generated foam of given quality is injected in the water-saturated core. While keeping constant the total flow rate of gas and water, the foam quality is changed from one experiment to another to cover a range from 0.25 to 0.91. The measured apparent foam viscosity ranges from 274 to 905 cP. The transition between the two regimes occurs at the optimal foam quality f * g = 0.78 for which the maximum apparent foam viscosity of 905 cP is reached, as reported in Fig. 4. Table 1 summarizes the rock and fluids properties of these foam displacement experiments.

Recently, Lotfollahi et al [START_REF] Lotfollahi | Comparison of implicit-texture and population-balance foam models[END_REF] performed an extensive treatment of the foam apparent viscosity data recorded during these experiments. They calibrated different SE models, as well as three PB models, among which Kam et al model [START_REF] Kam | Dynamic simulations with an improved model for foam generation[END_REF]. Specifically, adjustable parameters were c g /c c , c f , n, m and S * -PB w . The adjusted parameter values that fit the experimental data are reported in Table 2 (gas and water relative permeabilities in the absence of foam were adjusted with power laws, namely k rw = k 0 rw S n w and k rg = k 0 rg (1 -S) n g with S = S w -S wi 1-S gr -S wi ). In the following section, we first compute the local steady state foam texture using the Kam et al model parameters given in Table 2, since it is the first step for SE model identification to PB model according to Eq. (20).

Determination of foam texture at local equilibrium using Kam et al PB model: a posteriori validation

As explained in Section 4, one can solve Kam et al model at local equilibrium by equating Eqs ( 24) and [START_REF] Lawson | Alternate slugs of gas and dilute surfactant for mobility control during chemical flooding[END_REF]. Accordingly, two regimes should be distinguished for foam texture calculations. On the one hand, within the low-quality regime, the foam texture n f is known and is equal to the maximum value n max f so that only gas saturation remains to be calculated using Eq. ( 26). On the other hand, within the high-quality regime, one should first determine the gas saturation by solving Eq. ( 27), and after use Eq. ( 23) to calculate the corresponding foam texture value.

In practice, for any foam quality value f g ∈ ]0, 1[, the solution S g is obtained using the following algorithm:

1. Find the gas saturation by solving the equality [START_REF] Green | of SPE Textbook Series[END_REF]. 2. Calculate the corresponding foam texture from Eq. ( 23). 3. If the obtained foam texture is lower than the limiting foam texture n max f , a high-quality regime is obtained, and in this case the resolution process is iterated with another f g . Otherwise, a low-quality regime is achieved with a maximum foam texture that is equal to n max f , and the gas saturation has to be determined again from Eq. ( 26) before proceeding with another f g . This algorithm computes the exact fractional flow of gas in the presence of foam as a function of the gas saturation, hence the foam texture at local equilibrium. However, the uniqueness of the solution S g within the two regimes has to be checked. For this reason, we plot in Fig. 2 the foam textures n LE f (S g ) and n FF f (S g ) for several foam qualities and we verify the solution uniqueness. Their intersection yields the gas saturation that corresponds to a given foam texture within the high-quality regime. Intersection of n max f line and n FF f yields the gas saturation within the low-quality regime.

Fig. 2 shows that the intersection points between n max f line and n FF f on the one hand, and n LE f (S g ) and n FF f (S g ) on the other hand, are unique for each foam quality (see also Table 3 which reports a few saturation and texture solutions). Thus, the optimal set of Kam et al model parameters ensure the uniqueness of the solution (foam texture and gas saturation) for both low-and high-quality regimes.

Finally, to check the validity and accuracy of that PB model, the apparent foam viscosity data from Alvarez et al experiments were recalculated using above foam texture and gas saturation solutions. More precisely, for any given value of foam quality, the foam apparent viscosity is recalculated as 1

µ app f = k rw µ w + k rg µ PB g (derived from u = u w + u f g = k|∇P| µ app f = kk rw |∇P| µ w + kk rg |∇P| µ PB g ) where µ PB g = µ g + Φ u 1/3 c f n f ( f g /S g) 1/
3 is the effective gas viscosity given by Eq. ( 5), and where n f and S g are the foam texture and the saturation solutions determined at local equilibrium as before. Table 3 reports the detailed results of this calculation for a few qualities ranging from 0.1 to 0.9.

Fig. 4 reports the obtained results over the entire foam quality range: one can observe the very good agreement of the computed apparent foam viscosity from Kam et al model with the measurements of Alvarez et al [START_REF] Alvarez | Unified model for steady-state foam behavior at high and low foam qualities[END_REF], both at the low-and high-quality regimes. Furthermore, the obtained results point out an abrupt foam collapse that occurs at almost a single water saturation S * -PB w = 0.31. Thus, a PB model, such as Kam et al model, can be used to properly fit steady state foam mobility at both high-quality and low-quality regimes.

= 19.1 × 10 12 m -3 ). f g (-) S g (-) n f (m -3 ) n D (-) µ PB g (cP) µ app f (cP) 0.

Identification of SE model to calibrated PB model

The Kam et al PB model, as calibrated and validated before, is used to directly identify the SE model following the procedure of Section 3, with the relationship between fluid saturation and foam texture as determined in Section 4. First, we determine the limiting water saturation S * -SE w by inverting the computed dimensionless foam texture function n D . We found that S * -SE w = n -1 D ( 1 2 ) = 0.313. In addition, the parameter Θ is calculated ensuring F 2 S * -PB w = ε. A tolerance ε of 0.05 is applied in order to render at best the dimensionless foam texture variation, leading to a value of Θ = 2.1 × 10 3 , as reported in Table 4. Fig. 3(a) shows the resulting dry-out function F 2 and the dimensionless foam texture function n D as functions of water saturation.

The reference capillary number, given by Eq. [START_REF] Chambers | Capillary phenomena in foam flow through porous media[END_REF], is set at the lowest reached value during the experiments, namely N ref cg = (µ g f min g u)/(σ ΦS min g ) where f min g = 0.25 is the lowest foam quality. The lowest gas saturation S min g is calculated as follows. By application of the generalized Darcy law to the water phase one has u w = kk rw µ w |∇P| hence k rw = µ w u f w k|∇P| . With a k rw power law function such that k rw = k 0 rw S n w with S = S w -S wi 1-S gr -S wi , S w can be written as 4. The obtained shearthinning function F 4 is reported in Fig. 3 The SE model parameters, which are reported in Table 4, are then used to calculate the apparent foam viscosity µ app f which can be written in the SE model framework as 1

µ app f = k rw µ w + FM•k rg
µ , derived again from the expression of total velocity as the sum of water and gas velocities, with gas velocity from Eq. ( 8) where FM is given by Eq. ( 9), and only involves the product of F 2 and F 4 functions. The water saturation corresponding to each foam quality is obtained by solving the gas fractional flow equation 5 reports the detailed results of this calculation for a few qualities ranging from 0.1 to 0.9. The so obtained apparent foam viscosity, which is reported in Fig. 4(a), is very close to the measured one and to the one derived from the Kam PB model over the entire foam quality interval.

f g = 1/(1 + µ g µ w k rw FM•k rg ). Table
In practice, the SE parameters would be used to simulate the displacement for several foam qualities. The simulated steady state pressure gradient |∇P| = ∆ P L allows to compute = k|∇P| u by application of the one-phase Darcy law to the equivalent fluid. Table 6 reports the so obtained apparent foam viscosity derived from the simulated pressure gradient for several qualities, with our SE model. As expected, these simulated values are very close to the ones reported in Table 5, since most of the relative errors with respect to the reported values in Table 5 are of the order of 1%. The apparent foam viscosities computed from SE model and PB model at local equilibrium are almost equal over the whole range of foam qualities. The maximum apparent foam viscosity and the optimal foam quality are also successfully predicted by the SE model. We note a slightly smoother transition between low-and high-quality regimes than the sharp transition predicted by Kam et al model. A larger value of Θ parameter could have been used to refine the fit in the transition region around the optimal foam quality f * g . Table 7 reports the so obtained PB and SE models apparent foam viscosities and their relative errors with respect to the measurements

∆ PB = |(µ app f ) PB -µ app f |/µ app f and ∆ SE = |(µ app f ) SE -µ app f |/µ app f
. These relative errors are also reported in Fig. 4(b). Overall, they do not exceed 5% except for large qualities. Thus, both models adjust equally well the steady state measurements of the apparent foam viscosity.

To conclude, the main input for SE model identification to PB model is the relationship between fluid saturation and the foam texture predicted by PB model. As expected, the resulting SE model is found to fit foam quality-scan experiments as well as the original PB model. This opens the way for reliable use of foam flow models for reservoir applications.

Comparison with other calibration methods

Several adjustment procedures of the apparent foam viscosity measurements to a single scan of foam quality are available in the literature. They are reviewed and discussed hereafter. Boejie and Rossen [START_REF] Boeije | Fitting foam-simulation-model parameters to data: I. coinjection of gas and liquid[END_REF] proposed a simple manual procedure to fit SE models to the apparent foam viscosity measurements. The procedure is based on six steps dealing separately with the data in the high-quality and low-quality regimes. The latter approach assumes an abrupt transition between the two regimes such that large values of Θ are considered beforehand (Θ between 10 4 and 10 5 ). In addition, the pressure gradient is assumed to decrease linearly with the foam quality, in the high-quality regime. The other SE model parameters are calculated step by step. The proposed approach is simple as it can be carried out directly from the plot of the apparent foam viscosity versus foam quality. Nevertheless, the method cannot fit properly the measurements when the experimental data show a gradual transition between high-and low-quality regimes. In that case, the authors propose to use the results of their method to initiate a least squares fit of all foam model parameters to further tune the SE model.

Lotfollahi et al [START_REF] Lotfollahi | Comparison of implicit-texture and population-balance foam models[END_REF] used a least squares minimization method to fit different SE models and PB models to two experimental data sets from Alvarez et al [START_REF] Alvarez | Unified model for steady-state foam behavior at high and low foam qualities[END_REF] and Moradi-Araghi et al [START_REF] Moradi-Araghi | Laboratory evaluation of surfactants for CO 2 -foam applications at the South Cowden Unit[END_REF]. The fit is performed within the two regimes using all the measurements. The SE model parameters to be optimized are M ref , S * -SE w , Θ and e c . The authors indicate that the final fit may depend on the initial guess, that is the method does not guarantee to yield the optimal model. This non-uniqueness is also pointed out by Kapetas et al [START_REF] Kapetas | Effect of permeability on foam-model parameters -an integrated approach from coreflood experiments through to foam diversion calculations[END_REF], who obtained different sets of foam parameters when applying the two previous methods on a large set of experimental data.

To improve the model match, Farajzadeh et al [START_REF] Farajzadeh | Effect of permeability on implicit-texture foam model parameters and the limiting capillary pressure[END_REF] used data weighting and constraints in the least squares minimization method. In addition, they added liquid saturation values for each measurement to the definition of the objective function which reads in this case O(X) = Assuming no shear-thinning effect in foam quality-scan experiments, Ma et al [START_REF] Ma | Estimation of parameters for the simulation of foam flow through porous media. part 1: the dry-out effect[END_REF] proposed a graphical method with M ref and S * -SE w as variables, in order to match the transition foam viscosity between the low-and high-quality regimes. The Θ parameter is matched after with a least squares method applied to the other steady state apparent viscosity measurements. Later on, Ma et al [START_REF] Ma | Nonuniqueness, numerical artifacts, and parameter sensitivity in simulating steady-state and transient foam flow through porous media[END_REF] determined the same three model parameters (M ref , S * -SE w , Θ ) by using a least squares minimization method as other authors, but adding a constraint on S * -SE w to avoid non unique (M ref , S * -SE w ) solutions. Then, Zeng et al [START_REF] Zeng | Insights on foam transport from a textureimplicit local-equilibrium model with an improved parameter estimation algorithm[END_REF] converted the three-parameter (M ref , S * -SE w , Θ ) optimization problem (Ma et al [START_REF] Ma | Nonuniqueness, numerical artifacts, and parameter sensitivity in simulating steady-state and transient foam flow through porous media[END_REF]) to a single-parameter (M ref ) problem by using steady state apparent viscosity data and searching the M ref value that fits a linear function of the saturation with (S * -SE w , Θ ) as parameters. The same authors also adjusted the shear-thinning function separately from total-flow-rate-scan experiments.

To end with, this limited survey of SE models indicate that calibration methods are based on a direct identification of the parameters of the SE model relative permeability function, without using the physical representation of foam as a succession of lamellas. Many alternative methods of implementation have been developed to deal with the multi-variable nature of the identification problem and with the peculiarity of foam flow equations involving discontinuities in solutions. The method of this paper aims at overcoming these difficulties via the PB model representation as flow-restricting lamellas.

Conclusions and perspectives

The objective of this paper was to establish the equivalence between semi-empirical foam flow models, based on a modified relative permeability concept, and population-balance models that keeps track of the intrinsic physics of foam that explains mobility reduction in porous media. The identification of both foam flow models leads to a method for calibrating SE models from the PB model translation (i.e. foam-texture translation) of foam flow data. This physical approach of foam flow could circumvent some difficulties in the direct calibration of SE model from foam mobility (or apparent viscosity) data. The use of least squares method is not really adapted to the calculation of a flow model involving discontinuities as observed with foams. Lamellas-based models involve the foam texture physical variable whose evolution is precisely responsible of more or less pronounced foam flow discontinuities. The method developed and discussed in this paper calls for further developments regarding the impact of capillary pressure on the foam SE model parameters. In addition, the question has to be addressed for reservoir applications involving different rock permeabilities and pore structures.

Fig. 1 :

 1 Fig. 1: Limiting capillary pressure schematics [10, 12]. (a) Capillary pressure as a function of water saturation: capillary pressure increases with the gas saturation up to P * c , above which foam lamellas become unstable and coalesce; (b) fractional flow of gas as a function of water saturation and bubble size: when water saturation reaches S * w corresponding to P * c , foam becomes coarser, gas fractional flow increases and consequently gas mobility increases. At low-quality regime, for S w > S *w , the fractional flow of gas follows the green curve, then a transition occurs at S w = S * w (dotted line), and finally the fractional flow of gas follows the red curve at high-quality regime for S w < S * w .

1 2 + 1 π

 1 arctan Θ S w -S * -SE w , where the S * -SE w saturation corresponds to the middle of the saturation range [S * -SE w δ , S * -SE w

+ 1 π

 1 Alternatively, one can determine directly the limiting water saturation S * -SE w in SE models from the dimensionless foam texture function n D as S * -SE w = n -1 D ( 1 2 ), where n -1 D is the inverse of the n D (S w ) function. Once S * -SE w is calculated, the parameter Θ can be determined by setting the F 2 function to zero at S * -PB w , i.e. F 2 S * -PB w = 0. As such a value cannot be strictly attained, a tolerance ε should be permitted on F 2 (S * -PB w ) such that 1 2 arctan [Θ (S * -PB w -S * -SE

Fig. 2 :

 2 Fig.2: Derived foam texture at local equilibrium using Kam et al optimal model adjusted by Lotfollahi et al[START_REF] Lotfollahi | Comparison of implicit-texture and population-balance foam models[END_REF] on Alvarez et al experiments[START_REF] Alvarez | Unified model for steady-state foam behavior at high and low foam qualities[END_REF]. The points corresponding to the solutions for each foam quality are given by circular markers (•) for n f = n max f and square markers (2) for n f < n max f (n max f = 19.1 × 10 12 m -3 ).

S w = S wi + ( 1 -.= 1 -.

 11 S wi -S gr ) Replacing k rw by the above expression versus |∇P| we obtain S w = S wi + (1 -S gr -S wi ) µ w u f w k|∇P|k 0 rw 1/n w . Applying the one-phase Darcy law to the equivalent fluid yields µ app f = k|∇P| u , and as a consequence S w = S wi + (1 -S wi -S gr ) µ w (1f g ) lowest gas saturation S min g can be calculated as S min g S wi -(1 -S gr -S wi ) µ w (1f min g ) Specifically, one has f min g = 0.25 and S min g = 0.63 for Alvarez et al measurements, hence N ref cg = 1.27 × 10 -8 , as reported in Table

Fig. 3 :

 3 Fig. 3: (a) Dimensionless foam texture obtained explicitly from Kam et al model and implicitly from our model to fit Alvarez et al measurements; (b) Shear-thinning function F 4 used in our SE model.

Fig. 4 :

 4 Fig. 4: (a) Comparison of Alvarez et al apparent foam viscosity measurements with the apparent viscosity computed from SE model with parameters directly estimated from the calibrated Kam et al model; (b) Apparent foam viscosity relative errors with respect to the measurements.

f

  , j (X, S w, j )µ exp f , j ] 2 , where X = (M ref , S * -SE w ,Θ , e c ) is vector of parameters to be optimized, w j weights, µ app f , j computed apparent foam viscosity, µ exp f , j measured apparent foam viscosity and S w, j is liquid saturation for each data point which is calculated from the liquid relative permeability as S w, j = S wi + (1 -S gr -S wi )µ w (1f g, j ) µ exp f , j k 0 rw 1/n wfor power law k r functions, as explained in Section 5.3 [see Eq. (28)].

Table 1 :

 1 Rock and fluids properties of Alvarez et al foam displacements[START_REF] Alvarez | Unified model for steady-state foam behavior at high and low foam qualities[END_REF].

	Gas T ( • C)	N 2 25	u (cm/hour) f g range (-)	3.175 0.25 -0.91
	P (bar)	41	µ w (T, P) (cP)	0.7
	σ (mN/m) 30	µ g (T, P) (cP)	0.02
	k (mD)	530	ρ w (T, P) (kg/m 3 ) 998
	Φ (-)	0.18	ρ g (T, P) (kg/m 3 )	46

Table 2 :

 2 Kam et al model (2007) parameters adjusted on Alvarez et al experiments[START_REF] Alvarez | Unified model for steady-state foam behavior at high and low foam qualities[END_REF] by Lotfollahi et al[START_REF] Lotfollahi | Comparison of implicit-texture and population-balance foam models[END_REF].

	Gas and water relative permeabilities		Kam et al model parameters
	S wi (-) S gr (-)	0.2 0.2	S * -PB w n max f	(-) (m -3 )	0.31 19.1 × 10 12
	k 0 rw (-) k 0 rg (-)	0.2 0.94 c g /c c (-) c f (m 10/3 hour -1/3 cP)	2.77 × 10 -11 8 × 10 16
	n w (-)	4.2	n (-)		1.91
	n g (-)	1.3	m (-)		0.2

Table 3 :

 3 Steady state foam properties predicted by Kam et al PB model (n max f

Table 4 :

 4 SE model parameters directly estimated from the calibrated PB model.

	M ref (-)	S * -SE w	(-) Θ (-)	N ref cg (-)	e c (-)
	6.4637 × 10 4	0.313	2.1 × 10 3	1.27 × 10 -8	1 3

Table 5 :

 5 Steady state foam properties predicted by our SE model.

	f g (-)	S w (-)	F 2 (-)	F 4 (-)	FM (-)	µ	app f	(cP)
	0.1	0.4013 0.9983 1.3254 0.1169 × 10 -4	309
	0.3	0.3631 0.9970 0.9381 0.1654 × 10 -4	582
	0.5	0.3407 0.9945 0.8004 0.1944 × 10 -4	774
	0.7	0.3199 0.9780 0.7229 0.2188 × 10 -4	909
	0.9	0.3128 0.3620 0.6671 0.6407 × 10 -4	391
	the foam apparent viscosity µ	app f

Table 6 :

 6 Apparent foam viscosity, for several qualities, predicted by our SE model from the simulated steady state pressure gradient |∇P| = ∆ P L across the core. Relative errors ∆ S w and ∆ µ app f refer to steady state values reported in Table 5.

	app f	(cP) ∆ µ app f	(%)

f g (-) S w (-) ∆ Sw (%) ∆ P (bar) |∇P| (bar/m) µ

Table 7 :

 7 Apparent foam viscosities and relative errors with respect to the measurements, both for Kam PB model (∆ PB ) and our SE model (∆ SE ; see the text).

		Measurements f g (-) µ app f (cP)	Kam et al PB model µ app f (cP) ∆ PB (%)	µ	Our SE model app f (cP) ∆ SE (%)
		0.25	542.2	520.8	4	519.5	4
		0.30	607.4	594.3	2	592.6	2
	Low-quality regime	0.36 0.45 0.59	674.0 740.7 810.6	653.0 730.7 850.4	3 1 5	650.9 727.6 843.9	3 2 4
		0.67	876.5	910.1	4	896.9	2
		0.70	905.6	931.1	3	910.6	1
		0.78	881.4	839.7	5	836.2	5
		0.80	808.6	745.3	8	749.0	7
		0.82	743.7	685.9	8	689.6	7
		0.85	680.0	587.7	14	589.9	13
	High-quality	0.86	609.8	541.2	11	542.0	11
	regime	0.87	541.7	506.6	6	506.9	6
		0.88	476.1	470.6	1	470.0	1
		0.89	410.2	434.1	6	433.0	6
		0.90	339.7	387.9	14	385.8	14
		0.91	273.9	360.3	32	357.9	31


David Pichardie
email: pichardy@eleves.bretagne.ens-cachan.fr

Yves Bertot
email: yves.bertot@inria.fr

Formalizing Convex Hulls Algorithms

We study the development of formally proved algorithms for computational geometry. The result of this work is a formal description of the basic principles that make c o n vex hull algorithms work and two programs that implement c o n vex hull computation and have been automatically obtained from formally veri ed mathematical proofs. A special attention has been given to handling degenerated cases that are often overlooked by conventional algorithm presentations.

Algorithms to compute the convex hull of a collection of points in two or three dimensions abound. The complexity of this problem is known to be approximately the same as for a sorting algorithm in the worst case, but the average complexity depends very much on the distribution of data: some algorithms actually have a complexity of the form n h, where n is the total numb e r o f p o i n ts and h is the number of points on the convex hull. These algorithms will be more e cient than nlog(n) when few points are on the convex hull. This makes it useful to have s e v eral algorithms around.

We h a ve studied two algorithms. The rst is an incremental algorithm, where new points are progressively added to the input data. If the new point is outside the convex hull, then it is necessary to add two new edges to the convex hull and remove some others. The second algorithm is known as the package-wrapping algorithm: it follows the intuition of tying a band around the collection of points.

At the n th step, the band is already in contact with a few points and the algorithm proceeds by computing what the next point is going to be: if with the band turning clockwise around the points, the next point i n t h e h ull is going to be the leftmost remaining point.

All these algorithms rely on an orientation predicate that expresses whether the points of a triangle are enumerated clockwise or counter-clockwise. This orientation predicate is easily computed using the coordinates of the points, but it is meaningless when the points are aligned. Usual presentation of convex hulls algorithms assume that three points in the input data are never aligned.

The structure of convex hull algorithms rely on the fact that the orientation predicate satis es some properties: for instance the triangles built with four points have t o b e o r i e n ted in a consistent manner. Knuth 8] describes the minimal properties of the orientation predicates and calls them axioms.

Knuth's approach of axioms for convex hulls has the nice property o f s e p arating concerns about the properties of arithmetic expressions containing point coordinates and the control structure of algorithms. It has proved a good organization principle for our description and proof development.

Thus, our work contains distinct parts. The rst part is about the axioms, and showing that they hold for an implementation predicate. This part involves some numerical computation. The second part describes the main structure of the algorithms based on the axioms, with numeric computation completely avoided. In the last part, we revisit the algorithms to make them robust with respect to degenerated data.

All our proofs have been mechanically checked using the Coq system 2].

Related work

Automatic theorem proving and geometry have been in contact for some time.

The work in this domain that we are better aware of is that of Chou 5]. In this development, theorems about basic geometric gures like straight lines, triangles and circles are proved automatically, but there is no result about computational geometry per se, since algorithms are not the object of the study. Puitg and Dufourd 9], used the Coq proof system to describe notions of planarity. W ork on nding minimal axioms to describe geometric concepts has also been done in constructive proof theory by Jan von Plato 11] and formalized in Coq by Gilles Kahn 7].

Last, we should refer to all the work done on computational geometry, b u t this domain is far too large to be cited entirely, w e can only refer to the books we h a ve used as reference, 6] and 3], and 10].

Knuth's axioms for convex hulls

In what follows, we assume that points are taken from a set P and we describe the orientation predicate as a predicate over three points p, q, r in P, which w e denote d pqr.

Knuth's axioms describe the various ways in which triangles can be looked at and arranged on the plane. The rst axiom expresses that when enumerating the points of a triangle, one may start with any of them: Axiom 1 8p q r: d pqr ) d qrp:

As a corollary of this axiom one also has d pqr ) d rpq. This axiom is cyclic it can be repeated inde nitely on the same data. For this reason, it impedes automatic proof search procedures.

The second axiom expresses that if a triangle is oriented counter-clockwise, then the same triangle with two p o i n ts transposed is in fact oriented clockwise.

Axiom 2 8p q r: d pqr ) : d prq:

An immediate consequence of the rst two axioms is that three points are oriented only if they are pairwise distinct: 8p q r: d pqr ) p 6 = q ^q 6 = r ^r 6 = p Many of the lemmas appearing in our development therefore rely on the assumptions that the points being considered are pairwise distinct. We a l s o h a ve a predicate on lists of points expressing that they contain pairwise distinct points.

The third axiom expresses that a triangle is either oriented clockwise or counter-clockwise.

Axiom 3 8p q r: p 6 = q ) q 6 = r ) p 6 = r ) d pqr _ d prq:

This a rmation voluntarily overlooks the fact that points may be aligned. This is a tradition of computational geometry that algorithms can be studied without taking care of what are called degenerate cases. We w i l l h a ve to study this more carefully later.

The fourth axiom expresses that the four triangles obtained with four arbitrary points may not be oriented in an arbitrary manner: there has to be some sort of consistency, easily understood by observing gure 1. Axiom 4 8p q r t : c tqr ^c ptr ^c pqt ) d pqr:

The fth axiom expresses that the orientation predicate may be used to sort points in some way, with a notion that is similar to transitivity. Still one has to be careful to avoid going in circles, this is the reason why this axiom uses 5 points, see gure 2. Clever algorithms use this property o f l i k e-transitivity t o avoid looking at all points.

Axiom 5 8p q r s t : c tsp ^c tsq ^c tsr ^c tpq ^c tqr ) c tpr t s p q r Fig. 2. Axiom 5: sorting points in a half plane.

The rst three orientation properties express that the three points p, q, and r are in the same half plan. The other three orientation predicates express transitivity: from the point of view of t, i f q is left of p and r is left of q, t h e n r is left of p.

There is a variant t o this axiom, that can actually be proved using axioms 1, 2, 3, and 5, and which w e w i l l r e f e r t o b y the name axiom 5 0 . Axiom 5 0 8p q r s t : c stp ^c stq ^c str ^c tpq ^c tqr ) c tpr Knuth shows that these axioms are independent. For the case of general positions, these axioms are enough, but if we w ant our algorithms to be robust for possibly aligned points in the input, there are two possibilities. The rst one is to introduce notions that are meaningful for triplets of points that are aligned. In this case, we also use segments and the property that a point belongs to a segment given by its two extremities. This adds up to 9 more axioms, which we do not describe here for the moment. The second solution is to nd a w ay to extend orientation to aligned points in a way that still satis es the axioms. Perturbation methods as described in 1] provide a satisfactory approach t o t h i s .

Proving the axioms

When points are given as elements of R 2 , the orientation predicate is given by observing the sign of a determinant: d pqr x p y p 1 x q y q 1 x r y r 1 > 0:

In the following we will write jpqrj for this determinant. Knuth's axioms are not proper axioms anymore, they become simple consequences of the properties of addition and multiplication of real numbers. The proof of Axioms 1 and 2 is given by t wo equations: jpqrj = jqrpj jpqrj = jprqj:

From a mechanical proof point of view, the proofs of these equalities is done painlessly by the usual con uent rewriting system for ring structures, as implemented by the Ring tactic in Coq. Axiom 3 simply does not hold in general.

When points are aligned, the determinant is null. A common solution to this problem is to consider only data that are in general position, that is, assume that no three points given in the problem are ever aligned. We will rst follow this practice to get a better idea of the algorithms. In a later section, we w i l l s e e di erent solutions to the problem of considering degenerated cases.

In man-made proofs, Axiom 4 relies on an equality obtained from a 4-by-4 determinant that is obviously null, since it has 2 identical columns.

x p y p 1 1 x q y q 1 1 x r y r 1 1 x t y t 1 1 = jpqrj ; j tqrj ; j ptrj ; j pqtj = 0 :

Proving the right-hand side equality is also a simple matter of equality i n a r i n g structure. Now, if jtqrj, jptrj, jpqtj are all strictly positive, then jpqrj is.

Verifying Axiom 5 also relies on an equation, known as Cramer's equation, which has the following form: jstqjjtprj = jtqrjjstpj + jtpqjjstrj Here again, automatically proving this equation is a simple matter of ring based computations.

The conclusion of this section on the orientation predicate is twofold. First numeric computation can be done in a ring rather than in a complete eld. This is important as one of the reasons for problems in geometric algorithms is the loss of precision in oating point arithmetic. For convex hulls, oating points are not required. The points may b e g i v en by i n teger coordinates in a small grid. If oating points are used, one may notice that, since the determinant b a s i c a l l y i s a polynomial of degree two in the coordinates of the three points being considered, it is only necessary to compute the orientation predicate with a precision that is double the precision of the input data.

The speci cation

What is the convex hu l l o f a s e t o f p o i n ts S? One de nition is to describe it as a minimal subset S 0 S, s u c h that all points in S can be described as positive linear combinations of points in S 0 .

An alternative description of the convex hull uses oriented edges. A possibly open line S 0 encompasses a set of points S if for any p o i n t p in S and any t wo consecutive p o i n ts t i t i+1 in S 0 , \ t i t i+1 p or p 2 f t i t i+1 g holds. We consider that a list l is a c o n vex hull of S if l is included in S, l is non-empty if S is not, l contains no duplicates, and if t is the last element o f l, t h e n t l, formed by adding t at the head of l, encompasses S. In other terms, the convex hull is de ned as a minimal intersection of half-planes.

When a data set is the empty set, a singleton, or a pair, then the convex hull is simply the list enumerating the set's elements.

Knuth chooses the second de nition in his monograph. We will also, for several reasons. First, lists are a natural data-structure in the functional programming language we u s e i n our proofs. Second, this de nition also relies on the orientation predicate, thus making it possible to forget everything about numerical computation in the main structure of algorithm. Third, most algorithms to compute convex hulls naturally rely on the fact that the intermediary data they construct already is a list of points included in S.

The three parts of the speci cation are given using inductive de nitions on lists (actually we also consider that the input data is given as a list of points).

Still, we h a ve made an e ort to establish a bridge between the two de nitions.

We h a ve mechanically veri ed a proof that if S 0 is a convex hull of S and p is an arbitrary point i n S, then there exists three points t i , t j , a n d t k 2 S 0 such t h a t p is inside the triangle formed by t i , t j and t k . This is expressed by the following three orientation predicates: d t i t j p t j t k p t k t i p The proof by induction on the length of a convex polygon S 0 encompassing p works as follows. First, S 0 cannot be empty, since it encompasses p. If S 0 is a singleton or a pair, the degenerate p p p will do. If the convex polygon contains at least three points, rst check whether p is equal to t 1 , t 2 , o r t 3 , then compute t 1 t 3 p. I f a n y of these conditions holds, then take t 1 t 2 t 3 . If none of these conditions holds, then p is actually inside the convex polygon de ned by t 1 t 3 : : : You can use an induction hypothesis on this smaller convex polygon to reach the result.

Once you have the three points encompassing p it is possible to show t h a t t h e coordinates of p are a positive linear combination of the coordinates of the three points. This is simply because the coordinates of p verify the following equation: p j t i t j t k j = t i j t j t k pj + t j j t k t i pj + t k j t i t j pj This equation is easily veri ed using the ring decision procedure.

If d t i t j p, t j t k p, a n d t k t i p hold, then the three determinants on the right-hand side are positive a n d b y Axiom 4 the determinant on the left-hand side also is.

Dividing by this determinant is legitimate and p really is obtained by a p o s i t i v e linear combination of the t's. The cases where the convex hull contains less than three points are also easy to prove. Notice that this proof requires that we w ork in a eld.

Finding initial data

Some algorithms require that one should be able to nd one or two points that are already present in the convex hull to start with. The package-wrapping algorithm is among them. The usual technique is to look for the minimal point using a lexical ordering on the coordinates.

If t is the minimal point for the lexical ordering on coordinates, and if p, q, and r are points of a triangle encompassing t, where the rst coordinates of p, q, and r are such p 1 r 1 then it is necessary that r 1 q 1 . Using this result again, we h a ve p 1 r 1 ) r 1 p 1 . T h us all three rst coordinates must be equal. But if they are, the points are aligned and cannot constitute a proper triangle.

Once we have the point t 2 whose coordinates are minimal for the lexicographic order, we can nd its left neighbour by c hoosing the only point that has the same second coordinate as t 2 . I f n o s u c h p o i n t exists, take the point t h a t i s maximal for the order de ned by: p q , c tpq That is transitive is ensured by Axiom 5 0 , using the point (t x t y + 1 ) as point s.

Proving algorithms

In this section we review the descriptions and proofs for the abstract setting: positions are general (no three points are aligned) and all numerical computation is hidden behind the orientation predicate and the function producing initial points. The incremental algoritm works by taking the points of the input data one by one and constructing the convex hull of all the points seen so far. At each step, either the new point already lies inside the current c o n vex hull and one can directly proceed to the next, or some edges need to be removed and two e d g e s need to be added. All tests rely on the orientation predicate.

The incremental algorithm

The edges from the previous convex hull are sorted into two categories, called blue and red edges (see gure 3). Blue edges represent edges, for which the new point lies on good side. In our description we a l s o s a y that an edge is blue if the new point is one of the extremities. Red edges are the other. When all edges are blue, the point is inside. All edges cannot be red at the same time and red edges are contiguous.

When looking for red edges there may be four cases:

1. no red edges (including no edges at all), 2. the edge between the last point of the list and the rst point of the list (the closing edge) is red, but the last edge of the list is not. 3. the closing edge is red and the last edge too, 4. the closing edge is blue. In the rst case, nothing needs to be done. In the second case, the list of edges can be viewed as the concatenation of two lists, l r and l b , s u c h that l r contains only red edges l b contains only blue edges, and the edge between the last element in l r and the rst element i n l b is red. This relies on the property that the list describing the convex hull has no duplicates. In this case, the result is the list p l b . T o prove this, we need to show t h a t p is not already in the list (easy), that the new point p is in the input data (easy), that p is left of all the edges in the convex hull (easy, since the edges that are kept are all blue, and the two n e w edges have p at the extremity), and that all the data considered so far is left of the edges tp and pt 0 , w h e r e t and t 0 are the purple points, the last and rst elements of l b , respectively. The point t is di erent from t 0 if there is a red edge.

Let us detail part of the proof for the other cases (see gure 4). Since t and t 0 are distinct, there exists a t k distinct from t that is last-but-one in l b . L e t q be the other extremity of the rst red edge, and let r be an arbitrary point among those that have already been processed. By The two other cases are solved by rotating the convex hull to return to the second case. Rotating a convex hull is simply decomposing this convex hull as the concatenation of two lists. In the following, we will denote concatenation of two lists l 1 and l 2 as l 1 l 2 . Since concatenation is a simple structural recursive function (recursive o ver the rst argument ) , i t i s a s i m p l e m a t t e r o f s t r u c t u r a l recursion over l 1 to show that rotation preserves inclusion in another list, that it preserves the absence of duplicate elements, and that it preserves the data encompassed. Thus, if l 2 l 1 is a convex hull of l, l 1 l 2 also is.

For the third case, one traverses down the list representing the convex hull until one nds the rst point t j such t h a t t j;1 t j is blue and t j t j+1 is red. In that case the rotation is given by l 2 = t j t and l 1 = t 1 t j . I f t i is the rst point in l 1 such that the edge t i;1 t i is red and the edge t i t i+1 is blue, then one has: l r = t j t (t 1 t i;1 ) and l b = t i+1 t j .

The fourth case is symmetric to the third one. One also needs to nd the two purple points that are boundaries between blue and red edges but they are in reverse order.

The main algorithm repeats this computation for each point. It is programmed as a simple structural recursion over the list of inputs.

The package-wrapping algorithm

The package-wrapping algorithm relies on the possibility to nd an initial edge that is known to be part of the convex hull, then it proceeds by constructing an open encompassing list l of pairwise distinct points, included in the data, nding the right neighbour of the rst element in the list at each s t e p . A t e a c h step one compares the right neigbour with the last element of the list. When the two points are equal, the algorithm terminates without adding the current list. When the two points are di erent, the right neigbour is added to the list and recursion restarts (see gure 5). During the recursion the list alway s h a s a t l e a s t future point t p Fig. 5. Intermediary iteration for the package-wrapping algorithm. two elements, let t be the head and s be the second element. Finding t's right neighbour is done by nding the greatest element for the relation de ned by p q , c tpq:

This relation is not re exive, but its re exive closure is a total order, thanks to Axioms 2, 3, and 5. Finding the least element of a list for an arbitrary order is an easy piece of functional programming, and it is also easily proved correct. If p is this greatest element, c tqp holds for every point in the input data, and by Axiom 1 c ptq also holds. Thus the input data is left of the edge pt and the list p l also encompasses the input data.

A k ey ingredient to proving the termination of the algorithm is to ensure that the list being constructed never contains duplicate elements and is included in the input data. The hardest part is to show t h a t c hecking that p is distinct from the last element o f l is enough to ensure that p does not belong to l. In fact, if t is the head of the list, t 0 is the last element a n d q is another element of the list, one needs to show that d tqt 0 holds. This is done by recursion over the length of list l but there are many i n tricacies with base cases, because the proof relies on Axiom 5, which needs 5 distinct point s . L e t u s h a ve a look at this proof.

First assume q is the second element of the list, d tqt 0 comes easily from the fact that the list encompasses the input data and t 0 is in the input data. Then if q is not the second element of the list, let t 1 be this second element. We d o t h e proof by structural recursion over the list q t 0 . If this list is reduced to 0 or 1 element, this is contradictory with our hypotheses, if this list has two elements, then qt 0 is an edge from the list, and d qt 0 t holds because t is in the input data. The result then comes from Axiom 1. If the list q t 0 has more than two elements, let q 0 be the second element, one has d tq 0 t 0 by induction hypothesis, d tqq 0 becauseqq 0 encompasses the input data and thanks to Axiom 1. Moreover, the properties d tt 1 q, d tt 1 q 0 , and d tt 1 t 0 hold, because q, q 0 , and t 0 are in the input data. This is enough to use axiom 5 and conclude. Now, the package-wrapping algorithm is not structural recursive. The termination argument is that the size of the list l is bounded by the size of the input data, because all its elements are pairwise distinct, and that this size increases at each step. For this reason, the function is de ned using Coq's solution for de ning well-founded recursive functions. This means that we h a ve to exhibit a set containing the principal argument of the recursive function and a well-founded order on this set and show that recursive calls are done only on smaller terms.

The set L S we exhibit is the set of lists that are included in the input data S and contain no duplicates. In Coq's type theory, this set is described as a type using dependently typed inductive types. We show that the length of a list in in this set is always smaller than the size N of the input data S. W e de ne the order < S to be : l < S l 0 , N ; length(l) < N ; length(l 0 ):

Proving that this order is well-founded is a simple application of theorems about composing well-founded orders with functions provided in Coq's libraries. In other theorem provers where the termination of recursive functions is ensured by exhibiting measure functions, the measure function to exhibit is the following one:

f : l 7 ! N ; length(l):

At each recursive step, the main argument is l and the argument for the next recursive c a l l is p l such that p l 2 L S , therefore one has length(p) < length(p l) N and f(p l) < f (l).

If we w ant t o a l l o w for aligned points in the input data, we need to review our technique. One solution is to change Axiom 3, since it is the one that assumes that a triple of points can always be viewed as a triangle. If we do that, any decision taken with respect to orientation in the algorithms must be reviewed and replaced by a n e w decision taking extra cases into account. The structure of the algorithm remains, but all the details of its control-ow are deformed by this constraint.

A second solution is to change the orientation predicate so that it satis es all ve axioms even in the presence of aligned triples. This solution somehow imposes that one changes the de nition of convex hulls, since our de nition relies on the orientation predicate. We shall see that both solutions actually require that we relax our de nition of a convex hull.

Segments

When three distinct points are aligned, one of them lies between the two other.

We i n troduce a new predicate, written epqrd to express this. Axiom 3 is changed to obtain the following formulation: Axiom 3 8p q r: p = q _ q = r _ p = r _ d pqr _ d prq_epqrd_eqrpd_erpqd:

In Knuth's work, Axioms 1-5 express a form of internal consistency for the orientation predicate. We need to have the same kind of internal consistency for the segment predicate and between segments and oriented triangles. The rst two axioms describe internal consistency for three aligned points.

Axiom 6 8p q r: epqrd)erpqd. Axiom 7 8p q r: epqrd) :eprqd.

The next four axioms describe consistency between segments and triangles.

Axiom 8 8p q r t :d pqr^eptqd) c ptr.

Axiom 9 8p q r t :d pqr^epqtd) c ptr.

Axiom 10 8p q r t :d pqr^eptqd) c tqr.

Axiom 11 8p q r t :d pqr^etpqd) c tqr.

The last two axioms describe internal consistency of segments for four aligned points.

Axiom 12 8p q r t :eqrtd^epqtd)eprtd. Axiom 13 8p q r t :eqrtd^epqrd)eprtd.

With these new notions, we c hange the speci cation of a convex hull for a set S to be a list of points t 1 : : : t k such that for every p in S and every t i , t j such that t j = t i+1 or t j = t 1 if t i = t k one has d t i t j p_et i t j pd.

The working horse predicate

In the new setting we use a new predicate that combines orientation and segment membership. For the rest of the paper, we will denote this predicate pqr, w i t h the following meaning: pqr , d pqr _ r 2 f p qg_eprqd:

Generalizing our algorithms to the degenerated cases works by replacing the orientation predicate by the new predicate wherever it occurs and adding a few equality tests to take care of degeneracies. What matters really is that equivalents of Axiom 5 still hold: 8p q r s t :tsp ^tsq ^tsr : qtp^qtr) ptr 8p q r s t :stp ^stq ^str : tqp ^tqr ) tpr:

Thus, if t is known to be on the convex hull, we k eep the property that the order 0 can be used to nd minimal and maximal elements, when this order is de ned as follows: p 0 q , tpq:

Perturbation

An alternative approach i s t o c hange the orientation predicate so that Axiom 3 always holds. This option seems awkward: if we c hange the orientation predicate how can we b e s u r e t h a t w e actually compute the same convex hull? The answer is to make sure that we c hange the orientation predicate in a manner that remains consistent with the general orientation predicate, so that the convex hull computed in non-degenerated cases is unchanged.

The solution we are going to present is based on the idea that if points are aligned, then moving one of the points only slightly will remove the degeneracy. This make it possible to compute an approximation of the hull. However, points are not actually moved, and all points that should appear on the convex hull will be guaranteed to appear in the approximation we compute. Imprecisions will occur only for points that are on the segment b e t ween two legitimate vertices of the convex hull. For these points, this method may o r m a y not include them in the resulting convex hull.

The solution is to consider that all points are moving in a continuous manner and that the con guration that is being studied is the snapshot at date t = 0 .

The determinants used to compute the orientation predicate also are continuous functions in t, that may h a ve t h e v alue 0 when t = 0 . The orientation predicate is continuous for all unaligned triples. However, if the perturbation function can be viewed as a polynomial function in t it is possible to ensure that lim t!0+ is either positive or negative. If the points are not aligned, the limit is the same as the value in 0, if there are aligned points, taking a value of t that is small enough will always return the same sign, so that we d o n o t e v en have to predict the value of t: the sign can be predicted from the signs of the polynomial coe cient.

We use a function taken from 1] to indicate how a l l p o i n ts move. The coordinates of the point p actually are given by the formula x p +t y p y p +t (x 2 p +y 2 p ).

The orientation predicate for points x, y, z then becomes the following term:

x p + t y p y p + t (x 2 p + y 2 p ) 1

x q + t y q y q + t (x 2 q + y 2 q ) 1

x r + t y r y r + t (x 2 r + y 2 r ) 1

= D 1 + D 2 t + D 3 t 2
where D 1 , D 2 , a n d D 3 are the determinants de ned as follows:

D 1 = x p y p 1 x q y q 1 x r y r 1 D 2 = y p x 2 p + y 2 p 1 y q x 2 q + y 2 q 1 y r x 2 r + y 2 r 1 D 3 = x p x 2 p + y 2 p 1

x q x 2 q + y 2 q 1

x r x 2 r + y 2 r 1: The proof that all three determinants cannot be zero is formalized using the sketch provided in 1] and relies very much on deciding equalities in a ring, plus a lot of painful reasoning on equalities. However, we w ork on polynomials and all our proofs are done in an ordered archimedian eld: the eld of rational numbers su ces so that our model is adequate for arithmetic numbers as they are implemented in computers.

If the rst derminant is 0, then the three points are aligned. There exists three numbersm 1 , m 2 , a n d m 3 such t h a t m 1 y = m 2 x + m 3 is the equation of the line containing the three points. As these points are distinct, m 1 and m 2 cannot both be null. There are two cases, depending on whether m 1 = 0 .

If m 1 = 0, then x p = x q = x r , the second determinant is equal to the following determinant: y p x 2 p + y 2 p 1 y q x 2 q + y 2 q 1 y r x 2 r + y 2 r 1 = y p y 2 p 1 y q y 2 q 1 y r y 2 r 1 = ( y p ; y q )(y q ; y r )(y r ; y p ):

(1) This value can only be zero if one of the factors is 0, i n w h i c h case two o f t h e points have both coordinates equal: they are the same.

If m 1 = 0 , D 2 = 0 if and only if m 1 D 2 = 0 . m 2 1 y p x 2 p + y 2 p 1 y q x 2 q + y 2 q 1 y r x 2 r + y 2 r 1

= ( m 2 1 + m 2 2 )
y p y 2 p 1 y q y 2 q 1 y r y 2 r 1 = ( m 2 1 + m 2

2 )(y p ; y q )(y q ; y r )(y r ; y p )

(2) Here again, at least two of the points must have the same vertical coordinate for the determinant t o b e 0.

The same reasoning applies to the third determinant, which is symmetric with the second one, replacing x coordinates with y coordinates. Now for every triple of points, the values D 1 , D 2 , and D 3 cannot be 0 at the same time. We can then show that for every triple p q r in the input data, there exists a number pqr such that the perturbed determinant n e v er becomes 0 in the interval ]0 pqr . Let abs(x) be the absolute value of x. I f D 1 is non null, then we can choose pqr to be the following value: ):

If D 1 = 0 and D 3 = 0 we can take a n y positive v alue for pqr . This concludes the proof of existence of pqr. T h i s p r o o f w as also made only using axioms that are ful lled by the eld of rational numbers.

We could now replace the orientation predicate used in the previous section by a new orientation predicate, noted pqr and de ned in the following manner: pqr = x p y p 1 x q y q 1 x r y r 1 > 0 _

x p y p 1 x q y q 1 x r y r 1 = 0 ^ y p x 2 p + y 2 p 1 y q x 2 q + y 2 q 1 y r x 2 r + y 2 r 1 > 0 _

x p y p 1 x q y q 1 x r y r 1 = y p x 2 p + y 2 p 1 y q x 2 q + y 2 q 1 y r x 2 r + y 2 r 1 = 0 ^

x p x 2 p + y 2 p 1

x q x 2 q + y 2 q 1

x r x 2 r + y 2 r 1 = 0

Note that the pqr are never computed, we o n l y s h o w their existence to show t h a t the computation of these two determinants is enough to predict the orientations for a general position that is arbitrarily close to the input data. It is not necessary to return a convex hull composed of perturbed points: the points from the input data are correct, since they can be shown to be arbitrarily closed to a legitimate convex hull (we h a ve not done this proof yet).

It is also interesting that the convex hull computation done in this manner returns exactly the same result as the computation described in previous section if the input data is not degenerated. If the input data is degenerated and the convex hulls are required to contain no aligned points, a trivial traversal of the list of points may be needed to remove the useless points. The cost of this traversal should be proportional to the total cost. It is remarkable that testing the sign of the D 2 can be done by simply combining the signs of the terms (y p ;y q ), (y q ;y r ), (y r ;y p ), according to equations (1) and (2). In fact, D 2 is strictly positive if and only if the line carrying the three points is not horizontal, the point r is outside the segment pq, and the points p and q are ordered in their second coordinate, similarly for D 3 . As a result, testing the sign of D 2 and D 3 is equivalent to testing whether r) is outside the segment pq and whether p and q are ordered lexicographically (using the second coordinate rst). Thus, we h a ve a program that computes (pqr) without computing the larger determinants above. This algebraic reasoning brings an elegant justi cation to the use of segment predicates and lexical ordering in combination with traditional orientation.

Conclusion

The development of certi ed algorithms is well suited to domains where it is easy to have a mathematical formulation of the data and the speci cations. Computational geometry is one of these domains. Euclidean geometry makes it easy to connect the various notions encountered in computational geometry to real analysis and algebra. In our experiment, the presence of a decision procedure for ring equalities has been instrumental, although it has also been quite annoying that this decision procedure was not complete, at least for the version we used. This decision procedure is the Ring tactic, a tactic based on a re ective approach 4].

The problem of numeric computation also deserves some thoughts. We h a ve performed our proofs with an axiomatized s e t o f n umbers for the coordinates. The axioms we have taken were carefully chosen to make them accurate as a representation of computer numbers used in exact computation. These numbers live in an archimedian eld. The model for the eld that is needed to express the correction of the algorithm is the eld of rational numbers. However, the algorithm does not actually use numbers in the whole eld. It only computes a polynomial expression of degree two, so that if the numbers describing input data are bounded, the set of numbers needed for the algorithm to work is much smaller. It could be useful to express the di erence between the concrete arithmetic used in the algorithm, and the abstract arithmetic used in the proof. This is actually the spirit of the work of Laurent th ry and his colleagues in a project to study computer arithmetics.

Perturbation methods are used in other parts of computational geometry, for instance for Delaunay triangulation, where degenerate positions occur when four points are cocyclic. They seem very well adapted as a method to make algorithms robust to degeneracy without changing their structure. Still, it seems a rare and lucky occurrence that perturbations are used in this paper to justify an algorithm that actually does not use any perturbation, since it only boils down to an extension of the orientation predicate with a segment predicate and lexical ordering.

We actually heard about the perturbation method after having done most of the proofs for the algorithm based on the extended axioms combining the traditional orientation and segment predicates. Still both developments are useful in their own right, since the rst development can be tuned to include or exclude colinear points from the convex hull at will. For the second method, we can predict that colinear points on edges of the convex hull that grow lexicographically will be included in the convex hull, while colinear points on edges that decrease lexicographically will not be included. This information that we h a ve not formally proved at the time of writing this article, may b e used when writing a second pass that removes colinear points, but it won't be easy to recover colinear points that have already been discarded by the convex hull computation.

We w ould like t o a c knowledge the help of several researchers from the computational geometry group at INRIA Sophia Antipolis. In particular, Mariette Yvinec pointed us to Knuth's work on axioms and Olivier Devillers described us the perturbation method.

Fig. 1 .

 1 Fig. 1. Axiom 4: Consistent orientation for four triangles.

Fig. 3 .

 3 Fig. 3. Red and blue edges in the incremental algorithm

Fig. 4 .

 4 Fig. 4. Using axiom 5 0 to justify adding an edge.

  or D 3 = 0 , just ignore the corresponding term in the minimal computation. If D 1 = 0 and D 3 6 = 0 , t h e n w e c a n c hoose pqr to be the following value: pqr = min(1 abs( D 2 D 3