domain: sde sdu

text:
Amri Abdellah
email: abdellahamri7@gmail.com

Zakaria Saâdi
email: zakaria.saadi@irsn.fr

Rachid Ababou
email: rachidababou@gmail.com

Modeling two-phase flow with hysteresis: comparative study of hysteresis models and application

Keywords: Porous media, Immiscible two-phase flow, Hysteresis, Capillary entrapment, Wetting-Drying cycles, iTOUGH2 code simulations, Radioactive waste disposal, Water-Gas flow, CO2 storage, Laboratory scale experiment

Drainage (drying) and imbibition (wetting) cycles are expected in many applications involving subsurface flows, such as underground CO2 storage, hydrocarbon extraction and production, and deep geological disposal of radioactive waste (particularly in the potential production of hydrogen due to corrosion and radiolysis). In the presence of imbibition/drainage cycles, the two-phase constitutive relations between water saturation, capillary pressure and relative permeability to liquid and gas, cannot be represented by simple one-to-one functions: they depend on the history of wetting/drying processes. Recent numerical and experimental studies highlight the importance of hysteresis and its effect on two-phase flows in porous media. In this study, an extensive critical review of conceptual hysteretic models is first presented. The performance of these models is evaluated and tested from a theoretical point of view based on their assumptions, and by confronting them with measured data for various porous materials extracted from the literature (cement, sands, Bentonite MX80). The advantages and disadvantages of the different hysteretic constitutive models are discussed. Finally, we implement the five hysteresis models in the iTOUGH2 code, and we develop a numerical simulation of a two-phase flow experiment with immiscible Modeling two-phase flow with hysteresis fluids at the laboratory scale (meters), in order to compare the flow results obtained with or without hysteresis, and also, with or without capillary entrapment. The results highlight the need to take into account hysteresis in the modeling of multiphase porous media flows, especially when they are subjected to drying-wetting cycles. It is suggested that the model of Beriozkin and Mualem 2018 should be preferred both from a theoretical standpoint and for its computational advantages.

Highlights

• Quantitative assessment of several significant and representative conceptual models for hysteretic two-phase flow in porous media. • Identification of best hysteresis model in terms of physical parametrization, accuracy, and ease of implementation. • Gas entrapment in hysteresis models is found to be essential for accurate simulations of gas-liquid flow in porous media.

Introduction

Accurate characterisation of the relationships of hydraulic properties is essential to understanding and modeling multi-phase flow in porous media. These relationships are subject to the hysteresis phenomenon during drying and wetting cycles. Various models have been used over the years to describe the hysteretic capillary-saturation relationships. A detailed overview of existing hysteresis models can be found in [START_REF] Mualem | General scaling rules of the hysteretic water retention function based on Mualem's domain theory[END_REF], [START_REF] Zhang | A review and statistical study of existing hysteresis models for cementitious materials[END_REF] and [START_REF] Albers | Modeling the hysteretic behavior of the capillary pressure in partially saturated porous media: a review[END_REF]. Many studies [START_REF] Lenhard | A model for hysteretic constitutive relations governing multiphase flow: 3. Refinements and numerical simulations[END_REF], [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF], [START_REF] Zhang | Modelling of sorption hysteresis and its effect on moisture transport within cementitious materials[END_REF], [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF]) have highlighted the necessity of taking into account hysteresis phenomena in the modeling of multiphase transport in porous media.

For example, Parker, Lenhard (1987) and Lenhard, Parker (1987) produced two successive articles dealing with hysteresis in the capillary water retention curve P c (S) and in the liquid-gas relative permeability curves k rL (S) and k rG (S). A hysteretic model based on the concept of non-wetting phase entrapment was proposed, using the [START_REF] Scott | Hysteretic effects on net infiltration // Advances in infiltration[END_REF] scaling technique. After a series of publications, [START_REF] Lenhard | Measurement and modeling of three-phase saturation-pressure hysteresis[END_REF] showed that there is a significant difference in fluid saturations between hysteretic and non-hysteretic simulations of NAPLs (Non-Aqueous Phase Liquids) distribution in the vadose zone. [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] introduced the concept of non-wetting phase entrapment saturation in the hysteresis model of [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF] in a simple manner, introducing the formula of [START_REF] Land | Calculation of imbibition relative permeability for two-and threephase flow from rock properties[END_REF] in the unsaturated Van Genuchten-Mualem ("VGM") formulation [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1[END_REF]). The model was extended to gas and water permeabilities using the Parker, Lenhard (1987) formulation. The [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] hysteretic model was implemented in the iTOUGH2 code, and illustrated with an example simulation of injection of supercritical carbon dioxyde in a saline aquifer, where the effect of hysteresis is well demonstrated. [START_REF] Zhang | A review and statistical study of existing hysteresis models for cementitious materials[END_REF] conducted a comparative study of several models, including the empirical models of [START_REF] Rubin | Numerical method for analyzing hysteresis-affected, post-infiltration redistribution of soil moisture[END_REF], [START_REF] Feng | Hysteretic influence associated with thermal conductivity sensor measurements // Proceedings from Theory to the Practice of Unsaturated Soil Mechanics in Association with the 52nd Canadian Geotechnical Conference and the Unsaturated Soil Group[END_REF], [START_REF] Braddock | Application of a soil water hysteresis model to simple water retention curves // Transport in porous media[END_REF], [START_REF] Li | Modelling of hysteresis response for arbitrary wetting/drying paths // Computers and Geotechnics[END_REF], and the conceptual models of [START_REF] Mualem | A conceptual model of hysteresis[END_REF] and [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF]. They showed that the conceptual models are more accurate than the empirical ones for predicting wetting and drying cycles of moisture dynamics in cementitious materials.

On the other hand, [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF] developed a hysteresis model based on relating macroscopic variables to the microscopic properties of the porous media, including void distribution and connectivity function. Using an intermediate-scale laboratory experiment of non-wetting fluid injection in a porous medium initially saturated with a glycerol-water mixture, they concluded that their hysteresis model predicts well the non-wetting plume migration and its distribution during and after injection, compared to the nonhysteretic model. More recently, [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] showed that the Parker, Lenhard (1987) model for the capillary water retention curve P c (S l ) has the undesirable property of "non parallelism" . That is, the two scanning curves departing from the main drying curve and the first drying curve are not "parallel". Consequently, they are not accompanied by the same increase in saturation. In other words, the two primary wetting processes, predicted by the Parker, Lenhard (1987) model are described by non-parallel effective saturation curves, displaying different increments of water content during the same changes of capillary pressure. Thus, the Parker, Lenhard (1987) model introduces an inconsistency into the apparent saturation model, although on the other hand, it is practical because it uses an available computational tool (i.e. Land's 1968 formula). On the other hand, [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] used the linear law for non-wetting phase entrapment applied to the [START_REF] Mualem | A conceptual model of hysteresis[END_REF] model to derive 6 equations of scanning curves stemming from the main curves.

In conclusion, choosing the appropriate model of hysteresis (and of nonwetting phase entrapment) appears to be a challenging task in modeling multiphase flow in porous media.

In the present study, we perform a review and comparison between the well-known physically-based conceptual models in the literature. The models are tested on a wide range of porous media. The representation of non-wetting phase entrapment is discussed, and a general algorithm for modeling hysteresis is provided, with a discussion of the important hyper-parameters. Finally, the intermediate-scale laboratory experiment of Soltrol-220 injection into a slab of sandy porous medium fully saturated with a glycerol-water mixture [START_REF] Trevisan | Investigation of mechanisms of supercritical CO2 trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions[END_REF]) is simulated with and without hysteresis.

2 Brief description of some conceptual hysteresis models

Hysteresis models can be classified into two categories, empirical and conceptual. The empirical models are based on fitting hyper-parameters from the two scanning curves in order to predict any scanning curve. These models are not taken into account in this study due to pumping errors (PE) resulting from the non-closure of scanning hysteresis loops associated with spurious cyclic oscillations as reported in [START_REF] Zhang | A review and statistical study of existing hysteresis models for cementitious materials[END_REF]. These empirical models may require an additional measured scanning curve of higher order for their validation.

The conceptual models are generally based on independent and dependent domain theory which has been first introduced in porous media flow by [START_REF] Poulovassilis | Hysteresis of pore water, an application of the concept of independent domains[END_REF] . This theory has been developed by [START_REF] Preisach | Über die magnetische Nachwirkung[END_REF] to describe ferromagnetic hysteresis. Conceptual models have been extended in many recent studies to take into account the capillary entrapment phenomena. In this section we present an exhaustive critical review of a wide range of conceptual hysteresis models based on [START_REF] Mualem | A conceptual model of hysteresis[END_REF], [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF], [START_REF] Haverkamp | Soil water hysteresis prediction model based on theory and geometric scaling // Environmental mechanics: Water, mass and energy transfer in the biosphere[END_REF], [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF], and [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF].

Empirical hysteresis models without non-wetting phase entrapment

The classical models of hysteresis are based on interpolating values of scanning curves from main curves. In this review we focus on hysteresis models derived from two theories. The first one is the domain theory which combines models based on the dependent domain theory from [START_REF] Mualem | A conceptual model of hysteresis[END_REF] and the independent domain theory from [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF]. The second theory is based on the rational extrapolation of [START_REF] Parlange | Capillary hysteresis and the relationship between drying and wetting curves[END_REF], also adopted later by [START_REF] Haverkamp | Soil water hysteresis prediction model based on theory and geometric scaling // Environmental mechanics: Water, mass and energy transfer in the biosphere[END_REF].

Mualem 1974 model

According to the domain theory in soil physics, the porous media can be presented by a bi-variate pore water distribution function f (r, ρ) where : r = r(P c ) -r(P c,min ) r(P c,max ) -r(P c,min ) and ρ = ρ(P c ) -ρ(P c,min ) ρ(P c,max ) -ρ(P c,min ) are two normalized pore radius parameters varying from 0 to 1 as capillary pressure P c varies from P c,min to P c,max . They can be interpreted respectively as the radius of pore body and the radius of pore neck, following representation of hysteresis phenomena as the ink-bottle effect as described in M 'Jahad (2012). Following this representation of pore size distribution, the water saturation is obtained by integration of the pore water distribution function, as follows:

S(P c ) = r∈I r(Pc ) ρ∈I ρ(Pc ) f (r, ρ) dρ dr
(1) [START_REF] Mualem | A conceptual model of hysteresis[END_REF] assumed that each group of pores is independent of its neighbors. More precisely, this assumption stipulates that pore neck size is totally independent from pore body size and the bi-variate function can be expressed as the product of two independent functions f (r, ρ) = h(r) × l(ρ). The water saturation S(P c ) at a given capillary pressure can be obtained by integration of pore water distribution function over specific intervals, I ρ(Pc) and I r(Pc) , as follows:

S(P c ) = r(Pc) 0 ρ(Pc) 0 f (r, ρ) dρdr = r(Pc) 0 ρ(Pc) 0 h(r) × l(ρ) dρdr (2)
In order to illustrate the drying and wetting processes, and to derive scanning curves equations, [START_REF] Mualem | A conceptual model of hysteresis[END_REF] In these two diagrams, the process of drying and wetting of pores is governed by two assumptions:

• For the drying process from

P c ( R) to P c ( R -d R) only the groups of pores characterized by R ≤ r ≤ 1 and R -d R ≤ ρ ≤ R are drained, • For the wetting process from P c ( R) to P c ( R + d R) all the groups of pores characterized by R ≤ ρ ≤ R + d R are filled.
Based on these assumptions, [START_REF] Mualem | A conceptual model of hysteresis[END_REF] derived an explicit formulation of any scanning curve order (first drying, second wetting, third drying..). Let S N,d (P c ) (resp. S N,w (P c )) be the unknown saturation for the Nth drying (resp. Nth wetting) curve from the Nth turning point S tp,N (P c,N ), then S N,d (P c ) (resp. S N,w (P c )) can be expressed by:

S N,d (P c) = S tp,N (P c,N ) - 1 r(Pc) ρ(Pc) ρ(P c,N ) h(r) × l(ρ) dρdr (3) S N,w (P c) = S tp,N (P c,N ) + 1 0 ρ(P c,N ) ρ(Pc) h(r) × l(ρ) dρdr (4)
From the integration of the pore water distribution based on equations ( 3) and ( 4), the Nth scanning path for the drying and wetting process can be written explicitly as :

S N,d (P c ) = S tp,N (P c,N ) -(S w (P c,N ) -S w (P c ))( Smax,w-S d (Pc) Smax,w-Sw(Pc) ) if P c ≤ P c,N -1 S N -2,d (P c ) if P c ≥ P c,N -1 (5) S N,w (P c ) = S tp,N (P c,N ) + (S w (P c ) -S w (P c,N ))( Smax,w-S d (P c,N ) Smax,w-Sw(P c,N ) ) if P c ≥ P c,N -1 S N -2,w (P c ) if P c ≤ P c,N -1 (
6) where S 0,w = S w (resp S 0,d = S d ) is the main wetting (resp drying) curve, and S max,w is the maximum water saturation for a complete wetting process starting with an initially dry porous medium.

Equations ( 5) and ( 6) define scanning curves of any order as illustrated in figure (3), where a series of scanning curves were plotted based on the hysteretic model of [START_REF] Mualem | A conceptual model of hysteresis[END_REF] and the Van Genuchten (VG) formulation (Eq. ( 7), Van Genuchten (1980)) for the main curves with parameters n w = n d =2.145 , α d =2.61×10 -4 [P a -1 ], α w =9.16×10 -4 [P a -1 ], S ls =1.0, and S lr =0.

S(P

c ) -S lr S ls -S lr = [1 + (αP c ) n ] -m (7)
The

Mualem 1984 model (dependent domain theory)

The previous model of [START_REF] Mualem | A conceptual model of hysteresis[END_REF] assumed that the drying and wetting processes in pores are independent of the state of saturation of neighboring pores, which is not the case in a real porous medium. To take into account this dependency, [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF] introduced different factors in order to take into account the portion of drainable pores in a new model, called the universal model, corresponding to the 'dependent domain' theory of [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF]. This new model integrates (involves) only one dependent factor P d which can be related to the pore blockage against the non-wetting phase (air).

The model proposes the following formulation for the N-th wetting (resp, N-th drying) curve S N,w (P c ) (resp S N,d (P c )) departing from the N-th turning point S tp,N (P c,N ). First, we have:

S N,w (P c ) = S tp,N (P c,N ) + P d (S tp,N (P c,N )) × (S w (P c ) -S w (P c,N ))( Smax,w-Sw(P c,N ) Smax,w ) S N -2,w (P c ) (8)
where the first equation is for P c ≥ P c,N -1 , and the second equation is for P c ≤ P c,N -1 . And secondly, we have:

S N,d (P c ) = S tp,N (P c,N ) -P d (S N,d (P c )) × (S w (P c,N ) -S w (P c ))( Smax,w-Sw(Pc) Smax,w ) S N -2,d (P c) (9)
where, again, the first equation is for P c ≤ P c,N -1 , and the second equation is for P c ≥ P c,N -1 . In these equations, P d (S) is defined as:

P d (S) = (Smax,w-S)×Smax,w (Smax,w-Sw(P * c )) 2
and P * c is defined by the relation:

S d (P * c ) = S.
Note that equation ( 9) provides an implicit solution for the N-th drying curve; therefore, an iterative method is required to compute S N,d (P c ). Comparing the independent and dependent domain theories of [START_REF] Mualem | A conceptual model of hysteresis[END_REF] and [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF] as illustrated in figure (3) and figure (4), respectively, it can be seen that, for the same turning points, the dependent domain theory leads to a higher desaturation than the independent domain theory.

Haverkamp et al. 2002 geometrical scaling model

Based on the theoretical model of [START_REF] Parlange | Capillary hysteresis and the relationship between drying and wetting curves[END_REF], [START_REF] Haverkamp | Soil water hysteresis prediction model based on theory and geometric scaling // Environmental mechanics: Water, mass and energy transfer in the biosphere[END_REF] derived geometrical scaling conditions for defining a hysteresis model. The model assumes that all the main curves and the scanning curves have the same normalized Van Genuchten formulation (Eq. ( 7)). For the Nth scanning curve, the normalized water saturation can be defined as:

SN,p (P c ) = S N,p (P c ) -S lr,N,p S ls,N,p -S lr,N,p = [1 + (α N,p × P c ) n N,p ] m N,p (10)
where SN,p (P c ) (resp S N,p (P c )) is the normalized water saturation (resp. the actual water saturation) for the Nth scanning curve with respect to the drying path (p=d) and wetting path (p=w). Note here that the symbol "p" is used as a path label (drying or wetting). Parameters S ls,N,p , S lr,N,p , n N,p , m N,p , α N,p are the VG-parameters corresponding to the N-th scanning curve with respect to the path "p". The determination of the Nth scanning curves (S ls,N,p , n N,p , m N,p , α N,p ) from the main curves parameters (S ls,0,p , n 0,p , m 0,p , α 0,p ) is based on the constitutive relationships as follows...

• All wetting and scanning curves have identical shape parameters:

n = n 0,w = n 0,d = n N,w = n N,d (11a)
m = m 0,w = m 0,d = m N,w = m N,d (11b)
• The pressure scale parameters of the main wetting and the main drying curves are related by :

α 0,w = k × α 0,d where k = 2 if m × n > 1 [1 + m × n] 1 m×n if 0 < m × n < 1 (12)
• The closure of scanning curves implies that each curve joins the point of departure of its parent curve, which leads to the specific conditions S ls,0,d = S ls,0,w = S ls,1,w = S ls and S lr,0,d = S lr,0,w = S lr,1,d = 0.

• The pressure scale parameter for the N-th scanning curve is given by:

α N,p = α 0,p × [ S ls S ls,N,p -S lr,N,p ] -1 m×n (13)
• The N-th scanning curve is defined in the range of the actual turning point S tp,N and the previous turning point S tp,N -1 . For example, for the first wetting process we have S tp,N is the turning point from the main drying curve and S tp,N -1 = 1. For the Nth wetting process departure from (S tp,N ,P c,tr,N ) and arrival to (S tp,N -1 ,P c,tr,N -1 ) for the path p process, where p=d for drying and p=w for wetting, the VGM parameters S lr,N,p and S ls,N,p are derived from the constitutive equations :

S ls,N,p -S lr,N,p = S tp,N -S tp,N -1 S * N,p (P c,tr,N ) -S * N,p (P c,tr,N -1 ) (14a) S * N,p (P c ) = [1 + k ϵp × ( S ls,N,p -S lr,N,p S ls,N -1, p -S lr,N -1, p ) 1 m × (α k-1, p × P c ) n ] -m (14b) S lr,N,p = S tp,N -1 × S * N,p (P c,tr,N ) -S tp,N × S * N,p (P c,tr,N -1 ) S * N,p (P c,tr,N ) -S * N,p (P c,tr,N -1 ) (14c)
where ϵ p =1 for the wetting path (p=w) and ϵ p =-1 for the drying path (p=d).

The p parameter designates the previous path of a scanning curve p (if p=d then p=w and if p=w then p=d). The system of equations ( 14 The geometrical scaling of [START_REF] Haverkamp | Soil water hysteresis prediction model based on theory and geometric scaling // Environmental mechanics: Water, mass and energy transfer in the biosphere[END_REF] gives an explicit expression for any scanning curve order. However, the main wetting and drying curves must verify equation ( 13), which is not always justified. Furthermore, an additional computational cost is induced by this model, due to the need for saving the Van Genuchten parameters for each scanning curve, as well as for all turning points.

Non-wetting phase entrapment with hysteretic models

In multiphase porous media, subject to wetting and drying cycles, the nonwetting fluid may be present as a disconnected phase due to capillary forces. Several experiments highlight this phenomenon and its impact in applications like CO 2 trapping: Hesse, Woods (2010), [START_REF] Krevor | Capillary heterogeneity trapping of CO2 in a sandstone rock at reservoir conditions[END_REF][START_REF] Trevisan | Investigation of mechanisms of supercritical CO2 trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions[END_REF]. In this subsection, we discuss the different representations of non-wetting phase trapping and their implementation in different hysteresis models.

Non-wetting phase entrapment saturation

As reported through several experimental studies (e.g., [START_REF] Raeesi | Capillary Pressure Hysteresis Behavior of Three Sandstones Measured with a Multistep Outflow-Inflow Apparatus // Vadose[END_REF]) the amount of non-wetting phase entrapment saturation S ∆ nw,r is highly correlated to the initial non-wetting phase saturation S nw,i . Land's formula (equation ( 15)) is the most widely used in existing hysteresis models:

S ∆ nw,r = 1 1 S ∆ nw,max + 1 1-Stp + 1 1-S lr (15)
where "∆" stands for entrapment phenomenon and "nw,r" means non-wetting phase residual saturation. However, recent studies [START_REF] Suzanne | Experimental relationships between residual gas saturation and initial gas saturation in heterogeneous sandstone reservoirs[END_REF], [START_REF] Pentland | Measurement of Nonwetting-Phase Trapping in Sandpacks[END_REF]) reveal that the linear law of [START_REF] Aissaoui | Etude théorique et expérimentale de l'hystérésis des pressions capillaries et des perméabilités relatives en vue du stockage souterrain de gaz[END_REF] performs better to represent measured data. Aissaoui's formulation is an adaptation of the linear law of non-wetting phase entrapment. The non-wetting phase entrapment saturation S ∆ nw,r increases linearly with the initial gas saturation S g,i before stabilizing (leveling of) to a maximal value of non-wetting phase entrapment saturation S ∆ nw,r,max :

S ∆ nw,r = S ∆ nw,r,max Snw,c S nw,i if S nw,i ≤ S nw,c S ∆ nw,r,max if S nw,i ≥ S nw,c (16)
The main parameter to represent the relation between initial now-wetting phase saturation and residual non-wetting phase entrapment is the maximal value of the non-wetting entrapped saturation S ∆ nw,r,max . This parameter depends on the type of fluid used (oil/air-water) and on the micro-structure of the porous medium. [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] considered that S ∆ g,r,max is a constant material property that depends only on the porosity ϕ according to Holtz's formula taken from the petroleum literature: S ∆ g,r,max = -0.3136 × ln(ϕ) -0.1334 (17) [START_REF] Suzanne | Experimental relationships between residual gas saturation and initial gas saturation in heterogeneous sandstone reservoirs[END_REF] conducted a series of measurements on the maximum residual gas saturation for different porous media, and they concluded that the maximum gas entrapment saturation S ∆ g,r,max cannot be predicted using porosity and permeability.

Lenhard and Parker 1987 models

In this subsection, we present the models of Parker, Lenhard (1987) and Lenhard, Parker (1987) for a two-phase system liquid-gas as described in [START_REF] Lenhard | Comparing simulated and experimental hysteretic two-phase transient fluid flow phenomena[END_REF]. A detailed description of hysteresis in three phase systems can be found in [START_REF] Lenhard | A model for hysteretic constitutive relations governing multiphase flow: 3. Refinements and numerical simulations[END_REF] and [START_REF] Lenhard | Measurement and modeling of three-phase saturation-pressure hysteresis[END_REF]. In order to account for the non-wetting phase entrapment, Parker, Lenhard (1987) introduced the notion of apparent liquid (water) saturation S * l (P c ), which can be expressed as the sum of the effective saturation of liquid ( Sl ) and of the effective saturation of trapped gas ( Sent g ), as follows:

Sl * (P c ) = Sl (P c ) + Sent g (P c ) (18
) with:

Sl (P c ) = S l (P c)-S lr S ls -S lr Sent g (P c ) = S ent g (Pc) S ls -S lr
Based on a series of previous publications [START_REF] Lenhard | Measurement and modeling of three-phase saturation-pressure hysteresis[END_REF], [START_REF] Kaluarachchi | Multiphase flow with a simplified model for oil entrapment // Transport in porous media[END_REF]), [START_REF] Lenhard | Modeling multiphase fluid hysteresis and comparing results to laboratory investigations[END_REF] presented a hysteresis model that was aimed at avoiding the effect of "pumping errors" in the hysteretic loops. The scanning curves were predicted based on the hysteresis model of [START_REF] Scott | Hysteretic effects on net infiltration // Advances in infiltration[END_REF]. This new model applies the hysteresis model of [START_REF] Scott | Hysteretic effects on net infiltration // Advances in infiltration[END_REF] to the apparent liquid saturation S * l (P c ) and a linear representation of gas (non-wetting phase) entrapment in order to derive the effective liquid saturation Sl (P c ). The equation of the N-th scanning path "p" is determined by the most recent reversal point (S * tp,N , P * c,N ) and the preceding reversal point (S * tp,N -1 , P * c,N -1 ) as expressed by:

S * (P c ) = S * tp,N -1 + (S * tp,N -S * tp,N -1 ) S * p (P c ) -S * p (P * c,N -1 ) S * p (P * c,N ) -S * p (P * c,N -1 ) (19)
Lenhard, Parker (1987) assume that hysteresis in relative permeability relations is a result of non-wetting phase entrapment. The classical model of water and gas permeability of [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF] is expressed as follows:

k r,l ( Sl ) = Sτ l Γ(0, Sl ) Γ(0, 1) 2 (20) k r,g ( Sl ) = (1 -Sl ) τ ′ Γ( Sl , 1) Γ(0, 1) 2 (21)
where Sl = S l -S l,r 1-S l,r is the effective water saturation, τ and τ ′ are dimensionless parameters related to the pore water connectivity (for liquid and gas phase, respectively), and Γ is the incomplete beta function.

In addition, Lenhard, Parker (1987) imposes a correction to the original Mualem formulation of the relative permeability to liquid (Eq.( 20)) in order to take into account the non-wetting phase entrapment in the mobile fluid phases:

k r,l ( Sl ) = Sτ l Γ(0, Sl ) -Sgt 1-Stp × Γ( Stp , Sl ) Γ(0, 1) 2 (22) k r,g ( Sl ) = (1 -Sl ) τ ′ Γ( Sl , 1) Γ(0, 1) 2 (23)
where Sgt) is the effective trapped non-wetting phase saturation,

• Sgt = Sgt(S l -Stp) (1-S l,r )(1-Stp-
• Stp = Stp-S l,r
1-S l,r is the effective turning point saturation,

• Sl is the effective trapped non-wetting phase saturation, defined as the sum of the effective water saturation Sl and the effective trapped non-wetting phase saturation Sgt .

We note that equation ( 22) is obtained in a general form and may be used in the next section in order to derive an expression for water permeability by using Land's formula for non-wetting phase entrapment. Adaptation may also be made for another use of S ∆ g,r -S g,i relationships. The Lenhard, Parker (1987) model is a more realistic hysteresis model, it takes into account the phenomenon of non-wetting phase entrapment and its effect on the relative permeabilities. However, it involves intermediate operations to switch from apparent to effective water saturation. [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] applied the dependent domain theory of [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF] (see subsection 2.1.2) to a new formulation of Van Genuchten-Mualem (VGM) P c (S l ) relationship, accounting for entrapped gas saturation. The new formulation of the capillary pressure -liquid saturation relationship is described as follows :

Doughty 2007 model

P c (S l ) = P p 0 S l -S l,r 1 -S ∆ nw,r -S l,r -1 m p -1 1-m p (24)
where P p 0 , m p , S l,r are the VGM parameters of the main wetting (p=w) and the main drying (p=d). The amount of entrapped non-wetting phase saturation S ∆ g,r saturation depends on the the recent reversal point from drainage to imbibition S tp and on the input parameter of maximum residual gas saturation S ∆ g,r,max , which is a characteristic of the porous medium, calculated from the main wetting curve and linearly given by Land formulation (Eq. ( 15)). . [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] integrates equations ( 22) and ( 23) by using equation ( 24) in order to obtain explicit expressions for relative liquid permeability k l,r and relative gas permeability k g,r , respectively:

k r,l ( Sl ) = Sτ l 1-(1- Sgt 1 -Stp )×(1-( Sl + Sgt ) 1 m ) m -( Sgt 1 -Stp )×(1-( Stp ) 1 m ) m 2 (25) k r,g ( Sl ) = (1 -( Sl + Sgt )) τ ′ × (1 -( Sl + Sgt ) 1 m ) 2
(26) This model was implemented by [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] in the iTOUGH2 code (Finsterle et al. ( 2017)), with three scanning curves, first wetting, second drying, and semilog straight lines to approximate the third wetting curve as suggested by [START_REF] Niemi | Preliminary capillary hysteresis simulations in fractured rocks, Yucca Mountain, Nevada[END_REF].

The [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] model uses a clever technique to represent the nonwetting phase entrapment saturation. This model as presented in [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] and in [START_REF] Doughty | User's guide for hysteretic capillary pressure and relative permeability functions in[END_REF] has the disadvantage of approximating (linearizing) the P c (S l ) scanning curves near full water saturation (which is an important state in our studies on clay rock for radioactive waste disposal). In addition, the model is not based on empirical aspects; it only uses the [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF] model applied to a new scaling of the wetting curve after the first reversal point from a main drying path to the first wetting path. [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF], with two parallel sequences of scans departing from the main wetting curve (MWC) and from the first drying curve (FDC). [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] concluded that the hysteresis model with nonwetting phase entrapped suggested by Parker, Lenhard (1987) is based on invalid presumptions, and that the use of the entrapped non-wetting phase formula of [START_REF] Land | Calculation of imbibition relative permeability for two-and threephase flow from rock properties[END_REF] implies non-parallelism between two primary wetting curves departing from the main drying and first drying curve as substantiated in [START_REF] Mualem | A conceptual model of hysteresis[END_REF]. The non-parallelism means that two paths of imbibition that start from the main drying and the first drying curve under an equal change of capillary pressure do not imply the same variation in saturation. [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] affirmed that the model of Parker, Lenhard (1987) presents a superfluous intermediate operation of predicting effective saturation from apparent saturation. Using the scaling technique of [START_REF] Mualem | General scaling rules of the hysteretic water retention function based on Mualem's domain theory[END_REF], [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] proposed a fully explicit equation to predict a desired scanning curve up to the third curve and expressed as follows:

Beriozkin and Mualem 2018 model

• The primary wetting curve from the first turning point (S tp,1 , P c,1 ):

S 1,w (P c ) = S max,w - S max,w -S tp,1 S max,w -S d (P c,1 ) (S max,w -S w (P c )) (27a)
• The secondary drying curve from the second turning point (S tp,2 , P c,2 ) is :

S 2,d (P c ) = S tp,1 + S max,w -S tp,1 S max,w -S w (P c,1 ) S w (P c,2 ) -S w (P c,1 ) S d (P c,2 ) -S d (P c,1 ) (S d (P c )-S d (P c,1 ))

(27b) • And the tertiary wetting curve departing from the third turning point (S tp,3 , P c,3 ) is :

S 3,w (P c ) = S max,w - S max,w -S tp,1 S max,w -S w (P c,1 ) × [(S max,w -S w (P c,2 ))+ S d (P c,2 ) -S d (P c,3 ) S w (P c,2 ) -S w (P c,3 ) S w (P c,2 ) -S w (P c,1 ) S d (P c,2 ) -S d (P c,1 ) (S w (P c,2 ) -S w (P c ))] (27c)
For scanning curves departing from the first drying curve (S ls = 1), a parallel shift is used by adding to each scanning equation (Eq. (27a), Eq. ( 27b), and Eq. ( 27c)) a value of 1 -S max,w .

For a sequence of scanning curves departing from the main wetting curve, the scanning curves are given as follows...

• The primary drying curve from the first turning point (S tp,1 , P c,1 ):

S 1,d (P c ) = S tp,1 S d (P c,1 ) S w (P c ) (28a)
• The secondary wetting curve departing from the second turning point (S tp,2 , P c,2 ) is :

S 2,w (P c ) = S tp,1 -S tp,1 S d (P c,1 ) -S d (P c,2 ) S tp,1 -S w (P c,2 ) S tp,1 -S w (P c ) S d (P c,1 ) (28b)
• And the tertiary drying curve departing from the third turning point (S tp,3 , P c,3 ) is :

S 3,d (P c ) = S tp,1 S d (P c,1 ) × [S d (P c,2 )+ S d (P c,1 ) -S d (P c,2 ) S tp,1 -S w (P c,2 ) S w (P c,3 ) -S w (P c,2 ) S d (P c,3 ) -S d (P c,2 ) (S w (P c ) -S w (P c,2 ))] (28c)
The explicit equations (28) are applied here to normalized water saturations, with S max,d = 1.0 and S l,r = 0 (the maximum drying saturation is unity and the residual liquid saturation is zero). This hysteresis model of [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] gives an explicit formulation of the scanning curve without the need of iterative solution to compute the value in the scanning curves at each grid block and at each iteration. This may be very promising for reducing the computational cost of hysteresis, as noted for instance by [START_REF] Doughty | User's guide for hysteretic capillary pressure and relative permeability functions in[END_REF]. However, this model [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF]) has not been verified and implemented so far.

Recap and conclusions on hysteretic models

In this section, we have described and reviewed hysteresis models with detailed equations and graphics. From the above review of models, we conclude that hysteresis models of [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] and [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] are respectively modifications of the original models of [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF] and [START_REF] Mualem | A conceptual model of hysteresis[END_REF], where the modifications take into account the non-wetting phase entrapment in the hysteresis loops. The non-wetting phase entrapment can be represented in a hysteretic model by different formulas, such as Land's formula [START_REF] Land | Calculation of imbibition relative permeability for two-and threephase flow from rock properties[END_REF]), the linear law of Mualem [START_REF] Mualem | A conceptual model of hysteresis[END_REF]), or Aissaoui's formula [START_REF] Aissaoui | Etude théorique et expérimentale de l'hystérésis des pressions capillaries et des perméabilités relatives en vue du stockage souterrain de gaz[END_REF]). We think that the representation of the non-wetting phase saturation only by the first turning point is based on invalid hypotheses. In fact, this representation does not take into account a higher order of scanning curves and how they can induce release of a fraction of trapped nonwetting phase. For example, taking a primary wetting curve departing from the main drying curve produces a residual gas saturation of S ∆ gr,1 ; after the second drying curve a third wetting curve may produce a residual gas saturation of S ∆ gr,2 ≤ S ∆ gr,1 , which is not taken into account. Even a portion of entrapped non-wetting phase was released due to the second drying path, the model kept the first S ∆ gr,1 entrapment value in order to avoid the pumping errors. The empirical model presented above did not take into account the rate of the process, which is an important factor for fluid transport in porous media as described by [START_REF] Topp | Comparison of water content-pressure head data obtained by equilibrium, steady-state, and unsteady-state methods[END_REF]. [START_REF] Hannes | Revisiting hydraulic hysteresis based on long-term monitoring of hydraulic states in lysimeters[END_REF] proposed three classifications of temporal dynamics of soil water retention: the first class is the non-equilibrium process where fast changes may occur in the water retention; the second class is a quasi-equilibrium process, where relatively slow changes in the water retention are pursued by the fluid redistribution; the third class is the temporal dynamics of soil structures, where the variation of pore structure may occur due to a chemical or mechanical process. Using laboratory data of the soil water retention at different flow processes, they showed that the conceptual models of Parker, Lenhard (1987), [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF], and [START_REF] Poulovassilis | A note on calculating hysteretic behavior[END_REF], failed to reproduce the hysteresis of a non-equilibrium process, especially when the rate of changes in water saturation is high.

In case the capillary entry pressure (P c,e ) is considered in the constitutive relations S l -P c -k r,l -k r,g of the Van Genuchten-Mualem model (Eqs. ( 24), (20), and ( 21)), the representation of hysteresis in the iTOUGH2 code depends on the second primary variable chosen for solving the compositional two-phase flow equations [START_REF] Amri | Parametric Sensitivity to Capillary Entry Pressure in Two-Phase Water-Gas Flow Models: Deep Geologic Disposal of Radioactive Waste // Transport in Porous Media[END_REF]). This second primary variable is either gas saturation or capillary pressure. When gas saturation is used as a primary variable, the release of trapped non-wetting phase in iTOUGH2 is activated by assuming numerically a high saturation at a specified time step. Otherwise, when the capillary pressure is used as a primary variable, the release of trapped non-wetting phase is straightforward (release when P c > P c,e ). However, the present study does not take into account P c,e , and its influence on non-wetting phase trapping is beyond the scope of this paper.

Comparative study of water retention hysteresis models

In the previous section we have analysed and reviewed different hysteresis models one by one; in this section we proceed with inter-comparisons of these hysteresis models in terms of hysteretic water retention curves. Intercomparisons introducing measured hysteretic relative permeability to liquid and gas were not possible because these permeability data were missing for the five materials studied (only saturation-capillary pressure data were available).

Several studies have focused on comparisons of hysteresis models [START_REF] Pham | A study of hysteresis models for soilwater characteristic curves[END_REF], [START_REF] Zhang | A review and statistical study of existing hysteresis models for cementitious materials[END_REF]Terleev et al. (2018)). These studies have showed that the conceptual hysteresis models, which are based on the main wetting and drying curves (and not on the intermediate scanning curves) for their calculations, are simple and more accurate for predicting the scanning curves than the empirical models. Compared to these conceptual models, the empirical models require more difficult calibrations of additional parameters, using more detailed experimental data on intermediate turning points.

On the other hand, the review studies cited just above did not cover all conceptual hysteresis models, and they did not include more recent concepts like non-wetting phase entrapment phenomena. These were taken into account in the more recent studies of [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] and [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF]).

In this section, we test the performance of the five hysteresis models discussed in Section 2 in reproducing the hysteresis data of different porous materials (Sands, Cement, Bentonite) collected from the literature.

Data preparation and comparison criteria

The characteristics of hysteresis data of five porous materials collected from the literature are presented in Table (1). The experimental data of water content θ(P c ) are normalized to S l (P c ) via equation ( 29):

S l = θ -θ r θ s -θ r (29)
The experimental data of the water retention curve measured for incremental relative humidity (RH), i.e. S l (RH), are transformed into S l (P c ) by using Kelvin's equation ( 30):

P c = ρ l RT M l ln(RH) (30)
where R = 8.314 J × K -1 J × mol -1 is the universal gas constant, T (K) is the absolute temperature and M l is the molar mass of liquid water molecule (kg × mol -1 ). First, the Van Genuchten -Mualem model [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1[END_REF]) is fitted to experimental main wetting and drying curves. We make the choice of S ls,d = S ls,w in the fitting procedure of the main curves by equation ( 7) for the hysteresis models [START_REF] Mualem | A conceptual model of hysteresis[END_REF][START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF][START_REF] Haverkamp | Soil water hysteresis prediction model based on theory and geometric scaling // Environmental mechanics: Water, mass and energy transfer in the biosphere[END_REF]. Secondly, the hysteresis models are used to predict the scanning curves based on an iterative solution for scanning curves equations of each model. The goodness of fit of measured scanning points is evaluated using the Mean Absolute Error (M AE) defined by equation ( 31) :

M AE = 1 n n i=0 |S i -Si | (31)
where S i , Si are the measured and predicted values of water saturation, and n is the number of measured data points. We define M AE(N, p -k) as the mean absolute error of the kth measurement curve of the Nth scanning curve at the path p ∈ {w, d}.

Results and discussions

Figures (8,9,10,11,12) present the scanning curves predicted by the five hysteresis models, showing the M AE value for each scanning curve above each plot. The global M AE for each material and for each hysteresis model is presented in Table (2) in order to analyze the different mismatches between the models and data.

From the comparison of the values of M AE, we can conclude that all hysteresis models represent well the tendencies of scanning loops, except the [START_REF] Haverkamp | Soil water hysteresis prediction model based on theory and geometric scaling // Environmental mechanics: Water, mass and energy transfer in the biosphere[END_REF] 2), looking only at the first 2 digits; but more importantly, the Beriozkin and Mualem 2018 model gives an explicit formulation of the saturation of scanning curves, without the need for iterative solution of the equations representing the scanning curves, which is an advantage over the other models [START_REF] Mualem | A conceptual model of hysteresis[END_REF][START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF][START_REF] Haverkamp | Soil water hysteresis prediction model based on theory and geometric scaling // Environmental mechanics: Water, mass and energy transfer in the biosphere[END_REF][START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF]. Modeling two-phase flow with hysteresis To study the effect of hysteresis on the behavior of two-phase flow systems, we have tested a constitutive model with and without hysteresis on the intermediate-scale laboratory experiments conducted by [START_REF] Trevisan | Investigation of mechanisms of supercritical CO2 trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions[END_REF].

The experiments involve the injection of surrogate fluids in order to mimic the behavior of the supercritical scCO2-brine system at the reservoir scale. Despite some differences in the input parameters between our modeling and theirs, we will compare our simulation results with some experimental data, and we will then discuss the effect of hysteresis on non-wetting plume migration.

4.1 The laboratory experiment of [START_REF] Trevisan | Investigation of mechanisms of supercritical CO2 trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions[END_REF] The sand aquifer was initially saturated with water, the wetting fluid. The experiment involved injection of a non-wetting fluid (Soltrol-220) into a rectangular two-dimensional slab of Sand 40/50 (the "aquifer"), with internal dimensions (L×W×H) equal to (91.4× 5.6× 61 cm 3 ) with an inclination angle of 2 • (Figure 13). The right boundary of the aquifer was at equilibrium with a hydrostatic pressure. The well injection was placed on the left boundary of the aquifer. It consisted of a vertical PVC pipe with an internal diameter of 2.3 cm and a vertical length of 11 cm. For more details on this laboratory experiment, see [START_REF] Trevisan | Investigation of mechanisms of supercritical CO2 trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions[END_REF] and [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF].

In our simulation of the experiment, we used a uniform grid with mesh size of 0.5 cm × 0.5 cm to represent the 2D domain shown in Figure ( 13). Our simulation consisted of two stages: a first stage of Soltrol-220 injection at a rate of 0.7ml/min during 5.5 h, followed by a second stage of spontaneous redistribution of the injected fluid during 15 more days.

Hydraulic properties and hysteretic models used

In the present work, the measured data of Sand 40/50 of water retention curves are fitted using the Van Genuchten model as shown in Figure ( 14). We tested two scenarios in the fitting procedure. The first one consists in taking into account the non-wetting phase entrapment phenomenon by considering S w,max as a free parameter in the fitting procedure. The second one consists in neglecting the non-wetting phase entrapment saturation by constraining S w,max to 1.0. The fitted parameters with the two scenarios are given in Table (3).

The hysteresis models of [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] and [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] are used here to describe the mechanism of non-wetting phase entrapment saturation (Figure ( 14)). The [START_REF] Mualem | A conceptual model of hysteresis[END_REF] model is used to describe the hysteresis phenomenon without non-wetting phase entrapment. The intrinsic permeability is k = 6.43×10 -11 [m 2 ] and the porosity is ϕ = 0.35. 2007). (c) [START_REF] Mualem | A conceptual model of hysteresis[END_REF].

Note that there are no data points for the scanning curves, and they are entirely predicted by the afore-mentioned models. The main difference between the models of [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] and [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] is the amount of non-wetting phase entrapment saturation at the end of the first wetting scanning curve. This difference is mainly related to the representation of the non-wetting phase entrapment saturation. For example, the [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] Figure 17 presents a comparison between simulation results without a hysteretic model and with the hysteretic [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] model. The amount of non-wetting phase saturation is clearly different between the two models. The hysteresis model simulates a higher non-wetting phase at the vicinity of the well injection. The imbibition process near the injection well affects the permeability of the non-wetting phase and leads to a higher non-wetting phase saturation than the non-hysteretic model which consists of much higher movement of the injected fluid during time. By the end of the simulation (at time t = 77 hours) a significant difference in the non-wetting phase saturation between the two models is observed. This difference is due to the non-wetting phase entrapment which leads to small permeability in the imbibition region. The same observation was reported by Zhao et al. (2014) where the re-imbibition curves corresponding to the first wetting curve of relative permeability where under the main drainage curve. 3). By the end of the post-injection period, the [START_REF] Mualem | A conceptual model of hysteresis[END_REF] and [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] models simulate the same evolution of the non-wetting fluid. The difference between the two models becomes more significant as time increases (t > 6.94 h). In fact, according to the [START_REF] Mualem | A conceptual model of hysteresis[END_REF] model a small non-wetting phase remains in the system compared to the Beriozkin, Mualem (2018) model, especially near the region where the imbibition process occurred. This may be due to the higher permeability of the hysteretic model without Soltrol-220 entrapment, and also, to the easier migration of the non-wetting phase in the absence of entrapment. In the light of these simulation results, hysteresis may play a significant role in the migration of the non-wetting phase. This result is coherent with the shapes of scanning curves compared to the main drainage (Figure 16.b). In fact, the scanning curve presents a lesser capacity of wetting phase retention compared to the main drainage. This increases the retention of the non-wetting fluid. The non-wetting phase entrapment phenomenon accentuates the high retention of the non-wetting phase, leading to higher quantities of non-wetting fluid in the system.

Discussions on simulated vs. experimental saturations with and without hysteresis

In this subsection, the simulation results are discussed and compared with saturation measurements from the laboratory experiment of [START_REF] Trevisan | Investigation of mechanisms of supercritical CO2 trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions[END_REF]. Due to uncertainties in the input parameters, the present simulation results may not be in agreement with the experimental observations, but nevertheless, the impact of the hysteresis model is clearly observed in our simulations, as in those of [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF]. Qualitatively, the observed differences between [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF] and our simulated results include the following:

• The obtained spatial distribution of the non-wetting fluid S nw (x, y) indicates that the injected non-wetting fluid does not go upwards in our simulation, whereas it does in the simulation of [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF]; • We have indications that, probably for the same reason, the vertical position z c (t) of the center of mass of the non-wetting plume does not go upwards in our simulation, whereas it does in the simulation of [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF].

These differences between our simulations and those of [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF] may be due to the following facts:

• The inclination of the flow domain is possibly different;

• The hydraulic properties for our simulation are possibly different from those of [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF]; • The densities of the fluids are not exactly the same in our simulations and those of [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF].

Figures ( 19) and ( 20) compare the picture taken during the end of the postinjection period and the simulation results. The hysteresis model presents the observed migration plume better than the non-hysteretic model. Uncertainties in the input parameters may affect the simulation results but the trends are clear enough to confirm the accuracy of the hysteretic model to represent the migration of non-wetting phase during a post-injection period.

In order to quantify the amount of the difference between hysteretic and non-hysteretic models, [START_REF] Trevisan | Investigation of mechanisms of supercritical CO2 trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions[END_REF] used experimental saturation data measured by the X-Ray attenuation method to calculate spatial moments of the non-wetting plume, based on the non-wetting saturation distribution S nw (x, z). In particular, the center of mass coordinates x c (t) and z c (t) were analyzed. However, because experimental data are not entirely available from the experimental study of [START_REF] Trevisan | Investigation of mechanisms of supercritical CO2 trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions[END_REF] and the simulation study of [START_REF] Cihan | Investigation of representing hysteresis in macroscopic models of two-phase flow in porous media using intermediate scale experimental data[END_REF], this type of spatial moments analysis was not possible in our simulation study due, seemingly, to lack of information on some physical parameters of the experiment. From the 2D simulations results, we can conclude that the hysteresis models are more accurate to predict the non-wetting plume migration in the post-injection period. The comparison of hysteresis models also showed that neglecting the non-wetting phase entrapment may affect the simulation results at large times.

Conclusions

In this study, we have presented a detailed review, and a quantitative assessment, of several of the most important and representative conceptual models for hysteretic two-phase flow in porous media. Model equations and hypotheses were analyzed in detail, with comments on the required solution steps in the presence of turning points due to wetting-drying cycles. Globally, our survey indicates that the hysteretic models of two-phase flow agree fairly well with the experimental data, compared to the non-hysteretic models.

Given the goodness of fit of the [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] model to experimental data, and its explicit calculation of saturation scanning curves, and also, its taking into account of the non-wetting phase (gas) entrapment, it can be considered as a good candidate for implementing hysteresis in 2-phase flow simulations. The models of [START_REF] Mualem | A conceptual model of hysteresis[END_REF] and [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF] also give good results on experimental data but without representing the gas entrapment phenomenon. Finally, the model of [START_REF] Haverkamp | Soil water hysteresis prediction model based on theory and geometric scaling // Environmental mechanics: Water, mass and energy transfer in the biosphere[END_REF] gives the least accurate results among the models that we have selected for this review with respect to the experimental data that we have analyzed and simulated.

Our simulation results highlight the ability of hysteretic models to reproduce the post-injection non-wetting plume distribution in the intermediatescale laboratory experiment, compared to non-hysteretic simulations. Neglecting the non-wetting phase entrapment may also affect the simulation results. The modeling of this experiment, where a non-wetting fluid is injected in a reservoir initially saturated with a wetting fluid, confirms the necessity of considering hysteresis in modeling such multiphase flow systems. The [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] model is found to be the most adequate to model hysteresis due to its advantages of reproducing well the experimental data, and relative simplicity of implementation. Furthermore, it is expected that this model is the best one in terms of CPU times, because it solves explicit equations for calculating hysteretic scanning curves.

Further developments should take into account, in combination with hysteretic effects, the phenomenon of capillary trapping of the non-wetting phase (e.g., hydrogen gas), and also, the existence of a capillary entry pressure which can induce a gas entry pressure threshold: see [START_REF] Amri | Parametric Sensitivity to Capillary Entry Pressure in Two-Phase Water-Gas Flow Models: Deep Geologic Disposal of Radioactive Waste // Transport in Porous Media[END_REF] for a recent sensitivity study of water-gas flows with and without entry pressure.

Figure 1 :

 1 Figure 1: Schematic representation of wetting/drying processes based on the independent domain theory diagram (Mualem (1974)): (a) main wetting; (b) main drying.

Figure 2 :

 2 Figure 2: Schematic representation of wetting/drying processes based on the independent domain theory diagram (Mualem (1974)): (a) first wetting scanning curve; (b) first drying scanning curve.

Figure 3 :

 3 Figure 3: Schematic illustration of soil water hysteresis with the independent domain theory of Mualem (1974): (a) continuous sequence of scanning curves departing from the main drying curve; (b) continuous sequence of scanning curves departing from the main wetting curve.

Figure 4 :

 4 Figure 4: Schematic illustration of soil water hysteresis with the dependent domain theory of Mualem (1984): (a) continuous sequence of scanning curves departing from the main drying curve; (b) continuous sequence of scanning curves departing from the main wetting curve. Idem Van Genuchten parameters of main curves in figure (3).

Figure 5 :

 5 Figure 5: Schematic illustration of soil water hysteresis with the geometrical scaling model of Haverkamp et al. (2002): continuous sequence of scanning curves departing from the main drying curve. Idem VG-parameters in figures (3) and (4).

Figure 6 :

 6 Figure 6: Schematic illustration of water-gas hysteresis using the Doughty (2007) model. The parameters of Van Genuchten model are the same as in the figure (3)

Figure 7 :

 7 Figure 7: Schematic representation of the hysteresis model of Beriozkin, Mualem (2018), with two parallel sequences of scans departing from the main wetting curve (MWC) and from the first drying curve (FDC).

Figure 8 :Figure 9 :

 89 Figure8: Fitted main wetting and main drying curves to the Dune Sand (material N • 1) data fromTerleev et al. (2018), and three predicted primary wetting scanning curves departing from the main drying branch, according to the five hysteresis models (one model per subfigure, respectively, from (a) to (e)).

Figure 10 :Figure 11 :

 1011 Figure 10: Fitted main wetting and main drying curves to the Compacted Sand (material N • 3) data from Huang et al. (2005), and predicted scanning curves according to the five hysteresis models (one model per subfigure, from (a) to (e)).

Figure 12 :

 12 Figure12: Fitted main wetting and main drying curves to the Hostun Sand (material N • 5) data from[START_REF] Lins | Physical Modeling of SWCC for Granular Materials // Theoretical and Numerical Unsaturated Soil Mechanics[END_REF], and predicted first wetting curve and second drying scanning curve according to the five hysteresis models (one model per subfigure, from (a) to (e))..

Figure 13 :

 13 Figure 13: Schematic representation of the 2D domain and boundary conditions corresponding to the laboratory experiment of Trevisan et al. (2014) as simulated in the present paper.

.

Figure 14 :

 14 Figure 14: Illustration of the hysteresis effect on retention curve properties of Sand 40/50. Scanning loops are calculated with different models. (a) Beriozkin, Mualem (2018). (b) Doughty (2007). (c)[START_REF] Mualem | A conceptual model of hysteresis[END_REF].

  model uses a linear model to represent the relation S gr -S tp which gives a non-wetting phase entrapment saturation of 0.105 when the turning point saturation S tp = 0.41. However, the Doughty (2007) model uses the formula 4.3.1 Effect of different hysteresis models on the distribution of injected non-wetting fluid during post-injection As presented in section 3, the Beriozkin, Mualem (2018) hysteretic model presents a good agreement with experimental data of different porous materials types, and we focus mainly on this hysteretic model for simulations of the twophase flow laboratory experiment. We also compare with simulations obtained with another hysteresis model (Doughty (2007)), and with simulations without hysteresis. We choose here to implement these two-phase hysteretic models, as well as the non-hysteretic model, with the iTOUGH2 code (Finsterle et al. (2017)). Figure (16) shows the 2D distribution of the non-wetting phase (Figure 16.a) and the curve number or index number (Figure 16.b) at different postinjection times. At the end of the injection, a new index curve equal to 2 starts to appear in the vicinity of the injection well, corresponding to the first imbibition path in the P c -S l relationships. At the same time the region of drainage (icurve = 1) starts migrating towards the other boundary.

Figure 16 :

 16 Figure 16: (a) Simulated 2D distribution of the non-wetting fluid saturation at different times by Beriozkin, Mualem (2018) hysteresis model. (b) Flow reversal index numbers where icurve = 0 designates fully saturation state (S l =1), icurve = 1 designates main drying curve, and icurve = 2 designates the first wetting curve.

  From the results shown in figures (16.a) and (17.b), no big differences are observed between Beriozkin, Mualem (2018) and Doughty (2007) models.

Figure 17 :

 17 Figure 17: Simulated 2D distribution of the non-wetting fluid (Soltrol-220) saturation at different times (t = 5.5, 6.94, 22.4, 77 hours). (a) Results from the model without hysteresis. (b) Results with the hysteretic model of Doughty (2007).

Figure 18 :

 18 Figure 18: Simulated 2D distribution of the non-wetting fluid (Soltrol-220) saturation at different times (t = 5.5, 6.94, 22.4, 77 hours). (a) Results with the hysteretic model of Mualem (1974). (b) Results with the hysteretic model of Doughty (2007).

Figure 19 :

 19 Figure 19: Picture taken during the experiments at the end of the recovery period t = 15 days from Cihan et al. (2017).

Figure 20 :

 20 Figure 20: Simulated 2D distribution of non-wetting fluid saturation at the end of the post-injection period t = 15 days. (a) Result with the hysteretic model of Beriozkin, Mualem (2018). (b) Result with the non-hysteretic model.

Table 1 :

 1 Hydraulic properties of different porous materials studied.

	Material References	Type of soil	Scanning
				curves
	1	Terleev et al. (2018)	Dune sand	1st wetting
	2	Baroghel-Bouny (2007)	Paste(CEM I)	1st wetting
				2nd drying
	3	Huang et al. (2005)	Compacted Sand 1st wetting
				2nd drying
				3rd wetting
	4	Qiao et al. (2021)	Bentonite MX80	1st wetting
	5	Lins et al. (2007)	Hostun sand	1st wetting
				2nd drying
				3rd wetting

Table 2 :

 2 model which fails to predict the main wetting branches M AE results for the scanning curves fit of each material's experimental data by the five hysteresis models. The column N • refers to the number of the porous material as presented in Table (1). for all materials (except material N • 4). The Mualem 1974 and Mualem 1984 models share the same first drying curve, but the Mualem 1974 model exhibits better results than Mualem 1984 for higher order scanning curves (e.g., materials N • 2, 3). The Doughty 2007 model is a modification of the Mualem 1984 model by introducing a new wetting curve, determined after a reversal point from main drying to first wetting. The Doughty 2007 model presents a high error compared to the Mualem 1984 model, and therefore the introduction of nonwetting phase entrapment saturation enhances the Mualem 1974 model, but not the Mualem 1984 model. Finally, the Beriozkin and Mualem 2018 model is slightly better, particularly compared to Haverkamp et al. 2002 and Doughty 2007 based on mean M AE values in Table (

	N •	Mualem 1974 Mualem 1984 Haverkamp 2002 Doughty 2007 B&M 2018
	1	0.05737	0.05738	0.04148	0.05656	0.06040
	2	0.02511	0.03221	0.24811	0.05723	0.02397
	3	0.02889	0.05733	0.23982	0.08832	0.02799
	4	0.04069	0.04066	0.04322	0.04062	0.04062
	5	0.02472	0.01998	0.15304	0.05887	0.02156
	Mean 0.03535	0.04151	0.14513	0.06032	0.03491

Table 3 :

 3 Adjusted hydraulic parameters for the Sand 40/50. Scenario 1 is for S w,max as a free parameter. Scenario 2 is for fixed S w,max = 1.0.

	Scenario	nw	αw [P a -1 ]	Sw,max	Sw,r	n d	α d [P a -1 ]	S d,max	S d,r
	1	2.95 3.89×10 -3	0.8431	0.0073 3.91 1.76×10 -3	1.00	0.0097
	2	2.56 5.26×10 -3	1.0	0	3.91 1.76×10 -3	1.00	0.0097

(a) First wetting curve departing from main drying curve (b) First drying curve departing from main wetting curve

50, 9. 7067-7081.

Acknowledgments

The authors would like to thank the anonymous reviewers for their insightful comments and editorial suggestions which enhanced the content of the paper. The first two authors also gratefully acknowledge the funding provided by EURAD, the European Joint Programme on Radioactive Waste Management for DONUT (Grant Agreement N°847593).

of [START_REF] Land | Calculation of imbibition relative permeability for two-and threephase flow from rock properties[END_REF], which gives a non-wetting phase entrapment saturation of 0.147 for the same turning point S tp = 0.41. This difference between the two formulae concerning the S gr -S tp relationships is presented in Figure ( 15) in comparison with experimental data of Sand 40/50 from [START_REF] Trevisan | Investigation of mechanisms of supercritical CO2 trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions[END_REF]. The fit of experimental data of S gr -S tp reveals that the Aissaoui's formulation [START_REF] Aissaoui | Etude théorique et expérimentale de l'hystérésis des pressions capillaries et des perméabilités relatives en vue du stockage souterrain de gaz[END_REF]) is in good agreement with measured data. This observation is also in agreement with [START_REF] Suzanne | Experimental relationships between residual gas saturation and initial gas saturation in heterogeneous sandstone reservoirs[END_REF] and [START_REF] Pentland | Measurement of Nonwetting-Phase Trapping in Sandpacks[END_REF].

Simulation results and discussion

In this section, we present a comparison between the simulation results obtained respectively with (i) a non-hysteretic model, (ii) the hysteretic Beriozkin, Mualem (2018) model, (iii) the hysteretic Doughty (2007) model, and (iv) the hysteretic [START_REF] Mualem | A conceptual model of hysteresis[END_REF] model. First, the effect of hysteresis is discussed by using the [START_REF] Beriozkin | Comparative analysis of the apparent saturation hysteresis approach and the domain theory of hysteresis in respect of prediction of scanning curves and air entrapment[END_REF] model. Secondly, the effect of non-wetting phase entrapment is discussed by comparing [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of nonhysteretic and hysteretic characteristic curves // Energy Conversion and Management[END_REF] and [START_REF] Mualem | A conceptual model of hysteresis[END_REF] models. Finally, a comparison between the simulation results and the measured saturations from the laboratory experiment is presented.