Introduction

The modelisation of the impregnation of porous solid skeleton by liquid is a research topic with a wide scope of applications. In industry, this process is involved in a large variety of contexts, from corrosion of refractory materials [START_REF] Blond | Effect of Slag Impregnation on Thermal Degradations in Refractories[END_REF] to resin flow in Resin Transfer Molding (RTM), composite process [START_REF] Govignon | Simulation of the reinforcement compaction and resin flow during the complete resin infusion process[END_REF] via leaching of concrete in civil engineering [START_REF] Mainguy | Modelling of leaching in pure cement paste and mortar[END_REF] or thermohydromechanical effect of nuclear waste storage [START_REF] Mainguy | Modelling of leaching in pure cement paste and mortar[END_REF]. Often, the capillary suction which drives the impregnation process is coupled with chemical reactions between the liquid and the solid skeleton, like refractory corrosion [START_REF] Blond | Effect of Slag Impregnation on Thermal Degradations in Refractories[END_REF] and concrete leaching [START_REF] Mainguy | Modelling of leaching in pure cement paste and mortar[END_REF]. Furthermore, in the resin for RTM process, there is a reaction (and/or a phase change) that take places only in the liquid.

The numerical modelling of such reactive impregnation requires a multiphysics model taking into account the material properties such as porosity and pore size distribution of the porous medium, density and viscosity of the impregnation liquid, the transport properties such as permeability and capillary suction, as well as the type of chemical reactions such as phase changes, dissolution, precipitation, layers growth … All these physics must be strongly coupled as the geometric and chemical properties may change because of the reactions. In addition, mechanical properties may also evolve and new phases may induce interactions with the solid skeleton.

To perform such strongly thermochemical -thermomechanical coupled simulations, two main ways seem to be relevant: (i) coupling dedicated codes by Mesh-based Parallel Code Coupling Interface (MPCCI)-type approach which are more developed in combustion and energy field [START_REF] Refloch | CEDRE Software[END_REF] and (ii) implementing the most significant part of thermochemistry in codes dedicated to simulate the thermomechanical part of the reaction [START_REF] Merzouki | Numerical study of the effects of refractory lining geometries on the swelling induced by oxidation[END_REF][START_REF] Merzouki | Modelling of the swelling induced by oxidation in SiC-based refractory castables[END_REF]. Whatever the chosen approach is, the computational cost is so high that it often requires important computing facilities.

The aim of this thesis is to propose a new method, which does not use classical modelling by partial differential equations and the associated numerical methods, allowing to drastically reduce the computational cost. For the sake of clarity, this thesis presents the basic lines of our approach applied to the (quasi) one-dimensional non-reactive unsaturated impregnation.

The "traditional" way of describing impregnation is through a non-linear Partial Differential Equation (PDE) which corresponds to the mass balance combined to Darcy's law, known as the "Richards' equation" (here for isotropic case) [START_REF] Coussy | Mechanics and Physics of Porous Solids: Coussy/Mechanics and Physics of Porous Solids[END_REF]

: = ( (- + )) (1)
where is the local volumic saturation (dimensionless), defined by the ratio of local liquid volume to local porous volume, is the intrinsic permeability ( ), is the dynamic 10/103 viscosity (Poi); is the porosity, defined by the ratio of the volume of pore to the total volume (dimensionless); is the capillary pressure, i.e. the sorption pressure driving the capillary rising (Pa), and is a function correcting the permeability according to saturation (dimensionless), is the mass density of the wetting liquid ( • -), is acceleration of gravity ( • -), and is the transient height of the liquid (m). Often, this equation is rewritten as:

ϕ = ( + ) (2)
This formulation allows to solve the PDE using the saturation as principal unknown. The function and the evolution of capillary pressure with saturation and more precisely its derivative become the parameters. Because of the lack of experimental data, the function is often assumed to be equal to unity, a cubic function or more complicated functions of saturation [START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1[END_REF]. Conversely, the curve is experimentally obtained for the steady state (i.e., equilibrium under only gravity field) and is denoted as the Capillary Pressure Curve (CPC). Capillary pressure is the equilibrium pressure difference between the non-wet and wet phases [START_REF] Coussy | Mechanics and Physics of Porous Solids: Coussy/Mechanics and Physics of Porous Solids[END_REF]. This curve is usually directly recorded as a tabular data for the simulation or is reproduced by appropriate well-known fitting functions (e.g., Brook's model or van Genuchten's model [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF][START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1[END_REF].

It is difficult to obtain analytic solutions of Richards' equation (non-linear) except for special cases [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF][START_REF] Zhang | Finite Analytic Method for Solving the Unsaturated Flow Equation[END_REF]. In order to obtain optimal numerical solutions of Richards' equation, the Capillary Pressure Curve (CPC) needs first to be determined. At a local scale, it is possible to assimilate the local pore space to a vertical capillary. Capillary pressure is thus created by the curvature of the meniscus and the surface tension of liquid, causing to the liquid to rise. As the liquid rises, capillary pressure acts against gravity applied to the liquid column [START_REF] Extrand | Forces, pressures and energies associated with liquid rising in nonuniform capillary tubes[END_REF]. Main methods to determine the CPC are:

(i)
The phenomenological method established and developed for the prediction of the CPC using closed-form analytical expressions [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF][START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1[END_REF]); (ii)

The Boltzmann transformation, under the assumption that the CPC is the mirror image of the extended Boltzmann transformation function [START_REF] Prevedello | New Analytic Solution of Boltzmann Transform for Horizontal Water Infiltration into Sand[END_REF][START_REF] Prevedello | The Similarity Hypothesis and New Analytical Support on the Estimation of Horizontal Infiltration into Sand[END_REF]); (iii) The morphological pore network model implying that local capillary pressure is a function of the radius of a hemisphere (which represents the air-water interface) by employing the Young-Laplace equation [START_REF] Dodds | Capillary Pressure Curves of Sphere Packings: Correlation of Experimental Results and Comparison with Predictions from a Network Model of 99/103 Pore Space[END_REF][START_REF] Sweijen | Capillary pressure-saturation relationships for porous granular materials: Pore morphology method vs. pore unit assembly method[END_REF]; (iv)

The fractal pore-space model implying that the CPC is a function of pore-solid interface fractal dimension determined by the percolation theory (Ghanbarian-Alavijeh and [START_REF] Ghanbarian-Alavijeh | Unsaturated hydraulic conductivity in porous media: Percolation theory[END_REF][START_REF] Hunt | Unsaturated hydraulic conductivity modeling for porous media with two fractal regimes[END_REF].

11/103

Assuming the CPC is known, several numerical methods have been used to solve Richards' equation [START_REF] Bouchemella | Numerical simulation of water flow in unsaturated soils: comparative study of different forms of Richards's equation[END_REF][START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF][START_REF] Szymkiewicz | Modelling Water Flow in Unsaturated Porous Media, GeoPlanet: Earth and Planetary Sciences[END_REF][START_REF] Wu | A finite-element algorithm for modeling variably saturated flows[END_REF][START_REF] Zhang | Finite analytic method based on mixed-form Richards' equation for simulating water flow in vadose zone[END_REF][START_REF] Zhang | Finite Analytic Method for Solving the Unsaturated Flow Equation[END_REF]. For all of these methods, small timestep and fine space discretization are required to ensure both convergence and accuracy [START_REF] Szymkiewicz | Modelling Water Flow in Unsaturated Porous Media, GeoPlanet: Earth and Planetary Sciences[END_REF], which limit speed of convergence due to a problem of mass-balance errors [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF][START_REF] Wu | A finite-element algorithm for modeling variably saturated flows[END_REF].

In order to minimize this problem, the modified Picard's iteration, also known as the quasi-Newton method, used with both Finite Element Method (F.E.M) and Finite Different Method (F.D.M) has been proposed by M. A. [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF]. The problem related to FEM is then settled by using the modified head method [START_REF] Wu | A finite-element algorithm for modeling variably saturated flows[END_REF] or by incorporating the local analytic solution in the computational scheme [START_REF] Zhang | Finite analytic method based on mixed-form Richards' equation for simulating water flow in vadose zone[END_REF]. The speed of the convergence is bigger, while the existence of the diagonal matrix for the one-dimension case for the FDM [START_REF] Bouchemella | Numerical simulation of water flow in unsaturated soils: comparative study of different forms of Richards's equation[END_REF] and the Finite Volume Method (FVM) [START_REF] Lai | A mass-conservative finite volume predictor-corrector solution of the 1D Richards' equation[END_REF] guarantee the simulation accuracy [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF]. However, small time-increment still leads to high computational cost and, furthermore, a very small time-increment can lead to spurious oscillations [START_REF] Vermeer | An accuracy condition for consolidation by finite elements[END_REF] that impact the accuracy of the results.

This synthesis underlines that it is not trivial to perform accurate and low time-consuming resolution of the Richards' equation. This observation is at the origin of the idea to developing a new numerical method free from the described numerical drawbacks.

The objective of this thesis is to develop the proposed approach and the associated new numerical methods to predict the time evolution of the capillary pressure profiles without spurious oscillations and with reduced computational cost.

The thesis is organized in the following way. In Chapter I, theoretical backgrounds of (nonreactive) impregnation in porous media are first introduced to understand the mechanisms of the fluid transport in porous media in section I.1. Also, in section I.2, the numerical analyses of the phenomenon are reviewed; and the associated drawbacks are presented as goals. In Chapter II, an introduction to the traditional percolation model, to the gradient percolation model, and finally to the Self-organized Gradient Percolation (SGP) model is proposed. This theory, which is a random probability model of site percolation where the probability for each site varies along fixed direction [START_REF] Nolin | A Modified van Genuchten-Mualem Model of Hydraulic Conductivity in Korean Residual Soils[END_REF], has already been applied to the simulation of diffusion [START_REF] Sapoval | The fractal nature of a diffusion front and the relation to percolation[END_REF] and growth of cities [START_REF] Nguyen | Simulations of the growth of cities[END_REF]. Mathematical background, the performance of the SGP algorithm, and characterizations of input data of this algorithm are also presented. Then, in sections II.3 to II.3.4, the new numerical algorithm based on the gradient percolation theory is proposed to produce numerical simulation of the impregnation of a porous media without solving the Richards' equation. The initialization of the algorithm is driven by an analytic solution of the homogeneous diffusion equation, which is a convolution between a Probability Density Function (PDF) and a smoothing function [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF]. The evolution of the capillary pressure profile with time is reproduced by the self-evolution of the standard deviation of the PDF. This model is therefore called Self-organized Gradient Percolation model (SGP). The proposed algorithm is evaluated by comparison with numerical results from FEM and from experimental measurements. Finally, results, conclusions and outlooks are drawn in chapters III and IV.

12/103 13/103

Chapter I Bibliography In this chapter, the essential backgrounds of impregnation in porous media are provided. In particular, researches both on theoretical and numerical analyses in the field of the fluid transport in porous media are reviewed.

The first study is the general concepts of impregnation in porous media at the (macroscopic) scale of the skeleton and pore space [START_REF] Chandler | Capillary displacement and percolation in porous media[END_REF] in Section I.1. This leads to the introduction of porous medium definition, capillary pressure and saturation, and capillary rise which includes the significant variables of the finite volume of a porous media such as the porosity and the saturation in Sections I.1.1 to I.1.3.We will then present the formulation of its governing equation in Sections I.1.4 and I.1.5. To be solved, the governing equation demands a list of properties to be measured. Herein, the Capillary Pressure Curve (CPC) is one of the parameters, which plays an important role in the modeling of the phenomenon. We will therefore report several well-known methods of determining the capillary pressure curve.

The second study is the numerical methods (e.g., Finite Element Method "F.E.M", Finite Difference Method "F.D.M", etc.) having often been used to solve this governing equation in Section I.2. As the objective of the thesis, the discussions focused on the numerical difficulties of these associated numerical methods are also pointed out in order to think of a novel numerical method which allows for overcoming the drawbacks.

I.1 General concepts of impregnation in porous media

This section is devoted to present the fundamentals of the impregnation phenomenon in porous media.

I.1.1 Porous medium definition

According to [START_REF] Bear | Transport Phenomena in Porous Media -Basic Equations[END_REF] a porous material of concern in this work is a medium defined as a portion of space partly occupied by a solid phase (or solid matrix, or skeleton) and partly by a void space. The solid phase gives the cohesion to the material and must be thereby continuous and distributed throughout the whole medium domain. There are two types of void space: i) void space made of isolated spaces which are not connected each other; this space is not involved in the impregnation and will be included in the solid matrix. ii) the void space made of only one connected network through which fluids can flow. In this document, porous network, pore space, or even porosity will refer to the connected pore space. The pore space can be occupied by one or more fluid phases, gaseous or liquid. It can be heterogeneous as in most of refractory materials or may show periodic pattern as in composite 14/103 reinforcements at mesoscopic or microscopic scale. In addition, it must be possible to associate a Representative Elementary Volume with a porous medium. This way, the porous material can be defined by macroscopic properties and state variables in a continuum approach. Among all the macroscopic properties, the first one is the porosity. The porosity (dimensionless) of a porous medium is quantified for a representative volume of material by the fraction of the volume of the porous volume to the sum of porous volume and skeleton volume (called bulk volume ). The formulation of the porosity is then given as:

= = + (3)

I.1.2 Capillary pressure and saturation

A porous space can be saturated or unsaturated. When saturated only one fluid phase saturates the space. When unsaturated at least two fluids share the porous space. In this work, we will consider unsaturated conditions with only two fluids where one will be the non-wetting phase, generally a gas, whereas the other one will be the wetting phase such as water, resin, slag…. In Fig. 1, two domains are saturated. At the bottom, all the pores are completely filled by a phase which is supposed to be the wetting phase. The upper part of the porous medium is also saturated by the non-wetting phase. Note that unsaturated state does not exclude that a part of the porous domain be saturated. Likewise, some pores may be saturated although they lie in an unsaturated domain. In this study, we will consider the wetting phase to define the saturation:

• The porous space is saturated when it is exclusively occupied by the wetting phase.

• The porous space is not filled when it is saturated by the non-wetting phase.

15/103

Hence, at a large enough scale corresponding to the representative elementary volume, it is assumed that non-wetting phase and wetting phase entirely occupies together the pores in the unsaturated domain. This leads to:

= + (4)
where and are the volumes occupied by non-wetting and wetting phases, respectively.

It is then possible to define the state of occupation of the pores by the wetting phase in a given representative volume by the so-called saturation: = ;

(5)

And it is obviously possible to define the saturation state from the non-wetting phase view point:

= ; (6)
with

+ = (7)
It is well known that, at microscopic scale, ⁄ phases interface form a meniscus resulting from energy minimisation involving the pressures of the phases, their wettability and the geometry of the pore or throat. Assuming cylindrical pore (Fig. 2), the Laplace equation takes the form:

= - = -cos (8)
where and are the pressures of the non-wetting and wetting phases, respectively, is the cylinder diameter,

-the surface interface energy of the interface between the fluids, and the contact angle between the wetting phase and the wall. The pressure difference is usually called the capillary pressure .

Additionally, Jurin's law defines the liquid height in the capillary at the equilibrium [START_REF] Coussy | Mechanics and Physics of Porous Solids: Coussy/Mechanics and Physics of Porous Solids[END_REF].

ℎ = - ( 9
)
where ℎ is the height of the penetration zone , is the density of the wetting fluid ⁄ , and is the acceleration of gravity -⁄ .

16/103

Lucas Washburn equation determines the penetration rate of the liquid [START_REF] Bell | The Flow of Liquids through Capillary Spaces[END_REF][START_REF] Washburn | Note on a Method of Determining the Distribution of Pore Sizes in a Porous Material[END_REF] and the height with time:

ℎ = √ -• • • ∝ √ (10)
where is the dynamic viscosity of the invading fluid and is the time.

A real porous network is made up of pores connected with throats making up a very complex geometry very different from the simple capillary tube and these equations do not allow for an accurate and realistic simulation of the impregnation in unsaturated conditions. Upscaling now to the porous medium and according to thermodynamic considerations, the capillary pressure = is actually the conjugate state variable of the wetting phase saturation SW [START_REF] Coussy | Mechanics and Physics of Porous Solids: Coussy/Mechanics and Physics of Porous Solids[END_REF]. Capillary pressure is a state function of the saturation decreasing from a maximum value when only the non-wetting phase is present ( = ) to a given pressure when the wetting phase saturates the porous space ( = ). This pressure, called entry pressure , is the pressure drop required for the non-wetting phase to start entering the porosity. Capillary pressure is usually modelled by an empirical function:

< - > = Π • Π (11)

17/103

where < -> is the positive difference between the capillary pressure and the entry pressure, Π is a reference pressure, and Π is a dimensionless function with Π = = (Fig. 3).

The macroscopic function Π results from the interactions between the fluids and the pore walls and can be empirical or derived from physical considerations; it will be presented in detail in section I. 1.5.3. Due to the interactions between the wetting phase, the non-wetting phase, and the porous medium, the wetting phase maybe "trapped" inside some small pores and the minimum (or residual) value of the saturation (denoted by ) is strictly greater than zero. On the other hand, the wetting phase may not fill the pores completely. The saturation is therefore limited by a maximum value (denoted by

). It is usual to use then effective saturation defined as follows [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF]):

= - - ; (12)
Fig. 3: Example of capillary pressure function, after [START_REF] Albers | Main Drying and Wetting Curves ofSoils: On Measurements, Prediction and influence of Wave Propagation[END_REF].

I.1.3 Capillary rise

Before addressing the general impregnation process and its simulation, capillary rise phenomenon is studied in this section to make the required properties appear. In this simple case, a dry porous medium with a cylindrical shape comes into contact with a liquid at its bottom (Fig. 4). The sample has been coated on its cylindrical surface to prevent flow through this surface. The problem can thereby be considered as 1D problem. The gas (air) within the porous 18/103 space can freely leave the porous space through the upper surface. The container is covered by a lid to insulate the liquid and apply a given pressure P0. The air outside the sample is at the atmospheric pressure. The impregnation is supposed non-reactive. The air is the non-wetting phase and the liquid is the wetting one. When equilibrium is reached, it is observed that a part of the sample is fully impregnated with the liquid at the bottom of the sample. An unsaturated zone lies above the saturated zone. The height of the saturated zone is noted h. Non-reactive impregnation modelling is derived from mass balance and momentum balance applied to both phases:

( ) = -( ⃗⃗ ) (13) ( ) = -( ⃗⃗ ) (14)
where and are the density of the air and liquid phase, respectively. ⃗ stands for the flow density vector. The remove of the air is assumed independent of the liquid imbibition and the air pressure is assumed to be equal to the atmospheric pressure. The flows are therefore decoupled and only liquid mass balance is required to study the impregnation. This assumption is called Richards' conditions.

In 1856, [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon: Exposition et Application des Principes a Suivre et des Formulesa Employer dans les Questions de Distribution d[END_REF] first introduced the formulation of the flow density vector (nowadays also known as Darcy's law) [START_REF] Coussy | Mechanics and Physics of Porous Solids: Coussy/Mechanics and Physics of Porous Solids[END_REF]; and since it has been used widely and taken under account the basic equations of the fluid transport in porous medium. Its local expression is:

19/103 ⃗ = - (15)
where denotes the "permeability" tensor of the porous medium with respect to the fluid (m 2 /Pa/s) and is the total pressure of the wetting phase (i.e. the liquid). is the driving force to which body forces can be added. Generally, body forces reduce to the gravity force -⃗ acting in the direction.

In the isotropic case, reduces to scalar k in its diagonal and is still related to the properties of the fluid (i.e. viscosity), and then the Darcy's law can be given by the term of the capillary pressure [START_REF] Coussy | Mechanics and Physics of Porous Solids: Coussy/Mechanics and Physics of Porous Solids[END_REF] as:

⃗ = - ( + ) (16)
Assuming the porosity constant over the whole domain and since the liquid density is constant, the mass balance equation writes:

= ( + ) (17)
At equilibrium, the flow stops and all over the sample height the liquid pressure is determined from:

+ = (18)
• For ℎ, = where ℎ is the height of the saturated zone (Fig. 4) and = , leading to:

= - (19)
• For ℎ , and the capillary pressure depends on the saturation (eq. ( 9)). At = ℎ, = and Π = = while = (Richards' conditions). According to the definition of the capillary function:

< - > = < - - > = (20)
The wetting phase pressure at = ℎ is then:

ℎ = - = - ℎ (21) If = 20/103 ℎ = / (22) Else ℎ = + - (23)
Since ℎ , the inner pressure must verify -. If not, there will not be any saturated zone and unsaturated zone will start at the bottom of the porous medium.

At

ℎ, changing the z origin and combining Eqs. ( 8) and ( 9) with Eq. ( 19) leads to:

- = -ℎ = Π • Π (24)
The saturation profile can be then derived from the reciprocal capillary function:

= = Π - -ℎ Π (25)
In other words, the saturation profile observable in a capillary rising test is an homothety of the capillary pressure curve.

I.1.4 Unsaturated impregnation and list of properties to be measured

Previous section addressed the equilibrium state of 1D impregnation of a porous medium with a wetting phase in Richards' conditions. In general case Eqs. ( 14) and ( 16) still apply:

= ( ( + )) (26)
where k has been expressed in function of the permeability of the porous medium, (m 2 ), and the dynamic viscosity of the wetting phase,

• . One way to solve the PDE ( 26) is to express as a function of and :

= - ( ( - ⃗)) (27)
where is the intrinsic permeability of the material (m 2 or Darcy with 1 D ≈ 10 -12 m 2 ), is the relative permeability (dimensionless) which modifies the permeability to account for the partial disconnection of the wetting phase.

To be solved, Eq. ( 27) requires five data , , , , and . Porosity and dynamic

21/103

viscosity can be measured experimentally [START_REF] De Bilbao | Intrinsic permeability of refractories from gas permeability measurements: Comparison of results[END_REF]. This leads to that the rest macroscopic transport properties become the main unknowns to be identified as the following:

• intrinsic permeability ; • relative permeability ; • derivative of capillary pressure regards to saturation .

In the next sections, we will briefly present these transport properties, as well as the determination of their value.

I.1.5 Determination of the properties involved

in the unsaturated impregnation process

I.1.5.1 Intrinsic permeability

Intrinsic permeability is a material property which consists of a length scaling geometry of the porous network and the porosity. That is why the determination of its value strongly depends on the geometry of the porous network with the details following:

• For simple geometries of the porous network (e.g. sand, concrete, stone and clays), the order of magnitude of intrinsic permeability can be directly deduced from the simplest material properties. In particular, it is derived from the length and the function of the porosity characterizing the geometry and assumed in accordance with the Kozeny-Carman formula; • For more complex geometries, often, the intrinsic permeability is measured experimentally.

As a consequence, experimental methods have been often used to determine the intrinsic permeability. Indeed, there are two well-known methods of experimental measurements: the steady-state and unsteady state permeability measurement methods. For the refractories and even for materials of very low permeability, the steady-state permeability measurement method is more suitable both for gas and liquid. Although gas permeability differs from liquid permeability as shown below, gas permeability measurement is easier to perform than liquid one thanks to the gas compressibility. In steady-state conditions, gas and liquid permeability are derived from integrated form of Darcy's law provided laminar flow [START_REF] Klinkenberg | The Permeanility of Porous Media to Liquids and Gases[END_REF]:

Liquid Permeability Gas Permeability - = (28) - = (29)
where , designate inlet and outlet pressures, respectively, the sample length, and the velocity of the fluid, respectively (see in Fig. 5). Squared pressure account for the gas compressibility. Normally, both permeabilities and should be equal to the intrinsic permeability since it is supposed to depend only on the porous medium. This is not the case since gas permeability depends on the applied pressure. According to Klinkenberg, the pressure dependence is due to gas slippage. When the size of the material pore throat radius approaches the size of the mean free path of the gas molecules, molecules that hit the pore wall are accelerated. The gas flow is not laminar any more since the velocity is not zero at the pore wall. Gas-slippage effects lead to a pressure-dependent permeability larger than that one would measure with a liquid. Klinkenberg showed that the apparent gas permeability, , is proportional to the inverse of pressure [START_REF] Klinkenberg | The Permeanility of Porous Media to Liquids and Gases[END_REF]:

22/103

= ( + ) (30)
where is the slippage factor and is the average pressure + / .

When the mean pressure increases, tends to as the gas tends to behave like a liquid. One can assume that intrinsic permeability equals liquid permeability, as it should be the permeability one would measure with liquid. It is then possible to derive intrinsic permeability from gas permeability measurements by increasing step by step the mean pressure and calculating the gas permeability for each triplet ( , and ). Plotting the gas permeability versus the reciprocal mean pressure should show a straight line where intrinsic permeability is derived from the intersect and the slippage factor from the slope. This approach is usually applied with the outlet at atmospheric pressure while inlet pressure is increased progressively. However, doing this way leads to another problem. As the inlet pressure increases, pressure drop and flow rate increase and inertial flow with turbulences may arise locally at microscopic scale. As a result, a part of the total energy provided by the pressure drop is consumed under the form of inertial flow and leads to a macroscopic non-linearity, as proposed by [START_REF] Forchheimer | Forchheimer's Therapeusis Of Internal Diseases[END_REF]:

- = + + (31)

23/103

where is the gas density, herein a triple of parameters , , are evaluated, and is the inertial factor.

Until now, de Bilbao et al. [START_REF] De Bilbao | Intrinsic permeability of refractories from gas permeability measurements: Comparison of results[END_REF] have developed an improved gas permeameter which allows to control the flow conditions to accurately determine the three parameters , , and . Although only the intrinsic permeability is of interest for this work, its accurate evaluation also needs the determination of . This work allows for gas permeability measurements both in atmospheric and backpressure modes.

I.1.5.2 Relative permeability

In general, relative permeability (dimensionless) is defined as the function of saturation such that it varies from 0 to 1 (i.e.

). As a consequence, it is often assumed to be equal to unity i. e. = or to be a cubic function i. e. = (see Fig. 6). It also can be more complicated function associated to the capillary pressure curve from van Genuchten's model (van Genuchten, 1980) (see in Fig. 7):

= √ × [( -⁄ ) ] (32)
where is the empirical parameter of the capillary pressure curve from van Genuchten's model. As illustrated in Fig. 7, the relative permeability tends to zero when the saturation tends to zero, too. The convergence is however dominated by the value of .

24/103

Assuming that the relative permeability is known, then the capillary pressure curve plays a key role in securing the optimal numerical solution of the Richards' equation (27). Hence, the transport properties of the capillary pressure curve become the main unknowns of the Richards' equation. Additionally, these transport properties can be estimated by modeling the impregnation using Finite Element Method combined with inverse identification approach [START_REF] De Bilbao | Corrosion of high alumina refractories by Al 2 O 3 -CaO slag: thermodynamic and kinetic approaches[END_REF].

I.1.5.3 Capillary pressure curve

In this section, the finite volume of a porous medium is analyzed at microstructure/ pore scale (Fig. 2). Nevertheless, it is really a challenging task to examine the capillary pressure at this scale. We therefore propose to assimilate the local scale of a pore space with a smooth vertical tube.

To determine the formulation of the capillary pressure curve, there are various approaches both in the laboratory's experiments and theoretical methods [START_REF] Dodds | Capillary Pressure Curves of Sphere Packings: Correlation of Experimental Results and Comparison with Predictions from a Network Model of 99/103 Pore Space[END_REF][START_REF] Hall | Water sorptivity of mortars and concretes: a review[END_REF][START_REF] Hunt | Basic transport properties in natural porous media: Continuum percolation theory and fractal model[END_REF][START_REF] Hunt | Unsaturated hydraulic conductivity modeling for porous media with two fractal regimes[END_REF][START_REF] Weihermüller | Erratum to "Revisiting Vereecken Pedotransfer Functions: Introducing a Closed-Form Hydraulic Model[END_REF]. However, in the scope of this thesis, we present the most known methods that have been used to determine the function of the capillary pressure curve:

(i)
The phenomenological models propose that the closed-analytic form for the capillary pressure curve is predicted by observations in order to fit experimental measurements [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF][START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1[END_REF];

25/103

(ii)

The Boltzmann transformation method assumes that the capillary pressure profile is the mirror image of the Boltzmann transform function [START_REF] Prevedello | New Analytic Solution of Boltzmann Transform for Horizontal Water Infiltration into Sand[END_REF]; (iii) Morphological pore network model constructs the capillary pressure curve from the relationship between the radius of a local pore space and saturation by considering a porous sample at the scale of pore spaces and skeletons (microstructure) [START_REF] Hilpert | Pore-morphology-based simulation of drainage in totally wetting porous media[END_REF][START_REF] Sweijen | Capillary pressure-saturation relationships for porous granular materials: Pore morphology method vs. pore unit assembly method[END_REF](iv) Fractal pores-spaces model constructs the capillary pressure curve relied on the fractal dimension by the percolation theory (Ghanbarian-Alavijeh and [START_REF] Ghanbarian-Alavijeh | Unsaturated hydraulic conductivity in porous media: Percolation theory[END_REF][START_REF] Hunt | Basic transport properties in natural porous media: Continuum percolation theory and fractal model[END_REF][START_REF] Hunt | Unsaturated hydraulic conductivity modeling for porous media with two fractal regimes[END_REF].

(i) Phenomenological models

Phenomenological models (also known as the closed-form analytical models) are lying at observations from experimental data to propose the closed-form analytical expression for the capillary pressure curve. Often these equations have been validated by comparing with the experimental measurements. Some common phenomenological models are reported in Table 1.

= { [ + ( | | )] -, < , , , , (Oh et al., 2015) = { -× [ -⁄ -] ⁄ , < × [ - ], > < , = , , , ,
To ensure that the saturation varies from the residual saturation to the maximum saturation, here the effective saturation (Eq.( 12)) is used to describe the saturation of the capillary pressure curve.

26/103

The capillary pressure curves listed in Table 1 are only reproduced when experimental data is known. For example, we apply the most common model [START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1[END_REF] to reproduce the capillary pressure curve which is then validated by comparing with the "reference" curves taken from the literature [START_REF] Prevedello | New Analytic Solution of Boltzmann Transform for Horizontal Water Infiltration into Sand[END_REF]. Here we only estimate two unknowns ( , ) by fitting the "reproduced" capillary pressure curve to the experimental one.

Let us remark that a specified phenomenological model gives a reasonable description for some but not all cases. To illustrate that, we consider two examples as follows.

• First, van Genuchten's model has been often employed to fit well with the typical form of the capillary pressure curve, as displayed in Fig. 8 (b), for porous materials (e.g. loam, sandstone and clay [START_REF] Prevedello | New Analytic Solution of Boltzmann Transform for Horizontal Water Infiltration into Sand[END_REF]); • Second, Brook and Corey's model has been often used to fit well with the typical form of the capillary pressure curve, as displayed in Fig. 8 (a), for solid materials [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF]. Recently, the closed-form analytical expressions for the capillary pressurve curve has been developed to be only a general equation, but it can be applied both to the two typical forms above (Fig. 8 (a) and (b)), as can refer to the modified van Genuchten's models (Oh et al., 2015;[START_REF] Schaap | A Modified Mualem-van Genuchten Formulation for Improved Description of the Hydraulic Conductivity Near Saturation[END_REF]). Yet, there are more unknowns to be determined in these models, providing a more complicated in estimating these unknowns by matching the proposed capillary pressure curve with experimental data.

In particular, for (Brooks and Corey's model (1964) and van Genuchten's model ( 1980)), there are only two parameters (i.e. two unknowns ( , or , , respectively) to predict the capillary pressure curve. While for the modified VG's models, there are four unknowns ( , , , [START_REF] Schaap | A Modified Mualem-van Genuchten Formulation for Improved Description of the Hydraulic Conductivity Near Saturation[END_REF] or even five unknowns , , , , Oh et al., 2015) which are required to fit the experimental measurements. Moreover, there is no evidence to confirm that the formulation of capillary pressure curve of the modified VG's models is superior to the one of van Genuchten's model or Brooks and Corey's model.

(

27/103

Therefore, in the practice, the application of the modified van Genuchten's models is extremely difficult and maybe no suited for simulations. Also, to apply the phenomenological models for reproducing the capillary pressure curve, a simple model with two unknowns should be convenient.

(ii) Boltzmann transformation function Generally, the Richards' equation is a (non-linear) Partial Differential Equation (PDE). It is therefore very difficult to obtain its analytic solutions by means of the "pure" mathematics analysis. As a mathematical approach, the Boltzmann transformation thus proposed aims at reducing the order of the PDE. Actually, it is to transform the non-linear PDE into a non-linear Ordinary Different Equation (ODE). In the general case, the nonlinear ODEs are analytically solvable in terms of the known functions combined with integrals or the Taylor series for the solutions [START_REF] Corliss | Solving Ordinary Differential Equations Using Taylor Series[END_REF].

Here, Richards' equation is transferred into a non-linear ODE by a new proposed parameter (i.e. transformation of parameter into parameter ) as follows.

= √

(33)

where and designate the position coordinate and time, respectively.

Eq. ( 33) is called the Boltzmann transformation [START_REF] Prevedello | New Analytic Solution of Boltzmann Transform for Horizontal Water Infiltration into Sand[END_REF]. Assuming that the capillary pressure is the mirror image of the Boltzmann transformation with exponent (and for the case of the extended Boltzmann transformation [START_REF] Reichardt | The recent similarity hypotheses to describe water infiltration into homogeneous soils[END_REF] with respect to Eq. (34) (Fig. 9).

Fig. 9: Capillary pressure curve (

) is assumed as mirror image of Boltzmann transformation function .

28/103

• = • (34)
Eq. ( 34) shows that the derivative of the capillary pressure curve is deduced from the extended Boltzmann transformation function. The capillary pressure profiles are then obtained analytically by substituting the derivative of the capillary pressure curve (Eq. ( 34)) into Richards' equation. In order to validate the performance of the extended Boltzmann transformation, the comparisons between its solutions versus experimental data points were done [START_REF] Prevedello | New Analytic Solution of Boltzmann Transform for Horizontal Water Infiltration into Sand[END_REF]. Nevertheless, the comparisons do not support any evidences of the hypothesis above, with whatever the value of exponent is [START_REF] Reichardt | The recent similarity hypotheses to describe water infiltration into homogeneous soils[END_REF]. Until now, employing the Boltzmann transformation to acquire the capillary pressure curve still is an ongoing idea without any validations. As a consequence, this method is not a good choice for the optimal function of the capillary pressure curve at the moment.

(iii)

Morphological pore network model

For the morphological pore network, it has been proposed to determine the capillary pressure curve by considering the pore-size distribution of a porous medium at the macroscopic scale [START_REF] Hilpert | Pore-morphology-based simulation of drainage in totally wetting porous media[END_REF][START_REF] Sweijen | Capillary pressure-saturation relationships for porous granular materials: Pore morphology method vs. pore unit assembly method[END_REF]. Therefore, the capillary pressure is here defined as the difference between "non-wetting" phase pressure and "wetting" phase pressure (Eq. ( 8)). The difference between them, i.e.

-, is required for the "wetting" fluid to invade a liquid-saturated pore unit, where a pore unit is represented as hemisphere (Fig. 10).

Fig. 10:

The hemisphere approach stands for the interfacial area between the "non-wetting" fluid and the "wetting" liquid, where a spherical presents the non-wetting fluid ( = ) -the wetting liquid ( = ) interface and is the inside radius of the hemisphere.

The capillary pressure is here defined as a function of the inside radius by employing the Young-Laplace equation [START_REF] Sweijen | Capillary pressure-saturation relationships for porous granular materials: Pore morphology method vs. pore unit assembly method[END_REF] as following.

= - = (35)
where designates the "non-wetting" fluid -"wetting" liquid interfacial tension.

The location of a pore is then determined from Eq. ( 35) combined with the erosion algorithm.

29/103

Basing on the connectivity of the "non-wetting" fluid to the "non-wetting" pore reservoir (denoted as NWP Reservoir Fig. 11), the capillary pressure curve is established [START_REF] Sweijen | Capillary pressure-saturation relationships for porous granular materials: Pore morphology method vs. pore unit assembly method[END_REF]. [START_REF] Hilpert | Pore-morphology-based simulation of drainage in totally wetting porous media[END_REF].

For more information, the reservoir is the universal terminology in petroleum engineering [START_REF] Coussy | Poromechanics: Coussy/Poromechanics[END_REF] to indicate the liquid source container (Fig. 11).

The reliability of this method strongly depends on the type of particle size distribution. For instance, the hemisphere approach only works well for narrowly distributed granular materials.

(iv)

Fractal pores-spaces model

For the morphological pore network model, fractal pores-spaces model also approaches the pore-size distribution of a porous medium to obtain the capillary pressure curve, as before. In particular, the formulation of the capillary pressure curve regarding the Pore-Solid Fractal (PSF) approach (Ghanbarian-Alavijeh and Hunt, 2012) combined with the percolation theory is given as

= -( + ) [ - ] - ; < < ( 36
)
where indicates the pore-solid interface dimension.

Eq. ( 36) is related to the pore-solid interface fractal dimension , which is calculated by considering the fractal geometry of the "front" (i.e. the boundary between the liquid-filled pores and the skeletons).

On the other hand, this model has been developed by using a bimodal pore size distribution (which is relevant to two fractal regimes) without the PSF approach [START_REF] Hunt | Unsaturated hydraulic conductivity modeling for porous media with two fractal regimes[END_REF]. However, the capillary pressure curve only fits well with experimental data in the first fractal regime and tends to the underestimate in the second regime. Here the first and the second fractal regimes are two parts of the same capillary pressure curve. Each part is plotted thanks to two different equations of two pore-solid interface dimensions and .

Obviously, for the current study in 1D case, the fractal dimension cannot be calculated. Hence, the fractal pores-spaces model cannot be applied to determining the capillary pressure curve in our case.

30/103

I.2 Synthesis of numerical methods

To reproduce the impregnation phenomenon in a porous medium, its governing equation (i.e. Richards' equation) is considered. This equation has often been numerically solved by different methods. As such, the Finite Element Method (F.E.M) [START_REF] Bergamaschi | Mixed finite elements and Newton-type linearizations for the solution of Richards' equation[END_REF][START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF], the Finite Volume Method (F.V.M) [START_REF] Lai | A mass-conservative finite volume predictor-corrector solution of the 1D Richards' equation[END_REF], and the Finite Different Method (F.D.M) [START_REF] Bouchemella | Numerical simulation of water flow in unsaturated soils: comparative study of different forms of Richards's equation[END_REF].

To obtain the optimal solutions of the Richards' equation, a suitable form of the capillary pressure curve and the relative permeability needs to be determined.

If the relative permeability and the capillary pressure curve are known, then in the literature [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF][START_REF] List | A study on iterative methods for solving Richards' equation[END_REF][START_REF] Wu | A finite-element algorithm for modeling variably saturated flows[END_REF][START_REF] Zhang | Finite Analytic Method for Solving the Unsaturated Flow Equation[END_REF] the Richards' equation is often solved by numerical method. Indeed, the analytic solutions of this equation are very difficult to achieve, except in a few simple cases [START_REF] Zhang | Finite Analytic Method for Solving the Unsaturated Flow Equation[END_REF]. Thus, this section is devoted to the numerical methods employed to obtain approximate solutions of the Richards' equation.

Generally, the tasks and stages of numerical methods to solve non-linear PDE of parabolic type can be summarized into the three main steps following [START_REF] Szymkiewicz | Modelling Water Flow in Unsaturated Porous Media, GeoPlanet: Earth and Planetary Sciences[END_REF]:

Step 1: Discretization in space and then in time. First, the discretization in space is mandatory for the classical numerical methods such as F.E.M, F.D.M, and F.V.M in order to define a number of nodes along the spatial domain. Typically, the values of unknown function at each node will be solved. The output of the discretization in space is a system of Ordinary Differential Equations (ODEs) or algebraic differential equations including the unknown nodal values and their derivatives with respect to time. Second, considering the system of the ODEs, its specified time interval will be divided into time-steps where time increment is chosen (Fig. 12). The numerical solution at the next time-step will be then obtained through the explicit or implicit methods. Step 2: Linearization. It is to notice that after the step 1, the system is the constitution of the discrete nonlinear equations. This system will be therefore transformed into discrete linear equations by using an iteration method (e.g. Newton method or Picard method or, etc.). In particular, thanks to the iteration method, an exact term will be replaced by an approximate term to calculate the values of the unknown parameters by means of the assumption, and then a number of successive corrections are calculated.

31/103

Step 3: Resolution of the system of the discrete linear equations. This system, obtained at each time of the iteration, must be solved in order to calculate the corrections of the unknowns. That is why the solutions are often called numerically approximate solutions. Depending on the treated problem, the direct or iteration methods will be considered.

Let us remark that in the framework of the thesis we will not present the basic discretization methods. The objective of this section is to point out the numerical drawbacks arising from these steps. In particular, in the following section, two main numerical difficulties will be highlighted:

1) Convergence and accuracy: Considering step 1, we first are going to present the issue of the spurious oscillations. More precisely, there is an impact of a very small time increment on the accuracy in the numerical computations [START_REF] Vermeer | An accuracy condition for consolidation by finite elements[END_REF]; 2) Speed of convergence: Considering step 2, a suitable iteration method allows to increase a number of successful corrections, which characterizes the speed of the convergence.

I.2.1 Accuracy

Considering the step 1 of discretization in space and time, the time increment is chosen in order to take into account the requirement of accuracy. In general, the accuracy of the numerically approximate solutions increases as a chosen time increment decreases. It is, however, mathematically shown that there exists a lower limit of the time increment: i.e. if the chosen time increment is less than the value of the lower limit we call "critical time increment" herein after, which results in spurious oscillations [START_REF] Thomas | Modélisation de l'imprégnation réactive, cas des céramiques réfractaires par des oxydes liquides[END_REF][START_REF] Vermeer | An accuracy condition for consolidation by finite elements[END_REF]. Additionally, spurious oscillation describes a numerical phenomenon where a numerically approximated solution oscillates around the correct one. With the aim to ensure the accuracy and the convergence of the numerical (F.E.) solutions, it is necessary to avoid the spurious oscillations by identifying and then respecting the critical time increment. In addition, it is also reported that the form of the critical time increment depends on the type of nodal element.

Let us remark that there is the equivalence of the heat diffusion equation and Richards' equation. For the heat diffusion equation, its formulation in the case of 1D is given by = ∆ (37)

where , , , , and denote the mass density of material, the thermal capacity, the temperature, the thermal conductivity and time, respectively.

Besides, Richards' equation (26) in the case of 1D (here along -direction only) is re-written as follows: Eqs. ( 37) and ( 38) are observed to be analogous where the mass density of material, thermal capacity, and the thermal conductivity are equivalent to a function of the saturation level in Richards' equation [START_REF] Karthikeyan | Numerical oscillation in seepage analysis of unsaturated soils[END_REF][START_REF] Thomas | Modélisation de l'imprégnation réactive, cas des céramiques réfractaires par des oxydes liquides[END_REF]. The different formulations of the critical time increment from the forms of the heat diffusion equation therefore are deduced as investigated in Table 2 [START_REF] Karthikeyan | Numerical oscillation in seepage analysis of unsaturated soils[END_REF].

ϕ = ( ) (38
∆ ∆ ∆ ∆ In Table 2, = ⁄ / , /
is the surface tension of the wetting liquid, (dimensionless) is the porosity and / is the hydraulic permeability.

To show the sensitivity of the critical time increment as described before (Table 2), here we give an example as follows. A numerical simulation of impregnation in porous media implemented in the F.E. software Abaqus/standard 2D axisymmetric (v6.14, 2017) (Fig. 13). First, we implement the example in two different types of nodal element (4-nodes element and 8-nodes element). The respectively critical time increment is therefore different regards to Table 2 in the case of 2D. Second, we input the time increment which is greater than the critical time increment in the case of 8-nodes element and is less than the one in the case of 4-nodes element.

In the example ((Fig. 13), the values of the unknowns are: = , × -(m/s); = ⁄ ; = .

× - ⁄ and = . where ∆ and ∆ denote space step-size and time increment, respectively. We can thus calculate the value of the critical time increment for each case of the type of nodal element following: In particular, critical time increment ∆ equals 0.01 and 1.6 for the case of 4-nodes element and 8-nodes element, respectively.

In the case of 4-node element, the time increment is chosen to be 0.01, and then this choice does not respect the critical value ( = . < ∆ = .

). Otherwise, in the case of 8nodes element, the time increment is chosen to respect the critical value ( = . > ∆ = . ).

It is to note that only the first time-steps are affected by the appearance of spurious oscillations [START_REF] Blond | Modelling of Multi-physics couplings[END_REF][START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF] This leads that it is impossible to observe this

35/103

Obviously, Fig. 14 illustrates the sensitivity of the critical time increment. A spurious oscillation (the dashed blue curve) occurs since the time-step does not respect the critical value and there are no spurious oscillations (the solid red curve) when the time increment is greater than the critical value.

In conclusion, to guarantee the accuracy (especially, avoid the issue of spurious oscillations in the first time-steps) of the numerical results, the time increment should be chosen to be small but greater than the critical increment. This inequality depends on the type of nodal element in Table 2.

I.2.2 Speed of convergence

Several "traditional" numerical methods such as finite element method, finite volume method, and finite difference method have often been used to solve Richards' equation.

Regarding the step of discretization both in space and time, these methods always require a small-time increment to guarantee the accuracy and the convergence of numerical calculations. Nevertheless, this requirement makes the computational time prohibitively large. The speed of convergence thus becomes one of criterion to evaluate the efficiency of a numerical method. A suitable iteration method is then required to increase a number of successful corrections, which characterizes the speed of the convergence. Consequently, this section is devoted to present the development of the iteration method for each numerical method to be superior in the computational effect.

First, considering the step of linearization of the F.E.M, F.V.M and F.D.M, the speed of convergence is limited by a small time-increment due to the mass-balance error. This small time-increment is chosen in step 1 to guarantee both the accuracy and convergence.

In particular, Richards' equation in the one-dimension case is given by

= ϕ • (39) Substituting = ( • ) ⁄ into Eq. (39) gives = • (40)
Regarding the step 1 above, the discretization of Eq. ( 40) is given by

+ , + - = + , • + , + (41)
where superscripts and denote time level and iteration level, and

+ , = + ,
.

Richards' equation is then transferred into the system of nonlinear discrete equations. It is thus demanded to transfer this system into the system of the discrete linear equations by

36/103

employing an iteration method. Indeed, the iteration method is proposed to approximate the right-hand term of Eq.( 41) by Newtons or Picard methods or others [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF]. Hence, the right-hand term is approximated by a new term, and the difference between the exact and approximate terms is called the mass balance error. This error then becomes the criteria of evaluating the performance of the used iteration method, i.e. if the error is acceptable "small enough", the numerical solutions will close to the exact ones. Because the mass balance error is proportional to the time increment, a small time-increment is always required to guarantee the accuracy of the numerical solutions. Hence, the mass balance error drives the speed of convergence [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF][START_REF] Wu | A finite-element algorithm for modeling variably saturated flows[END_REF].

The iteration method is consequently demanded to adapt: (1) the ability of employing a larger time increment and (2) ensuring the "acceptable" mass balance error. To sum up, the goal of the step 2 of the linearization is lying at enhancing the iteration method to minimize the mass balance error using the largest time increment.

Thus, we here report some well-known iteration methods which have often been used for this aim as the following:

1) In the 20 th century, [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF] proposed the Modified Picard's Iteration (MPI)

(also known as the quasi-Newton method [START_REF] Allen | A Finite-Element Collocation Method for Variably Saturated Flow in Two Space Dimensions[END_REF]. For this method, the approximate linear equation is obtained by employing the iteration method as follows:

, = , -+ | , - ( , -, -) + (42)
It is then demonstrated that the mass balance error is excellent and small enough both for F.E.M and F.D.M [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF].

2) In 2010, [START_REF] Wu | A finite-element algorithm for modeling variably saturated flows[END_REF]) also worked on the same goal and then proposed the Modified Head Based scheme (MHB) for only the F.E.M. Instead of using the same governing equation (Eq. ( 40)), it is suggested to approximate this equation by using the modified head-based formulation by [START_REF] Wu | Numerical simulation of saturated-unsaturated transient flow in soils[END_REF]. Then, this approximated equation can be solved by the standard iteration techniques such as Newton or Picard iteration methods. In order to evaluate the performance of the MHB scheme, it was compared with the MPI scheme over a wide time-step size and it confirmed that the MHB scheme is slightly better than the MPI scheme regards to the criteria of the mass balance error. 3) In 2014, [START_REF] Bouchemella | Numerical simulation of water flow in unsaturated soils: comparative study of different forms of Richards's equation[END_REF] considered the simplest (one-dimensional) case to deal with the numerical simulation of impregnation in porous media. By using the Finite Differential Method (F.D.M), as well as the Finite Volume Method (F.V.M) [START_REF] Lai | A mass-conservative finite volume predictor-corrector solution of the 1D Richards' equation[END_REF], the diagonal matrix which represents the system of the linear discrete equations is obtained. The stability (i.e. smooth and no spurious oscillations) of numerical solutions were demonstrated through the comparisons against laboratory measurements, simulation results from the Hydrus-1D software, or analytic solution when available.

38/103

I.3 Discussions

This synthesis above underlines that it is heavy to perform the accuracy including there is no spurious oscillations and low time-consuming resolution of Richards' equation. This observation is at the origin of the idea of developing a new numerical method free from the described numerical drawbacks. Hence, the objective of this thesis is to propose a new numerical method predicting capillary pressure profiles with no spurious oscillations and allowing to reduce the CPU time cost. In a very natural approach, we would like to introduce the gradient percolation theory, a random probability model based on the site percolation [START_REF] Nolin | A Modified van Genuchten-Mualem Model of Hydraulic Conductivity in Korean Residual Soils[END_REF]. Until now, this theory has already been applied to many numerical simulations of "real" phenomena such as the diffusion phenomenon [START_REF] Sapoval | The fractal nature of a diffusion front and the relation to percolation[END_REF] and the growth of cities [START_REF] Nguyen | Simulations of the growth of cities[END_REF]. With the inspiration from its applications recently, it seems to become an innovative approach to our problem. That is why in the next chapter, we will apply the gradient percolation theory to a numerical simulation of the impregnation process in porous media.

39/103

Chapter II Self-organized Gradient Percolation (SGP) algorithm

This chapter is devoted to the development of new numerical approach proposed to overcome the numerical drawbacks presented in chapter I. There are two main limitations arising from the chosen time-increment of the "traditional" numerical methods for solving the governing Richards' equation. The first one is that the speed of the numerical calculations is restricted by a small time-increment. The second one is that the issue of the spurious oscillations can occur in the case of a very small time-increment [START_REF] Thomas | Modélisation de l'imprégnation réactive, cas des céramiques réfractaires par des oxydes liquides[END_REF][START_REF] Vermeer | An accuracy condition for consolidation by finite elements[END_REF].

The gradient percolation approach is novel, which takes the inspiration from its recent application to the growth of cities [START_REF] Bergamaschi | Mixed finite elements and Newton-type linearizations for the solution of Richards' equation[END_REF][START_REF] Nguyen | Simulations of the growth of cities[END_REF][START_REF] Pooyandeh | A Comparison Between Complexity and Temporal GIS Models for Spatio-temporal Urban Applications[END_REF]. It is the evolution of the building constructions (in the city) called the "cluster" in the concept of the percolation theory. In particular, the "cluster" simply describes a group of connected buildings considered as population and colored in red (on the right-hand side of Fig. 16) [START_REF] Brian | Percolation Theory and Network Modeling Applications in Soil Physics[END_REF]. On the other hand, considering a numerical simulation of the crucible test (Tran et al., n.d.), the evolution of the saturation levels is reproduced by the spectrum of colors from red to blue, as displayed in Fig. 16 on the left-hand side. Consequently, the gradient percolation theory (a branch of the probability) seems enable us to simulate the impregnation phenomenon in porous material [START_REF] Adler | Real Porous Media: Local Geometry and Macroscopic Properties[END_REF][START_REF] Nguyen | Selforganized gradient percolation method for numerical simulation of impregnation in porous media[END_REF].

The growth of town is a multiple non-driven constrain problem, as the liquid impregnation is. It is why the idea to use this method has appeared. However, it is obvious that it is mandatory to modify this method in order to embed some meaningful physic concepts to ensure the reliability of the results. asymmetric (v6.14, 2017). On the right-hand side: the growth of town is modelized using the gradient percolation theory.

40/103

This chapter is organized as follows. Section II.1 first gives a thorough description of the percolation theory [START_REF] Grimmett | Percolation[END_REF][START_REF] Nolin | A Modified van Genuchten-Mualem Model of Hydraulic Conductivity in Korean Residual Soils[END_REF] both of the traditional and gradient percolation models. It is then directed toward applying the gradient percolation theory to reproducing the physical phenomenon numerically. Although it has the promising potential benefits clearly [START_REF] Adler | Real Porous Media: Local Geometry and Macroscopic Properties[END_REF][START_REF] Nguyen | Selforganized gradient percolation method for numerical simulation of impregnation in porous media[END_REF][START_REF] Sapoval | The fractal nature of a diffusion front and the relation to percolation[END_REF], it also has several disadvantages detected from the physics point of view in section II.2. Hence, this section explains our motivation to introduce an adapted version of the gradient percolation theory to deal with these drawbacks. Section II.3 works on this novel approach we call "Self-organized Gradient Percolation (SGP) method". The purpose of the simulation is to predict the capillary pressure profile (i.e. a curve between saturation and capillary pressure at every timestep).

A novel algorithm based on the gradient percolation theory to predict capillary pressure profiles, called "the SGP algorithm", is therefore established in the robust stepwise as follows:

1) At the initial time, it is found that the governing equation is the homogeneous diffusion equation, so the initial capillary pressure profile is expressed in the analytical form, the convolution integral of the Probability Density Function (PDF) of Gaussian distribution with a smoothing function. This form is taken from original idea of the "imaging process" [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF]Tropea et al., 2007) where the smoothing function is suitably chosen for the continuity and the boundary conditions of the treated problem. Also, the PDF is suggested to be slightly modified to adapt the physics meaning of interest; 2) At the next time steps, to respect the physics point of view, capillary pressure profiles are suggested to be given by means of the self-evolution process. To do that, we presume that the evolution of the capillary pressure profile is reproduced by the self-evolution of the standard deviation of the probability density function. In addition, to bridge the percolation and the local physics, the self-evolution is integrated with the capillary rise and the rise rate of the liquid [START_REF] Extrand | Forces, pressures and energies associated with liquid rising in nonuniform capillary tubes[END_REF]; 3) At equilibrium, thanks to the link above, the algorithm will be automatically stopped as there is no longer evolution of the standard deviation in agree with the physics point of view, i.e. when the capillary pressure is equal to the hydrostatic pressure [START_REF] Coussy | Mechanics and Physics of Porous Solids: Coussy/Mechanics and Physics of Porous Solids[END_REF][START_REF] Extrand | Forces, pressures and energies associated with liquid rising in nonuniform capillary tubes[END_REF]. Then, the capillary pressure profile at this state is known as the capillary pressure curve.

Finally, most of the main steps will be elaborated in section II.3.4 of the flowchart of the SGP algorithm. Also, the link between initial data of the algorithm and the physics meaning will be explained in more detail below.

II.1 Basic lines of the percolation theory

This section will thorough sketch some basics of the "traditional" percolation theory [START_REF] Grimmett | Percolation[END_REF] and the gradient percolation theory (Pierre, n.d.). The purpose of this review is, first, to provide that the percolation theory is explicitly linked to the problem of interest through using the network model. Second, it is to examine some remaining issues and challenges in the application to the interest behind the percolation theory. Finally, some

42/103

In a very natural way, such representations, the bond and site percolation models, can be applied to reproduce the impregnation phenomenon. In particular, a square of the network model represents a local pore space of a porous medium, and thus the percolation theory allows to describe only two possible states, "empty" or "occupied". There is no universal terminology to mention these states. One can just as easily say "on" or "off" [START_REF] Brian | Percolation Theory and Network Modeling Applications in Soil Physics[END_REF]. The "occupied" or "empty" states are, respectively, colored in black at random with probability or in white with the inferred probability -, independently of all other squares [START_REF] Broadbent | Percolation processes[END_REF][START_REF] Grimmett | Percolation[END_REF][START_REF] Pierre | Near-Critical Percolation in two dimensions, and related models[END_REF].

Let us consider the site percolation model (Fig. 17). The model starts as follows. A large array of squares is assimilated to a very large sheet of graph paper [START_REF] Brian | Percolation Theory and Network Modeling Applications in Soil Physics[END_REF]. It includes the set of vertices and edges. The line intersections are called "sites" and each of them has only four nearest neighbor sites. This is called a square lattice. In the thesis, we decide to employ the site percolation model since it permits us to specify one respective value (local saturation) for a square of the lattice.

As explained before, the percolation model is closely stuck to a random network model. Its mathematical model is a random probability matrix ℝ × , where ( × ) stands for the size of the matrix, designates an element of the matrix, and is a random value uniformly distributed on the interval [ , ]. A site of the lattice is declared to be occupied by the liquid with constant probability or to be empty with the inferred probability -, defined as:

• A site is occupied when ; • A site is empty otherwise.

Obviously, in this way, the current version of the percolation model cannot reproduce the impregnation phenomenon with respect to the physical meaning, since here it only permits the saturation level to be represented in two states, or rather, "occupied" or "empty". Hence, a more suitable version, the so-called gradient percolation theory, will be presented in more details in the next section.

II.1.2 Gradient percolation theory

Let us turn back to the physics meaning of the impregnation phenomenon. The saturation level (the principal state of a site) is determined by the capillary pressure, which is not constant on the lattice. Let us stress out that probability of the traditional percolation model to decide a site's state is constant. Hence, the probability must be modified to relate to the capillary pressure. Remind that the capillary pressure is the function of the level (height) with respect to the liquid regarding Jurin's law [START_REF] Coussy | Mechanics and Physics of Porous Solids: Coussy/Mechanics and Physics of Porous Solids[END_REF]. That is why we propose that the probability has to be transformed into an induced function of the location of the site. This can be done by adapting gradient percolation theory developed by Pierre Nolin [START_REF] Nolin | A Modified van Genuchten-Mualem Model of Hydraulic Conductivity in Korean Residual Soils[END_REF].

43/103

In the gradient percolation model, each site of the lattice is occupied (or empty) at random with probability that varies along the fixed direction.

More generally, probability can be considered as a function of the position, as for instance, the distance from site to the origin or a potential on Cartesian coordinate system. This probability is denoted by , being a coordinate of the site. The state of a site for the gradient percolation model is then decided as follows:

• A site is occupied when ; • A site is empty otherwise. In the example above (Fig. 18), every site of the lattice is colored in black or in white with probability = ⁄ , [ , ]. Obviously, the probability is proportional to second coordinate . This leads to the ability to have an "occupied" site at the larger height in respect of will be lower.

In the next section, some preliminary results of interest behind the gradient percolation theory will be presented and discussed with respect to the physics point of view.

II.1.2.1 Application of the gradient percolation theory

Considering a ( × ) square lattice with its plane ℤ where and are the number of columns and rows, respectively. To construct the gradient percolation model, three main steps of the procedure are as following:

• First, a matrix of random values having an uniform probability distribution on the interval [ , ] is set up; • Second, the probability function, varying along to the second coordinate

[ , ], is given as for instance:

= - + (43)

44/103

for different values of parameter ;

• Finally, a simulation result using the gradient percolation theory to give the state of each site is displayed in Fig. 19. The gradient percolation model (Fig. 19) here has two principal states. In particular, the first is that the "occupied" state is colored in black when the random value of the square is greater than or equal to its probability function . The second is that the "empty" state is colored in white otherwise.

A group of the "occupied" squares and the expansion of this group are simply called "cluster" and "cluster expansion", respectively. Actually, we observe that the cluster expansion depends on the "height" while the traditional percolation model is uniform [START_REF] Grimmett | Percolation[END_REF].

Obviously, the behavior of the gradient percolation model (Fig. 19) differs from the one of the traditional percolation model (Fig. 17). The application of the current version of the gradient percolation method to the impregnation phenomenon in porous media arises several disadvantages from the physics point of view, in Section II.2.

II.2 Limitations

From the physics point of view, the images shown in Fig. 19 do not respect the liquid rise in a solid state. There are three main problems as following.

• Problem 1: The principal states of the gradient percolation model are restricted in two states, "occupied" or "empty", which do not give meaningful results with respect to the physics of liquid impregnation. • Problem 2: Black squares of the gradient percolation model do not form a connected cluster. Obviously, there exist several black/occupied sites surrounded by white/empty sites, which also is not meaningful from the physics point of view. Indeed, by considering the impregnation phenomenon in porous media, the liquid rise in a solid state should be continuous.

45/103

• Problem 3: The cluster expansion of black/occupied sites is determined by the probability function heavily depending on parameter , but not on time. Furthermore, to agree with the physical meaning of the phenomenon, this cluster growth must be achieved without the need of imposing parameter . This requires the selfevolution process.

II.2.1 Problem 1: Multistate gradient percolation model

It is proposed to transform the 0-1 state into a continuous value state, which will be considered as the local saturation of a site , defined as a Gaussian random function:

~( , ) (44)
where is the variance, is the expectation.

The parameters of the Gaussian distribution are thus related to the expected state of site and depend on the physics.

From the physics point of view, the expectation is the probability density function of a (one directional) statistics distribution by solving the homogeneous diffusion equation only considering at the initial timestep [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF], in Section II.3.1. This distribution can be interpreted as function of the local microstructure of the porosity network.

As a consequence, the continuous value state of the gradient percolation model is generated from random values following a Gaussian distribution. In particular, will be a random value of Gaussian distribution with mean value and variance whose Cumulative Distribution Function (CDF) is denoted .

The aim is to generate the values of , distributed following a given Gaussian distribution. To do that in our study, we use the inverse transformation method for its generality and simplicity [START_REF] Rick | Simulating Data with SAS[END_REF]:

We assume that there exists a strictly monotone transformation

: [ , ] → ℝ, such that
= where is a uniform distribution on the interval [ , ].

Then,

= ℙ = ℙ = ℙ( - ) = -( - ) (45)
Therefore, random values can be generated from -( - ) in two main steps:

46/103

• The first one is the setup of a matrix of random values, which is given from a uniform distribution on the interval [ , ]; • The second one is computing = -( - ) in accordance to , . In the numerical practice, this one can be solved using the dedicated algorithm "norminv" developed in Matlab computing software (The Mathworks Inc., Natick, Mass, USA). In the numerical practice, to apply the two steps above for solving the problem 1, we introduce an example below.

We now proceed to consider a × square lattice. Each square of the lattice represents a local saturation. In particular, the principal state varying from the maximum saturation to the residual saturation is colored in the spectrum of colors from black to white, respectively.

For instance, in Fig. 20, it is proposed a simulation result for the multistate gradient percolation model in the case where the probability density function only depends on the second coordinate of following:

= -+ 2 × -⁄ ; = (46)
where , , , are unknowns and depends on : = -⁄ .

47/103

We observe that the behavior of the spectrum of colors from black (at the bottom) to white (at the top) is characterized by the three parameters , , , Table 3.

Table 3:

The effect of the parameters to the multistate GP model.

Parameters Characterizations

Dimension of the square lattice.

The number of black squares depends on the value of . This parameter is related to the permeability of the wetting fluid Standard deviation of random values at the same "height" . It is strongly related to the homogeneity of the support.

In conclusion, to obtain the multistate gradient percolation model we use an inverse transformation approach; the probability density function is chosen to respect the physics point of view. This model allows to reproduce the local saturation level varying from to , and it therefore confirms that the first limitation is overcome.

Considering problem 2, let us remind that there exist several white/empty squares surrounded by black/occupied squares (Fig. 20), meaning that a local pore can become of maximum saturation while it is surrounded by empty pores. This is not suitable for the physical meaning of our application. It is thus necessary to address this issue.

II.2.2 Problem 2: Convolution operator

We will now address the issue of the continuity of the SGP model. From the physics point of view, pressure of a local pore space is influenced by the one of the closely surrounding local pore spaces. To include the continuity to our model, the convolution procedure from the original idea of "the imaging process" [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF] is employed.

Let us recall [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]

: = { , if , otherwise (50)
and is the number of sites .

Obviously, and are defined on the ( × ) square lattice, (i.e. the bounded domain).

Hence, the convolution product between a local saturation of the gradient percolation model and the smoothing function is well defined.

In the numerical practice, we put = , ℝ which is the index of column and row considered in the square lattice, respectively. The convolution is then given as following:

= * = * , = ∑ -′ , -′ • ′ , ′ ′ , ′ ℝ 2 (51)
We then employ the smoothing function in Eq.( 49) to construct a smooth approximation for (the unified expression of Eq. ( 51) and Eq. ( 49)):

= ∑ -′ , -′ ′ , ′ = ∑ -′ , -′ ′ • ′ = , | ′ |+| ′ |≤ ′ , ′ (52)
The convolution procedure (Eq.( 52)) can be described in the representation below (Fig. 21). For example, the "correlated" gradient percolation model (which is simply a continuous gradient percolation model) (Fig. 22 on the right-hand side) is derived from the "uncorrelated" gradient percolation model (Fig. 22 on the left-hand side) through the convolution operator with the smoothing function (Eq.( 49)). The convolution procedure in Eq. ( 48) ensures the connectivity of the "wetting" cluster. Indeed, Eq. ( 48) is the convolution of the local average saturation with the smoothing point-spread function . To implement/construct the convolution procedure, the different steps are orderly proposed below:

49/103

50/103

Step 1:

Compute the values of the local average saturation given in Eq.( 44) numerically. In general, the local averaged saturation is derived by Gaussian distribution with mean and variance . In addition, in this work, is the local average saturation profile.

For example, this profile (Eq. ( 53)) is the so-called exp-function as:

= + ( ⁄ ) ⁄ - (53)
where is the vertical direction, is time (seconds), and are constant and can be understood as the "reference" material parameters linked to the local microstructure.

Step 2:

Determine the smoothing point-spread function. The smoothing function is associated with the fact that the pressure of a site is affected by the pressure of the surrounding neighbors. This can allow the transition from the discrete gradient percolation model to the continuous one. To achieve this point, the smoothing function should be defined by Eq. ( 49).

Step 3:

Compute the convolution of the local average saturation with the smoothing pointspread function in a numerical simulation. The result is displayed in Fig. 22 on the righthand side:

= * = ∑ -′ ′ ′
, and ′ are sites (54)

The computational cost indicates how much computational effort has to be invested in the convolution procedure for each site and is divided at:

(i)
The first computational effort of determining all sites belong to set regarding Eq. ( 49). This step is therefore devoted to the determination of set , i.e. the neighbors of each site ; (ii)

The second computational effort aims at computing the convolution of each site with the surrounding sites.

Regarding the CPU time cost of the numerical simulation, it remains expensive, even impossible, to obtain the local saturation for all sites in the case of an industrial application which include a very large number of sites. This computational challenge is overcome by the gradient percolation method, being expected to reduce the computational time. To compute Eq. ( 54) avoiding that the number of sites to be prohibitive, both fast Fourier transformation ℱ and inverse Fourier transformation ℱ -are employed.

(capillary pressure (or height) versus saturation), the initialization uses the (mirror) curve that stands for the saturation versus height graph. We assume that the capillary pressure curve can be derived by a probability density function (not necessarily to be the Gaussian distribution) depending on parameter . Then, the standard deviation of the proposed density function quantifies the height distribution regarding the saturation. It follows that a same local saturation can correspond to different heights, according to a statistical distribution. 64).

Note that in the SGP model, the value of density function differs from local saturation at each site . Following the physical meaning for these notations, and can be understood as:

•

is the ratio of the porous volume, occupied by the liquid, without the effect of the surrounding sites;

• is the ratio of the porous volume, occupied by liquid, accounting for the effect of the surrounding sites. Additionally, and are also known to be the fractions of "allowed" pores and "accessible" pores at each site, respectively [START_REF] Zheng | A new capillary pressure model for fractal porous media using percolation theory[END_REF]. Thus, it can be assumed that the fraction of "accessible" pores represents the local saturation.

To identify the evolution of the density function in time, it is proposed to reproduce the evolution of the density function in time by the evolution of the standard deviation of the probability density function (graphic illustration of the initialization in Fig. 23).

At each timestep , the standard deviation function will be denoted by (

) and the evolution of the standard deviation will be denoted by ( , ). The standard deviation at timestep can be thus calculated from the one at the following timestep as:

= -+ ( , ) (65)
It is to remark that at each timestep, (constant) standard deviation by default differs from standard deviation function (

).

Fig. 25:

The capillary pressure and the hydrostatic pressure allow to calculate the maximal "height" (i.e. the height at the equilibrium state), denoted by ℎ ( , ), and the "height"

during the process with respect to the liquid at a given saturation, respectively.

We now turn to Eq. ( 76) that allow us to obtain the steady state as suggests the following example:

A test of the interest of Alumina 99% (AL25) and Glycerine, which is simulated both in the SGP and finite element methods. The initial values of the test are reported in Table 4.

Parameters

Values Units Eq. ( 76) induces that the simulation will stop when there is no difference between capillary pressure and hydrostatic pressure (i.e. ( ℎ -) = ), which implies the end of the evolution of the standard deviation. This should correspond to the steady state, in agreement with the physics of capillary rising. ).

60/103

This proposition has been applied to the SGP algorithm with the CPC obtained through van Genuchten's model (considering = (Pa) and = . (dimensionless)). Fig. 27 shows the evolution of the standard deviation with time increment. This figure shows that the standard deviation evolution tends toward , corresponding to the limit steady state.

We conclude that the assumption described by Eq. ( 76), allows to access a physically meaningful description of the evolution of the capillary pressure profiles in time.

Before computing the local saturation on the boundaries, we have to determine the type of boundary conditions applied; the convolution operator (more precisely, the associated smoothing function) is also used to define the boundary conditions. This is discussed in more detail in the following section.

II.3.3 Boundary conditions

We previously explained the first goal of the convolution concerning the continuity of the liquid impregnation in porous media. The second goal of convolution is to compute the local saturation both on the Dirichlet and Neumann type boundary conditions. Following the SGP model, it is assumed that the matrix/ lattice represents physical domain ⊂ ℝ and is bounded by , decomposed as:

= ( 78
)
where and are the boundary with the Dirichlet and Neumann type boundary condition, respectively; these boundary conditions are not overlapped, i.e.

= ∅.

Concerning the impregnation process in porous media, several boundary conditions related to the physics exchange made on the concerned interface can be considered. We distinguish two main conditions: liquid can or cannot flow out of the boundary, i.e. the constant air pressure and the draining condition or the undrained condition, respectively. Herein, the draining condition is divided into four different types [START_REF] Ahmed | Equilibrium of droplets on a deformable substrate: Influence of disjoining pressure[END_REF]:

• The first one is free drainage. This condition is to characterize that the liquid can flow out of the porous material because of the impact of the gravity and the capillary pressure.

• The second one is the wetting of the porous material. A dry porous sample is supposed in contact with the liquid at its base (i.e. the bottom or the top with respect to the test worked on).

• The third one is the forced drainage. The liquid is assumed to be added to the top of the porous column that produces a travelling wave.

• The final one is the pulse drainage. The purpose is to describe a small volume of the liquid which is injected at the top of the sample and left to evolve.

Here we are dealing with the following two main questions.

61/103

• What is the meaning of these boundary conditions regarding the physics point of view?

• And what is the kind of the boundary condition we should apply to the specific surface we are interested in? We suppose that the (cylindrical) porous sample is in contact with the liquid at the bottom surface. The liquid can start to impregnate through this surface and rises up. We consider a nondraining condition on the circumference of the cylinder. Finally, we impose air pressure at the top surface. Hence, if the liquid has enough kinetics, the liquid can drain out of the top surface.

We now proceed to consider two cases: (A) the liquid can flow out of the boundary surface; (B) the liquid cannot flow out of the boundary surface. And we are looking for the method that permits us to impose the right boundary condition using the SGP numerical simulation. The aim is to develop the method to a suitable boundary condition so that the liquid can drain out of the boundary surface. There are two possibilities: either the free drainage boundary condition can be employed or the outside pressure is imposed in order to generate a difference of pressure and control the liquid drain out.

Note that, the value of the outside pressure affects the result of the solution. A test where the outside pressure is imposed on the top surface of the cylinder sample is drafted in Fig. 28. We now proceed to analyze the impact of the outside pressure (here denoted by ) with different values on the mass gain curve.

63/103

The maximum value of capillary pressure (which will equal the reference pressure at equilibrium in Fig. 29 (b)) is compared to the outside pressure : if there is no difference between them, then the liquid cannot drain out of the top surface. Obviously, the outside pressure should be imposed to be the reference pressure of the CPC. For example, in Fig. 29 (a), the larger outside pressure is, the earlier the capillary pressure profile is stopped, the smaller will be the total mass gain. In other words, the outside pressure acts on the steady state of the solutions.

(B) The liquid cannot flow out of the boundary surface.

The goal is to impose a boundary condition denoted "the undrained boundary condition" such that the liquid cannot drain out of the boundary surface.

For the liquid transport equation, the model equation on the time interval , ] with > , is given as [START_REF] Aggestam | Numerical model reduction with error control in computational homogenization of transient heat flow[END_REF]:

+ ⃗ = , , , × , ] (79)
where is a parameter, ⃗ is the flow density vector, and , is a source term.

The boundary is then divided into:

64/103

Considering the treated problem from the physics point of view, the top and vertical boundary surfaces are in contact with the outside pressure. It is thus suitable to impose the outside pressure as the origin of the boundary condition on these surfaces. However, to ensure that the liquid in the interior region cannot flow out of the vertical surfaces, it is sufficient that we take the boundary condition of that flux equal to zero on these surfaces whatever the outside pressure is (Fig. 30).

It is now necessary to "translate" these boundary conditions in terms of the convolution procedure.

The undrained boundary condition is considered first. We impose that the liquid doesn't flow out of the boundary and, furthermore, no-impact of the pressure outside is considered. Thus, to reproduce this point by the means of the convolution, the saturation of each site on the boundary surface is defined by the mean value of that of neighbors inside the matrix. Then, the smoothing function is rewritten in the following way:

= ∑ (80) where = { ℤ × : | | }, ×
is the size of the matrix, and is the number of sites .

The second type of boundary condition corresponds to the case that the liquid can flow out of a boundary facing the impact of a constant outside pressure. To stay with the spirit of the convolution procedure, the saturation of a site on the boundary is computed as Eq. ( 54) combined with Eq. ( 80); a column/ row is thus added outside the boundary to bound the matrix. Furthermore, outside pressure here considered is a constant air pressure, meaning that the added column/row should have the same value for saturation [ , ].

Finally, the drained boundary condition is considered. It is imposed that liquid can flow out of a boundary without the impact of outside pressure. As before, a column/ row is added outside the boundary to bound the matrix. To reproduce the no-impact of outside pressure, it is assumed that the column/row has the same zero value for saturation.

Herein, a porous medium (that is represented by the matrix/ the lattice by the means of the SGP method) can be understood as physics domain with boundary (Fig. 31). To impose the different boundary conditions on each boundary surface, the boundary domain is partitioned as the following:

= (81)
corresponding to the three types of boundary conditions above. At the initial timestep, depending on the treated problem, it is considered that the top/bottom surface of the porous sample is in contact with the surface of liquid. This leads to impose that the last row/ first row (to fix ideas take the first one) of the matrix is fully filled with liquid, respectively. To this Three boundary conditions are applied on each different interface. Of course, the choice of the boundary conditions depends on the interpretation of the physics phenomena.

II.3.4 The SGP algorithm

The SGP algorithm is illustrated in Fig. 32. The computation starts by the calculation of the initial value of the local averaged saturation . Following the gradient percolation theory, site saturation is then computed. In order to obtain an homogeneous field, the convolution operator is used for computing local saturation .

The physics problem introduced in the subsection II.3.2 is here characterized by the parameter ( , ) which is a function of , , ℎ . The value of ( , ) is added to standard deviation -. Obtained standard deviation allows the calculation of ( , )

yields the calculation of ( , ). The loop on timestep continues until obtaining a null value of ( , ). This null value designates that there is no more evolution of the standard deviation and the algorithm will stop.

67/103

II.3.5 Identification of SGP model

The flow chart proposed on Fig. 32 depicts the different time steps of the SGP algorithm. It requires 5 input data ( , , , ) and 3 initialization data ( ℎ , , ). In general, these parameters, , , ), are linked to the geometry of the porous media, to the fluid properties and the Capillary Pressure Curve, respectively. Unfortunately, the experimental estimation of the values of parameters , and may is beyond our experimental capabilities.

Hence, these parameters, , , and of the input data, have been resumed in only one = 2 ⁄ . It has also been assumed that the initial standard deviation and hydrostatic pressure ℎ are very close to .

Consequently, in order to use the SGP algorithm, the number of unknowns is reduced to only three parameters ( , ,

). Note that:

(i) depends on the capillary diameter and the liquid properties; (ii) ( , ) characterize the form of the capillary pressure curve.

The appropriate values of these parameters have been determined by an inverse identification method: a vertical impregnation test has been modeled in the commercial F.E.M code ABAQUS/Standard (v6.14, 2017), and then the simulated mass gain curve has been compared to the numerical results of the SGP method.

The evolution of the standard deviation (Eq.( 76)) at a given saturation * is then given by

( , ) = * ℎ * - (83)
where:

• * is a given saturation and constant timestep ;

68/103

• * is the CPC of saturation * . It can be taken from van Genuchten's model as

(van Genuchten, 1980): * = ( * -⁄ -) (84)
where and are the pressure reference and an empirical parameter, respectively. Besides, parameter affects the form of the capillary pressure profile. Hence, we presume = ⁄ , a link of parameter with parameter (Eq.( 84)).

• ℎ * is the hydrostatics pressure of saturation * at timestep as [START_REF] Coussy | Mechanics and Physics of Porous Solids: Coussy/Mechanics and Physics of Porous Solids[END_REF][START_REF] Extrand | Forces, pressures and energies associated with liquid rising in nonuniform capillary tubes[END_REF]:

ℎ * = * , (85)
To identify unknowns , , , * , we employ the inverse identification method: the test (for the hidden variables) is simulated by the finite element method (Abaqus software). The simulated mass gain curve (Fig. 34) is then compared to the curve from experiment. The aim is to fit both (experimental and simulated) mass gain curves by optimizing these unknowns, as explained below The inverse identification method is divided into three main steps as follows.

Step 1: Identify the values of and . The purpose is to optimize the values of and by matching both values of the total mass gain. In particular, Eq. ( 84) combined with Jurin's law [START_REF] Coussy | Mechanics and Physics of Porous Solids: Coussy/Mechanics and Physics of Porous Solids[END_REF] is thus given by 69/103

ℎ = ( + ⁄ ) - (86)
The formulation of the total mass gain (denoted by

) is therefore obtained as

= ∫ ( + ⁄ ) - ℎ (87)
where is the mean radius of the pores determined by mercury intrusion porometry.

The values of and is fully reported in Table 5, which gives a good agreement.

Table 5:

The chosen values of and .

(Pa) = ⁄ and = (dimensionless)

Total mass gain (g) Eq.( 87) Given value .

. .

As a consequence, we conclude that step 1 characterizes the total mass gain. It is to note that the chosen values of and in Table 5 are fixed when the time step goes on.

Step 2: Identify the value of parameter . It is to optimize the value of parameter through matching the shape of both curves (Fig. 35).

70/103

The value = × - ⁄ gives a very good agreement. We therefore conclude that step 2 characterizes the kinetic of the mass gain curve.

Step 3: Identify the value of parameter * . It aims to optimizing the value of parameter * by matching both capillary pressure curves, i.e. the shape of the mass front (Fig. 36). The value * = .

gives a good agreement, though slightly different at the beginning of the curve. We thus conclude that step 3 characterizes the shape of the mass front (or rather, the CPC).

We here reported (Table 6) the optimized values of unknowns ( , , , * giving a fairly good agreement using the inverse identification method.

Parameter Value Unit Characterization

Step 1 .

Total mass gain

Step 2

• - ⁄ Kinetic of mass gain Step 3 * .

Shape of the mass front

In order to validate these unknowns, we now need to compare capillary pressure profiles of the SGP method (where the optimized values of Table 6 are taken into account) versus corresponding curves exported from the (F.E.) Abaqus (v6.14, 2017). Fig. 37 shows the results at the different time steps. The variables and represent time related to the self-organized gradient percolation method and Abaqus simulation, respectively. We observe that these optimized values of Table 6 give a very good agreement. This confirms the validation of the optimized values by employing the inverse identification method.

71/103

II.4 Discussions

In this chapter, we use the gradient percolation theory in order to reproduce the impregnation phenomenon in porous media. In the classical approach the computer simulations of this application do not respect the physical parameters. For instance, let us mention the cluster connexity issue (existence of several black/occupied sites surrounded by white/empty sites the badly suited cluster expansion). We hence modified the method to overcome these limitations. A new algorithm has been proposed we call "Self-organized Gradient Percolation (SGP) algorithm".

The output of the SGP algorithm, is the capillary pressure profiles, when the gradient percolation theory has been used to calculate the local saturation. Also, input data required for the algorithm is linked to the physics thanks to the identification between the initialization of the algorithm and the homogeneous diffusion equation. This leads to the initialization is the convolution of the probability density function of Gaussian distribution with a chosen smoothing function. We considered two associated functions in a meaningful way with respect to the physics of the problem. The first one is the probability function with time-dependent mean, variance term and exponent. The variance term is related to the driving force in order to satisfy the self-evolution process. The exponent enables us to switch the type of the associated capillary pressure curve from van Genuchten's model to Brooks and Corey's model or others.

72/103

Furthermore, a suitably chosen smoothing function allows to ensure the spatial continuity and the boundary conditions of the treated problem (as well). Indeed, we looked for the method that allows to impose the right boundary condition using the SGP numerical simulation. Finally, the SGP algorithm is constructed completely.

In the next chapter, we would like to validate the performance of the SGP algorithm by comparing analytical, experimental and numerical results. The purpose is lying at this validation to confirm that the CPU time cost is reduced and there are no spurious oscillations.

73/103

Chapter III Application and validation

This chapter deals with the validation of the performance of the SGP method. This is done by comparing its results with results obtained both from classical numerical and experimental approaches.

In consistence to the remainder of the thesis, the chapter is structured in five major sections. In first section III.1, the setup for the experimental tests is presented. In the section III.2 and III.3 the assumptions for the Finite Element Model and the SGP model of these tests are detailed. In the forth section III.4, a first qualitative comparisons between numerical results and experimental results are assessed. In the section III.5 the results from the two numerical methods are compared. Then, in the final section, a discussion of all the results and the conclusion regard the SGP performance is done.

III.1 Description of the experiments

In this thesis, focused on numerical and mathematical development, the experimental setup and results comes from the literature [START_REF] Coffignier | Mesures et simulation de l'imprégnation non réactive d'une alumie poreuse par un liquide[END_REF][START_REF] Davidson | The measurement and description of water flow through Columbia silt loam and Hesperia sandy loam[END_REF][START_REF] De Bilbao | Identification of Transport Properties of Refractories: Intrinsic permeability and Capillary Pressure Curve[END_REF][START_REF] Prevedello | New Analytic Solution of Boltzmann Transform for Horizontal Water Infiltration into Sand[END_REF]. All the tests are summarized in (Table 7).

The experimental setup for tests 1 to 3 is to perform the tests of the infiltration process into a homogeneous soil. In particular, the soils used in tests 1 to 3 are Marine sand, Columbia silt loam, and Hesperia sandy loam, respectively.

For test 1 [START_REF] Prevedello | New Analytic Solution of Boltzmann Transform for Horizontal Water Infiltration into Sand[END_REF], samples of Marine sand, made up of 98% sand, 1% silt, and 1% clay, were collected near Matinhos, PR, Brasil. For the infiltration tests in the laboratory, the sample of Marine sand is packed and then put into an acrylic plastic cylinder with 0.022 m in diameter and 2.5 m in height. Additionally, the reservoir was again and again fed by water. This leads to the fact that the excess was drained away, and the water level at atmospheric pressure remained at the level of the free sand surface of the semicircular acrylic container. The measuring time of the test was about 36h assuming that the steady-state is obtained. The recorded data concerns the soil's water profile evolution with time, also known as the capillary pressure profile.

For tests 2 and 3 [START_REF] Davidson | The measurement and description of water flow through Columbia silt loam and Hesperia sandy loam[END_REF], samples of Columbia silt loam and of Hesperia sandy loam were collected at the border of the Sacramento River and in Bakersfield area, respectively. For the horizontal and vertical water infiltration tests, the sample of Columbia silt loam/ Hesperia sandy loam is packed into a clear plastic cylinder with 3.2 cm in diameter and 1 cm in wide sections. The water flow was measured thanks to a graduated cylinder, as shown in Fig. 38. The pressure of water entering the soils was controlled by a fritted glass bead plate. The measuring time of the experiments of Columbia silt loam and of Hesperia sandy loam is about 74/103 24h45min and 13h30min, respectively. As before, the recorded data are soil's water profile evolution with time, also known as the capillary pressure profile. Fig. 38: Experimental setup for tests 2 and 3 in order to measure both horizontal and vertical water infiltration processes [START_REF] Davidson | The measurement and description of water flow through Columbia silt loam and Hesperia sandy loam[END_REF].

Capillary rising test (i.e. tests 4 and 5 in Table 7) is an alternative test for the determination of the associated unknowns of the capillary pressure curve. The experimental setup has been developed by Emmanuel de Bilbao at the laboratory CEMHTI, University of Orléans, France [START_REF] Coffignier | Mesures et simulation de l'imprégnation non réactive d'une alumie poreuse par un liquide[END_REF] who is a partner of this thesis. The sample is held with a specifically designed ring above a liquid contained in a crystallizer (see Fig. 39 a).

The sample is a cylinder with a diameter of 35 mm and a height from 40 up to 60 mm. The lateral surface of the sample is coated with wax to ensure the 1D assumption. The bottom of the sample is slightly immersed in the liquid. A crystallizer is placed on the pan of a balance connected to a PC and the mass is recorded during the test. The mass of liquid contained in the crystallizer decreases while the liquid penetrates the sample during the capillary rising. Neglecting the liquid evaporation during the test, the mass loss of the crystallizer is equal to the liquid mass that has impregnated the porous medium. The test is performed at room temperature and glycerin or engine oil are used to prevent evaporation. The immersion depth is defined from the sample porosity, the size of the sample and crystallizer diameter so that the sample stays immersed during the whole test (see Fig. 40). The crystallizer diameter has to be as large as possible, in comparison with the sample diameter, to assume that the pressure at liquid/sample interface and Archimedes' force are constant. A crystallizer with a diameter of 115 mm has been chosen.

75/103

Fig. 39: a-The capillary rising test [START_REF] De Bilbao | Identification of Transport Properties of Refractories: Intrinsic permeability and Capillary Pressure Curve[END_REF][START_REF] De Bilbao | Corrosion of high alumina refractories by Al 2 O 3 -CaO slag: thermodynamic and kinetic approaches[END_REF]; b-The principle of the mass gain curve.

The results for these tests are presented hereafter for porous sample made up of a pure alumina refractory , in particular, Degussit AL25 was synthesized by Umicore for test 4 [START_REF] De Bilbao | Corrosion of high alumina refractories by Al 2 O 3 -CaO slag: thermodynamic and kinetic approaches[END_REF][START_REF] Mayer | Oxide-Ceramic Products for High-Temperature Technology[END_REF], and Alumina LCC for test 5 [START_REF] De Bilbao | Identification of Transport Properties of Refractories: Intrinsic permeability and Capillary Pressure Curve[END_REF]. Assuming the steady-state is reach, the measurement have been done for 30,000 and 300,000 (seconds) for the tests 4 and 5, respectively.

For the test 4, the porosity of Alumina 99% (AL25) was around 20%. The mean diameter of the pores, measured by mercury intrusion Porometry, was 14

. The three transport properties , , have been determined by inverse identification method.

For test 5, the intrinsic permeability and open porosity values of Alumina LCC were experimentally identified to be . × - and 0.19%, respectively (thanks to master degree cooperation between RWTH Aachen, Germany and University of Orléans, France). The mean diameter of the pores is determined to be 0.75 by mercury intrusion porometry. Moreover, the experimental values of the permeability and the porosity have been checked by comparison with numerical data obtained through inverse identification [START_REF] De Bilbao | Identification of Transport Properties of Refractories: Intrinsic permeability and Capillary Pressure Curve[END_REF].

The experimental data from the capillary rising test is the mass gain versus time that permit to estimate the mean saturation of the pores as follows:

III.2 Finite Element models of the experimental tests

For the FEM model, we proceed to implement tests 1 to 3 in Hydrus-1D (v4.16) and tests 4 and 5 in Abaqus (v6.14). The most important parameters of these simulations are the associated parameters of the capillary pressure curve (i.e. pressure reference or and empirical parameter or ). To determine their values we employ well-known models from the 77/103 phenomenological method: van Genuchten's model [START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1[END_REF] for tests 1 to 4 and Brook's model [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF] for test 5. These choices are justified by the hydraulic behavior exhibited by the materials.

• Van Genuchten's model

The capillary pressure curve, derived from van Genuchten's model [START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1[END_REF] is described by:

= ( -⁄ -) (90)
where the quantities of the CPC are the pressure reference (Pa) and the empirical parameter (dimensionless).

Eq. ( 90) is then integrated with the F.E. ) code for tests 1 to 3 and with the F.E. Abaqus (v6.14) code for test 4. The values of the quantities for tests 1 to 3 are taken from the literature [START_REF] Davidson | The measurement and description of water flow through Columbia silt loam and Hesperia sandy loam[END_REF][START_REF] Prevedello | New Analytic Solution of Boltzmann Transform for Horizontal Water Infiltration into Sand[END_REF]; for test 4 the parameters of the model are estimated through an inverse identification method [START_REF] De Bilbao | Identification of Transport Properties of Refractories: Intrinsic permeability and Capillary Pressure Curve[END_REF].

• Brooks and Corey's model

The capillary pressure curve, derived from Brook and Corey's model [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF], is given by:

= (91)
where the quantities of the capillary pressure curve are the pressure reference (Pa) and the empirical parameter (dimensionless). This model was used only to simulate test 5. The identification of the parameters was carried out by inverse identification based on F.E. Abaqus/ Standard simulation.

III.2.1 The F.E. Hydrus simulation for the tests 1 to 3

The F.E. Hydrus (v4.16) package for the 1D-infiltration of water into soils consists of two main stages (Fig. 41):

1) Pre-processing. The goal of this stage is to take into account the material properties, capillary pressure curve, time increment, boundary conditions, etc.;

In the framework of this thesis, the numerical simulation for tests 1 to 3 have been implemented in the Hydus-1D Software (v4.16) as follows.

78/103

-The liquid for these tests is water. In particular, this has been automatically taken into account the Hydrus-1D software by selecting the "Water Flow" option of the "Main Processes" dialog window; -In the "Geometry Information" dialog window, main parameters are the dimensions of the porous sample (Table 9). Also, the horizontal or vertical direction of the fluid flow is determined through using 0-1 identification, respectively, in the "Decline from Vertical Axis" dialog; -In the "Print Information" dialog, the selected-output results of the simulation are the capillary pressure profiles; -The maximum time and time increment of the simulation are chosen to avoid the issues of spurious oscillations [START_REF] Vermeer | An accuracy condition for consolidation by finite elements[END_REF]; -In the "Soil Hydraulic Property Model" dialog window, the type of the capillary pressure curve deduced from van Genuchten's model has been used for these tests (Table 9); -For the boundary condition, the free drainage and variable pressure boundary conditions are imposed on the bottom and the top surfaces, respectively. The Hydrus-1D software is only used to solve the Richards-equation based on soil water flux. That is why we decide to use the numerical F.E. Abaqus/ Standard software for the numerical simulation of test 4 and 5 for which the liquid is not water. The FEM model consider a quasi-one-dimensional case (Fig. 42). In Fig. 42, designates the height of the porous sample under the gravity field. For tests 4 and 5, the porous sample (Fig. 42) is of 0.04 m in height. An homogeneous mesh size with a ratio close to one has been chosen to obtained comparable results with Hydrus computations. The values of the material properties (both for porous material and liquid) used for the F.E. simulation are reported in Table 8.

As illustrated before, the choice of the time increment with respect to the element type is very important to avoid spurious oscillations. Indeed, here the used finite element type is 4node elements (denoted "CAX4P"). As mathematically proved, the value of time increment must satisfy the following inequality [START_REF] Karthikeyan | Numerical oscillation in seepage analysis of unsaturated soils[END_REF]:

> (92)

80/103

where stands for relative permeability (dimensionless), a cubic function of the saturation in this simulation, ⁄ is the hydraulic conductivity of the liquid, ⁄ is the rate of the change of the saturation regarding pore pressure, and is the mesh size. For the boundary condition, two types are used in these simulations as follow. Here the fluid from the reservoir is always in contact with the bottom surface of the sample, so the first boundary condition corresponds to a free draining. On the top of the sample, there is a (constant) air pressure imposed.

For more information, the choice of these boundary conditions has been based on the example of a typical damn (Anon, 2011). As described in Abaqus manual, the draining boundary condition/ the drainage-only flow type label (denoted "QD") condition indicates the velocity on the treated surface. One is given by

= { -∞ ; when > ∞ ; otherwise (93)
where is the pore-fluid velocity directed outward normal to the surface, is the seepage coefficient, and the bi-parameters , ∞ indicate the current pore pressure and a reference pore pressure value, respectively. Considering the special case ∞ = the drainage-only flow type label (denoted "Q") condition is then given by = { ; when > ; otherwise (94) Note that only the first timesteps are affected by the issue of spurious oscillations [START_REF] Vermeer | An accuracy condition for consolidation by finite elements[END_REF], i.e. it is impossible to observe this problem at the steady-state. It is thus suggested to observe the results obtained through the SGP method at the first time-steps. As previously explained, the initialization of the SGP method is done by an analytical formulation. For the further time-steps, following the SGP assumptions, the evolution of the capillary pressure profiles is that of the standard deviation of the PDF, so the analytical formulation of the results is ensured. This point ensures that there are no spurious oscillations. To illustrate the ability of the SGP algorithm regarding the effect of small mesh size, the simulated capillary pressure profile at a timestep is given with different small mesh sizes in Fig. 47.

Fig. 47:

The measured and simulated capillary pressure profiles at the tenth time -step obtained through the SGP method with different small mesh sizes.

90/103

Fig. 47 shows the measured capillary pressure profile and the simulated one by using the SGP method at the tenth timestep. It is to notice that the capillary pressure profile from the SGP method is exact and has the analytical form of the PDF of a statistics distribution (i.e. it is not an approximate numerical solution deduced from the exact one (mesh elements)).

In the SGP method, obviously there is the discretization of the analytical equation of the capillary pressure profile along the abscissa direction, and thus it can be seen in Fig. 47 in the case of too small mesh sizes, the capillary pressure profiles are denser and smoother regarding the analytic form.

In conclusion, the SGP algorithm fulfills the objectives of this work: (i) The solutions of the SGP method are free of spurious oscillations and (ii) the computational time is reduced.

91/103

Chapter IV Conclusions and Outlook

IV.1 Conclusions

In this thesis, we have developed a novel algorithm we call "Self-organized Gradient Percolation (SGP) method" for reproducing the non-reactive impregnation process in unsaturated porous media numerically. The main motivation of the new method is to overcome the numerical difficulties in terms of computational time and accuracy.

In chapter I, we reviewed research both on theoretical and numerical analyses within the framework of the liquid transport in a porous medium at the macroscopic scale. The state variables of a porous medium at the local scale (e.g. porosity and saturation) were therefore introduced. Furthermore, its governing equation is yielded by instituting the Darcy's law in the unsaturated case into the mass conservation equation, also known as the Richards' equation. Herein, the main unknown is the optimal formulation of the CPC. We thus investigated several well-known methods from the literature, which confirms that the phenomenological model is the most suitable for our problem. We also reviewed the associated numerical methods which have often been used to solve the Richards' equations in order to highlight two important features. There are the analysis of the issue of spurious oscillations and the time consuming.

In chapter II, we proposed to use a probabilistic method based on gradient percolation. In particular, the gradient percolation theory has been employed to calculate the local saturation. The main idea is to directly use the capillary pressure profile as a probability density function in the gradient percolation method. The evolution of the standard deviation is linked to capillary pressure to ensure the self-evolution of the probability density function and so, the percolation evolution. The link between the mathematics and physics issue (i.e. the non-reactive impregnation phenomenon) is done thanks to the initialization of the algorithm by an analytic solution of the homogeneous diffusion equation. In addition, the characterizations of the input parameters were also pointed out, related to the type of the considered capillary pressure curve, following van Genuchten's model or Brook & Corey's model or others. Furthermore, the convolution operator allowed to ensure the spatial continuity of the wetting fluid and to respect the boundary conditions of the treated problem.

In chapter III, the visual comparisons confirmed that the results of the SGP method agree well with the "reference" results obtained numerically by FEM and experimentally. To quantify these comparisons, we employed a Root Means Square Error (RMSE) and a Mass Balance Error (MBE) to calculate "quantitative" errors (i.e. between an SGP value and a "reference" value), the small error levels. The results therefore confirm the reliability of the SGP. Regarding the goals, first, the table of the CPU time both of the SGP method and the F.E. method were fully reported, showing that the computational duration of the SGP is reduced at maximum a factor 12 in the test 5; second, there are no possibilities of the issue of spurious oscillations during the transient stage because an analytical solution is used locally.

92/103

In conclusion, the SGP method algorithm fulfills the objectives of this work: (i) the solutions of the SGP method are free of spurious oscillations and (ii) the computational time is reduced.

IV.2 Outlook

To conclude this work, we recap here some open problems to be addressed in the future.

First, there is no direct link between the timestep of the SGP algorithm we call "percolation timestep" and real time. Currently, we have employed Washburn's equation to analytically associate each percolation timestep with real time, which has been validated by fitting with the "reference" results, but only after the simulation has been done (or otherwise, it will not respect the physics point of view). Therefore, an arising question is that how determine a robust stepwise relation between percolation and real time without this post identification? This will make us more confident about the real time of each a capillary pressure profile.

Second, there is no explanation for the physical meaning of parameter . We actually observe that the unity ⁄ of is the same than the intrinsic permeability, even the same for their values. We therefore assume that the meaning of parameter in physics is the intrinsic permeability, which is essential to be explained in our next researches.

Lastly, the first results considering non-reactive impregnation are promising. However, some details have to be modified to be able to first take into account the non-reactive impregnation in the 2D and 3D cases. In general, the strategy of the 1D SGP algorithm is the same for the 2D/3D one. In particular, its initialization should be the initial analytic solution of the homogenous equation, but with two/three standard deviations. The output of the algorithm will be the evolution of the local total pressure considered to be the self-evolution of these standard deviations. In our perspective, the kinetic of the evolution should be integrated with the physics such as capillary rise, rate of the liquid rise, etc. However, for this case, we are dealing with two problems: (1) solutions of the algorithm will strongly depend on the geometry of the treated problem which can lead to the time consuming; (2) the "reference" (i.e. experimental and F.E. numerical) results should be first implemented and then validated successfully, which will permit us to qualify the reliability of the SGP algorithm.

Second, the probability of chemical reactions and their effects will be worked on. For this stage, it is necessary to work on the properties of the chemical reactions. Then, the mass balance of the chemical elements will be studied, providing the diffusion equation of the concentration of the chemical reactions. It seems to allow that the SGP algorithm is applied to numerically solve this equation and, of course, there will require a lot of modifications to adapt. Today, we work on these modifications. In the next future, the objective is to implement this method in a FEM code to be able to compute a thermo-chemo-mechanical coupling without the necessity of solving the fully coupled PDE of the reactive mass balance. The CPU cost will be reduced and it may become possible to consider many more chemical reactions simultaneously.

93/103

List of figures

94/103

Fig. 26: The capillary pressure and the hydrostatic pressure allow to calculate the maximal "height" (i.e. the height at the equilibrium state), denoted by ℎ ( , ), and the "height" during the process with respect to the liquid at a given saturation, respectively.

95/103

List of tables

Fig. 1 :

 1 Fig. 1: Microscopic fluids occupation in a porous medium. Wetting phase view point.

Fig. 2 :

 2 Fig. 2: Microscopic fluids occupation in a pore assimilated to a smooth vertical tube.

Fig. 4 :

 4 Fig. 4: Capillary rise in Richards' conditions.

Fig. 5 :

 5 Fig. 5: The steady-state permeability measurement method (de Bilbao et al., 2018).

Fig. 6 :

 6 Fig. 6: The unity and the cubic function for the relative permeability.

Fig. 7 :

 7 Fig. 7: Relative function for various values of ranging from 0 to 1.

Fig. 8 :

 8 Fig. 8: Typical forms of the capillary pressure curve regarding (a) Brooks and Corey's model (1964) and (b) van Genuchten's model (1980) (Szymkiewicz, 2013).

Fig. 11 :

 11 Fig. 11: Graphical description of the NWP reservoir. The liquid and skeletons are colored in red and black, respectively[START_REF] Hilpert | Pore-morphology-based simulation of drainage in totally wetting porous media[END_REF].

Fig. 12 :

 12 Fig. 12: Schematic of the discretization in time. Considering the interval time , ], > , it is divided into time-step , , , ⋯ , ⋯ where designates the time increment.

  of the saturation level.

  Fig. 14: A capillary pressure profile is obtained through: (a) two different types of nodal element, (b) the type of 4-nodes element, and (c) the type of 8-nodes element. The chosen time-step does (the red curve) and does not (the blue curve) in respect of the critical time increment, respectively.

Fig. 16 :

 16 Fig. 16: On the left-hand side: numerical simulation of the crucible test is implemented in Abaqus/ Standard 2Dasymmetric (v6.14, 2017). On the right-hand side: the growth of town is modelized using the gradient percolation theory.

Fig. 18 :

 18 Fig. 18: An example of the gradient percolation model, where = ⁄ , the sitecoordinate of two points ( , ) with the second coordinate [ , ], and = on the triangular lattice[START_REF] Nolin | A Modified van Genuchten-Mualem Model of Hydraulic Conductivity in Korean Residual Soils[END_REF].

Fig. 19 :

 19 Fig. 19: A simulation result by the gradient percolation theory for different values of .

Fig. 20 :

 20 Fig. 20: An illustration of a multistate gradient percolation model at fixed values of and , with different values of , , .

Fig. 21 :

 21 Fig. 21: A local saturation (on the right-hand side) is derived from the gradient percolation model (on the left-hand side) using the convolution operator.

Fig. 22 :

 22 Fig. 22: An illustration of the continuous gradient percolation model.

Fig. 23 :

 23 Fig. 23: A simplified representation of Eq. (64).

Fig. 26 :

 26 Fig. 26: Evolution of the capillary pressure profile. On the left-hand side: it is at the 2000 th timestep; On the right-hand side: it is at the 10000 th timestep.

Fig. 27 :

 27 Fig. 27: The evolution of the standard deviation with ∆ = ( ℎ -⁄).

  liquid can flow out of the boundary surface.

Fig. 28 :

 28 Fig. 28: Outside pressure is forced at the top surface (the red line) of the sample.

  } to define the outflow or no flow boundary condition. Herein, boundary condition { : ⃗ • ⃗⃗ = } means that the flux is equal to zero in order to describe the case of no drain out of the boundary surface.

Fig. 30 :

 30 Fig. 30: The appropriate boundary conditions of the problem of interest.

  and , to describe that liquid can flow out of them, we consider a Neumann boundary condition i.e. and ⊂ . In fact, the drained boundary condition is imposed on and a constant air pressure boundary condition is imposed on .

Fig. 31 :

 31 Fig. 31: Boundary conditions are imposed on each boundary surface zone of a porous medium.

Fig. 33 :

 33 Fig. 33: Each capillary pressure profile is assimilated to a probability density function, where* is a fixed value of saturation and is distance from a capillary pressure profile to * .

Fig. 34 :

 34 Fig. 34: The simulated mass gain curve exported from Abaqus (v6.14, 2017).

Fig. 35 :

 35 Fig. 35: Comparisons of the (whole) mass gain curves. The given mass gain curve (dashed blue line) is compared with such curves of the SGP method associated with four different values of .

Fig. 36 :

 36 Fig. 36: Comparisons of the capillary pressure curves. The experimental CPC (blue points) is compared with such curves of the SGP method associated with three different values of * .

Fig. 37 :

 37 Fig. 37: Comparisons between the capillary pressure profiles of the SGP model (red lines) versus such curves of the F.E. model (blue points).

Fig. 40 :

 40 Fig. 40: Scheme of the device setup of the impregnation test (Coffignier, 2013).

Fig. 41 :

 41 Fig. 41: The analysis process of the F.E. Hydrus-1D (v4.16) package.2) Post-processing. The output of the simulation, as selected previously, are the capillary pressure profile at each time increment and the CPU time cost of the simulation.

Fig. 42 :

 42 Fig. 42: Numerical simulation for tests 4 and 5 has been implemented in the F.E. code Abaqus (axisymmetric, v6.14) in the quasi one-dimensional case.

Fig. 1 :

 1 Fig. 1: Microscopic fluids occupation in a porous medium. Wetting phase view point. ........ Fig. 2: Microscopic fluids occupation in a pore assimilated to a smooth vertical tube. .......... Fig. 3: Example of capillary pressure function, after (Albers, 2015). ..................................... Fig. 4: Capillary rise in Richards' conditions. ......................................................................... Fig. 5: The steady-state permeability measurement method (de Bilbao et al., 2018). ............ Fig. 6: The unity and the cubic function for the relative permeability. ................................... Fig. 7: Relative function for various values of ranging from 0 to 1. ................................... Fig. 8: Typical forms of the capillary pressure curve regarding (a) Brook-Corey's model (1964) and (b) van Genuchten's model (1980) (Szymkiewicz, 2013). ............................. Fig. 9: Capillary pressure curve S( ) is assumed as mirror image of Boltzmann transformation function S( ) .......................................................................................... Fig. 10: The hemisphere approach stands for the interfacial area between the "non-wetting" fluid and the "wetting" liquid, where a spherical presents the non-wetting fluid ( = ) -the wetting liquid ( = ) interface and is the inside radius of the hemisphere. .. Fig. 11: Graphical description of the NWP reservoir. The liquid and skeletons are colored in red and black, respectively (Hilpert and Miller, 2001). .................................................... Fig. 12: Schematic of the discretization in time. Considering the interval time , ], > , it is divided into time-step , , , ⋯ , ⋯ where designates the time increment. ........................................................................................................................................... Fig. 13: An example of a porous material with the specific dimensions and the mesh of elements of 4 and 8 nodes on the right and left side, respectively. ................................... Fig. 14: A capillary pressure profile is obtained through: (a) two different types of nodal element, (b) the type of 4-nodes element, and (c) the type of 8-nodes element. The chosen time-step does (the red curve) and does not (the blue curve) in respect of the critical time increment, respectively. .................................................................................................... Fig. 15: The schema of numerical methods about solving Richards' equation numerically. . Fig. 16: On the left-hand side: numerical simulation of the crucible test is implemented in Abaqus/ Standard 2D asymmetric (v6.14, 2017). On the right-hand side: the growth of town is modelized using the gradient percolation theory. ................................................ Fig. 17: Two typical types of the traditional percolation model (Grimmett, 1999)................. Fig. 18: An example of the gradient percolation model, where = ⁄ , the sitecoordinate of two points ( , ) with the second coordinate , , and = on the triangular lattice (Nolin, 2008).......................................................................................... Fig. 19: A simulation result by the gradient percolation theory for different values of . ...... Fig. 20: An illustration of a multistate gradient percolation model at fixed values of and , with different values of , , . ........................................................................................ Fig. 21: A local saturation (on the right-hand side) is derived from the gradient percolation model (on the left-hand side) using the convolution operator. ......................................... Fig. 22: An illustration of the continuous gradient percolation model. ................................... Fig. 23: A simplified representation of Eq. (64). .................................................................... Fig. 24: A simplified graphic representation of the link between (constant) standard deviation and standard deviation function ( ). ......................................................................

  ........................................................................................................................................... Fig. 27: Evolution of the capillary pressure profile. On the left-hand side: it is at the 2000 th timestep; On the right-hand side: it is at the 10000 th timestep. ......................................... Fig. 28: The evolution of the standard deviation with ∆ = ( ℎ -⁄ ). ........................ Fig. 29: Outside pressure is forced at the top surface (the red line) of the sample. ....... Fig. 30: The mass gain curves. (a) The grey curve (where = reference pressure of the CPC) is optimized to fit with experimental data (de Bilbao et al., 2016), the blue and the red curves are obtained through F.E.M by imposing = and = on the top surface, respectively. (b) The CPC with the reference pressure . ........................................................................................................................................... Fig. 31: The appropriate boundary conditions of the problem of interest. .............................. Fig. 32: Boundary conditions are imposed on each boundary surface zone of a porous medium. ............................................................................................................................ Fig. 33: The flow chart of the SGP algorithm. ........................................................................ Fig. 34: Each capillary pressure profile is assimilated to a probability density function, where * is a fixed value of saturation and is distance from a capillary pressure profile to * .67 Fig. 35: The simulated mass gain curve exported from Abaqus (v6.14, 2017). ...................... Fig. 36: Comparisons of the (whole) mass gain curves. The given mass gain curve (dashed blue line) is compared with such curves of the SGP method associated with four different values of . ....................................................................................................................... Fig. 37: Comparisons of the capillary pressure curves. The experimental CPC (blue points) is compared with such curves of the SGP method associated with three different values of * . ...................................................................................................................................... Fig. 38: Comparisons between the capillary pressure profiles of the SGP model (red lines) versus such curves of the F.E. model (blue points). ......................................................... Fig. 39: Experimental setup for tests 2 and 3 in order to measure both horizontal and vertical water infiltration processes (Davidson et al., 1963).......................................................... Fig. 40: a-The capillary rising test (de Bilbao et al., 2016; De Bilbao et al., 2015); b-The principle of the mass gain curve. ...................................................................................... Fig. 41: Scheme of the device setup of the impregnation test (Coffignier, 2013). ................. Fig. 42: The analysis process of the F.E. Hydrus-1D (v4.16) package. .................................. Fig. 43: Numerical simulation for tests 4 and 5 has been implemented in the F.E. code Abaqus (axisymmetric, v6.14) in the quasi one-dimensional case. .................................. Fig. 44: The whole calculated domain of the cylindrical porous sample , is divided into many small segments: {[ , ], [ , ], ⋯ , [ -, ]}.. ................................................ Fig. 45: The measured and simulated mass gain curves: (a) for the test 4 and (b) for test 5. .. Fig. 46: The evolution of measured and simulated capillary pressure profiles at different times: (a) for test 1, (b) for test 2, and (c) for test 3. ......................................................... Fig. 47: Evolution of the simulated capillary pressure profiles: (a) for test 4 and (b) for test 5. ........................................................................................................................................... Fig. 48: The measured and simulated capillary pressure profiles at the tenth time -step obtained through the SGP method with different small mesh sizes. ................................

  effect of the parameters to the multistate GP model. ......................................... Table 4: Initial values for the impregnation test worked on. ................................................... Table 5: The chosen values of and . ..................................................................................

Table 1 :

 1 Several phenomenological models of the capillary pressure curve.

	Model	Closed-form analytical equation	Parameters
	(Gardner, 1958)	= [( +	)	(			)]	⁄	+	,
	(Brooks and Corey, 1964)	= (		⁄ )	⁄				,
	(van Genuchten,	= [ +	|	| ]	--⁄	,
	1980)
	(Kosugi, 1996)	= . × erfc [	(	⁄ √	.	)	]	ℎ . ,
	(Schaap and
	van
	Genuchten,
	2006)

Table 2 :

 2 The formulations of the critical time increment in the respect of the type of nodal element of Richards' equation both in the case of 1D and 2D.

	1D case		2D case
	2 nodes	3 nodes	4 nodes	8 nodes

Table 4 :

 4 Initial values for the impregnation test worked on.

Table 6 :

 6 The optimized values of , , , and * for the test of interest.

Table 7 :

 7 Materials and liquid for the impregnation tests.

	Test	Porous sample	Liquid
		Marine sand	Water
		Columbia silt loam	Water
		Hesperia sandy loam	Water
		Degussit AL25	Glycerin
		Alumina 99%
		Low Cement	Engine
		Content Castable	Oil
		Alumina 99%
		(Alumina LCC)

Table 8 :

 8 Material properties for the F.E. simulations of tests 4 and 5

	Property	Value		Unit
		Test 4	Test 5
	Mass density	3870	3870	• -
	Specific weight of the liquid	12400	8742	• -
	Viscosity of the fluid	1.02	0.355	•
	Mass density of the fluid	1260	892	• -
	Intrinsic permeability	. × -	. × -
	Porosity	.	0.18	dimensionless

Table 11 :

 11 Computational time of the SGP model and the F.E model.

	Test	CPU time (seconds)
	SGP model	F.E. model
	0.716		1.72
	0.58		0.69
	0.33		0.49
	6.52		25.5
	1.1		12.1

Table 1 :

 1 Several phenomenological models of the capillary pressure curve. .........................

Table 2 :

 2 The formulations of the critical time increment in the respect of the type of nodal element of Richards' equation both in the case of 1D and 2D. ........................................

Table 3 :

 3 The

Table 6 :

 6 The optimized values of , , , and * for the test of interest. ............................

Table 7 :

 7 Materials and liquid for the impregnation tests. ........................................................

Table 8 :

 8 Material properties for the F.E. simulations of tests 4 and 5 .....................................

Table 9 :

 9 Values of input data used for the SGP algorithm. .....................................................

Table 10 :

 10 Comparisons the RMSE and MBE of the SGP, FE models, and experimental data for tests 1 to 5. ...................................................................................................................

Table 11 :

 11 Computational time of the SGP model and the F.E model. ....................................

	96/103

Fourier transformation is first applied to Eq. ( 54), to obtain:

The local saturation is then obtained by the application of the inverse Fourier transformation to Eq. ( 55) as follows:

In this section, the issue of continuity of the uncorrelated gradient percolation is tackled. Hence, Fourier transformation and the inverse Fourier transformation speed up the CPU time cost.

II.2.3 Problem 3: Self-organized Gradient Percolation (SGP) model

For this problem, the definition of a self-evolution process is required. We do that in the two following steps:

First, at the initial timestep, the initialization , is taken to be the analytical solution of the homogeneous diffusion equation [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF]. (ii) Second, in the next time steps this analytical solution is no longer available to describe profiles , . We therefore presume that the evolution of the profile is a self-evolution process of its initialization step.

Hence, time-depending a gradient percolation is constructed in a self-organizing way. The concept will be called Self-organized Gradient Percolation (SGP) model. The details of this model will be presented in Section II.3.

II.3 Self-organized Gradient Percolation

(SGP) algorithm

II.3.1 Richards' equation

The general equation governing matter and momentum balance is also the continuity equation [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF]. In this study, the governing equation for the impregnation process in the bounded domain [ , ] and on the time interval , ], > , will be:

The diffusion coefficient is to determine. To represent the transient fluid-transport problem, the traditional diffusion equation (Eq. ( 57)) will be identified with Richards' equation (Eq. ( 26)) at the initial timestep [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF][START_REF] Nguyen | Selforganized gradient percolation method for numerical simulation of impregnation in porous media[END_REF]. The diffusion coefficient is then given by: = (58)

During the impregnation process, porous sample is first analyzed at the initial timestep: There is the same saturation level (also known as residual saturation) for the entire porous sample, except at least at one boundary. It is thus found that the initial solution of Richards' equation is that of the homogeneous diffusion equation, which means that Eq. ( 57) is considered in the case of a constant diffusion coefficient.

A simplified 1D analysis of the homogeneous diffusion equation is briefly presented here in order to introduce its general analytic solution (it is illustrated in more detail in [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF]Tropea et al., 2007).

Using the spatial Fourier transformation ℱ in a first attempt to solve the homogeneous diffusion equation we get

One can then solve the first-order differential equation to obtain

where ℱ ( , ) is the Fourier-transformation image at the time zero.

where = .

For the one-dimensional case, the capillary pressure profiles are obtained by the following convolution form:

where:

53/103

• ( , ) is a probability density function corresponding to the Gaussian distribution [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF][START_REF] Nguyen | Selforganized gradient percolation method for numerical simulation of impregnation in porous media[END_REF]:

at time ; • is the standard deviation of the distribution to be determined later;

•

, stands for the smoothing point-spread function ( , ), constant on .

Let us stress out that the 1D convolution is taken over all = , of same "height" .

To be physically meaningful, the use of Eq. ( 62) to describe the capillary pressure profile at initial timestep requires the density function satisfies the following:

(i) As a consequence, it is proposed to modify Eq. ( 63) both in terms of the variance and the exponent in the following way:

where characterizes the type of a distribution. In particular = or = correspond to the probability density functions for Gaussian or Laplace distribution by means of statistics, respectively.

By injecting Eq. (64) into Eq. ( 62), we see that the initialization of the SGP algorithm is a convolution between the density function and the smoothing point-spread function.

Obviously, there are two points of view in the description of the capillary pressure profile:

• The first one is the mathematical point of view in which the capillary pressure profile regarding Eq. ( 64) is the probability density function of a statistics distribution (lefthand side of Fig. 23); • The second one is the physics point of view where the capillary pressure profile is usually represented by capillary pressure/height along vertical direction and saturation along horizontal direction (right-hand side of Fig. 23).

For a given saturation, the capillary pressure depends only on the "height", i.e. the vertical distance to the origin.

Hence, the physics interpretation of the initialization (from the mathematics point of view) is given by the graph on Fig. 23. Considering the usual capillary pressure curve on the left Considering initial standard deviation , the standard deviation function at the first timestep (i.e. = ) can then be expressed as:

On the other hand, a simplified graphic representation of the capillary pressure profiles will be produced as the following: In Fig. 24, we display the relation between standard deviation and standard deviation ( ) which is observed as the trapezium:

Eq. ( 66) is then rewritten as following

where (

) is an unknown function of the capillary pressure.

56/103

To determine ( ), the Welge's method and the Buckley-Leverett theory [START_REF] Christiansen | Evaluation of methods for measuring relative permeability of anhydride from the Salado Formation: Sensitivity analysis and data reduction[END_REF] are employed. The local saturation at a small increment , denoted by + , can be expressed by

where is the liquid volume of the local pore space and can be also represented by a proportionality constant [START_REF] Christiansen | Evaluation of methods for measuring relative permeability of anhydride from the Salado Formation: Sensitivity analysis and data reduction[END_REF] , and( ) is a function of relative permeability .

The formulation of the relative permeability is determined by the Buckley-Leverett theory, as following [START_REF] Ef. Apostolos | Fundamentals of Fluid Flow in Porous Media[END_REF]:

where is the viscosity of the "wetting" fluid.

The unified expression composed of Eq. ( 68) and Eq. ( 69) then gives

which yields the general standard deviation function, i.e. = ,

Then, the capillary pressure profile with its evolution in time will be deduced from Eq. ( 62) combined with Eq. ( 64) and the evolution of the standard deviation of a distribution as Eq. (72).

To close the expression, the formulation of , is assumed to be proportional to the standard deviation function:

where is a positive constant.

Therefore, the evolution of the density function corresponds to the evolution of the standard deviation of the probability density function. To ensure the fully understandable link with the initial physics problem, it is thus necessary to point out the link between this evolution and the real time of the problem.

57/103

II.3.2 Time evaluation

The present method considers the physics at the local microscale of the pores. Then, locally, the pore is assimilated to a capillary. To make the bridge between the self-organized gradient percolation and the local physics, it is proposed to reproduce the evolution of the height of the liquid with timestep (i.e. the derivative of the height of the liquid over time , denoted by ⁄ ) by the evolution of the standard deviation.

Since Poiseuille's equation [START_REF] Extrand | Forces, pressures and energies associated with liquid rising in nonuniform capillary tubes[END_REF] also derives the time evolution of the capillary rising, it appears reasonable to link the evolution of the standard deviation with the Poiseuille's equation. This idea leads to:

where and ℎ are the total pressure and the hydrostatic pressure at timestep , respectively.

In the case of the vertical capillary tube, Eq. ( 74) can be rewritten in the following way:

Hence, Eq. ( 75) induces that the evolution is equal to zero when there is no difference between capillary pressure and hydrostatic pressure (i.e., ℎ = ) with respect to the physics point of view.

To do that in the numerical practice, Eq.( 75) should be rewritten as the following form:

where = ⁄ depends on capillary diameter and the liquid properties.

To reproduce the evolution of standard deviation ( , ), capillary pressure and hydrostatics pressure ℎ need to be numerically calculated as the following:

The first is related to the capillary pressure curve (at the steady state) to identify ℎ ( , ) which is given as:

The second is related to the hydrostatic pressure (during the impregnation process) to determine the standard deviation such that = , ℎ .

Eq. ( 94) points out that there are two different cases: the liquid flow occurs at the imposed surface when its pore pressure is greater than the "reference" pore pressure value, otherwise, no flow of the liquid occurred.

When the value of the seepage coefficient tends to the infinite, the corresponding value tends to the reference pore pressure value ∞ for Eq. ( 93) and to zero for Eq. ( 94). It has been established that the value of the seepage coefficient number is large enough to be considered as an infinite number when it satisfies (Anon, 2011):

where and are the permeability of the porous material and the characteristic length of the underlying element, respectively.

To determine the suitable boundary condition on each surface of the cylindrical sample, we consider each surface separately in respect to the physics point of view as the following:

1) For the bottom surface, its pore pressure is imposed to be equal to the atmospheric pressure. The goal is that the liquid flow can penetrate into the sample; 2) For the vertical surface, the velocity of the flow on this surface should be equal to zero.

The drainage only flow type label (denoted "QD") is thus imposed on;

3) For the top surface, the air pressure/ the drainage only flow type label (denoted "Q") is imposed on.

III.3 SGP model of the experimental tests

The values of the input data needed in the SGP algorithm for the five tests are reported in Table 9. The chosen value permits the same capillary pressure curve than for FEM for each tests.

In order to evaluate the performance of the SGP method its numerical results are compared with the solutions obtained through the FEM and experimental measurements. In particular, the reference solutions considered in the framework of the thesis are the capillary pressure profiles.

83/103

Then considering Eq. ( 97) on the calculated domain (Fig. 43) yields the transformed form as the following:

Thanks to Reimann approximation, we can rewrite the above sum (Eq. ( 98)) into the corresponding integral expression as follows:

Note that, it has been experimentally observed in [START_REF] Coffignier | Mesures et simulation de l'imprégnation non réactive d'une alumie poreuse par un liquide[END_REF] that the porosity of our samples is almost homogeneous. So, it leads to:

III.4 Comparisons between numerical and experimental results

The comparisons of numerical results from the SGP and F.E. M. with the experimental results are presented as:

For tests 1 to 3, capillary pressure profiles obtained by the SGP and F.E.M. are compared with the measurements, respectively, as displayed in Fig.

(ii)

For tests 4 and 5, it is to note that there is no experimental data for the capillary pressure profiles. Hence, the simulated mass gain curve from the SGP method is compared with the measured mass gain curve (de Bilbao et al., 2016), as displayed in Fig. 44.

In order to quantify the gap between the different methods and between numerical results and experimental ones, it is proposed to defined discrepancy measure based on the least square root error. For tests 1 to 3, The Root Mean Square Error (RMSE) of the capillary pressure profiles is given by:

84/103

where , and , are the capillary pressure profile at time , obtained by measurements and the considered numerical method (i.e. the SGP/F.E. method), respectively; denotes the number of the experimental data points.

For tests 4 and 5, the Mass Balance Error (MBE) of the mass gain curve is given by:

where and represent the total mass gain at initial and given time , respectively; the subscripts "Exp" and "SGP" denote the mass gain curve obtained experimentally and through the SGP method, respectively.

The RMSE and the MBE, associated with the evaluation of the performance of the SGP method for the five tests, are listed in the Table 10.

86/103

In summary, Fig. 44 and Fig. 45 show the capillary pressure profiles at different times and the mass gain curves, respectively. It can be seen that numerical results both from the SGP and F.E. results fit quite well with the experimental results.

III.5 Comparisons between SGP and F.E.M results

The comparisons of the capillary pressure profiles at different times between the SGP versus F.E. methods for tests 1 and 3 can be observed in Fig. 45 For tests 1 to 5, the RMSE is specified in the comparisons of the capillary pressure profiles obtained through the F.E. versus SGP methods, where the subscripts "Exp" and "SGP" of the Eq. ( 101) are transferred into "F.E." and "SGP", respectively.

, is a capillary pressure profile obtained by using the F.E.M at time . The RMSE will be used to evaluate the performance of the SGP method for each test, as reported in Table 10.

The FEM results for tests 4 and 5 are plotted for comparisons with the SGP results in Fig. 46. It is interesting to see that the numerical solutions from the SGP method are very similar to that from the F.E.M for test 4. For the test 5, the slight difference between the SGP and F.E. solutions, is due to the shape of the Laplace distribution for the capillary pressure profile (i.e. reference parameter m for the probability density function equal to unity). Indeed, the assumption of the parameter m equal to unity is assumed to reproduce the CPC predicted from

87/103

Brook and Corey's model. This choice is debatable, perhaps that a value slightly different from 1 will give a better result.

III.6 Results and discussions

The interest of the SGP algorithm will be study by analyzing the results of the Root Means Square Error (RMSE) between the simulated capillary pressure profiles (by using the SGP method for tests 1 to 3 and the F.E.M for tests 4 and 5) and the measured ones, as well as the results of the Mass Balance Error (MBE) between the simulated mass gain curve by using the SGP method and the measured one for tests 4 and 5.

The capillary pressure profiles from experimental data don't correspond to the same time for each test. Indeed, measurement have been done at 15 min, 60 min and 300 min for test 1 (Fig. 45 (a)); at 88 min, 344 min and 740 min for test 2 (Fig. 45 (b)); at 158 min, 620 min and 1476 min for test 3 (Fig. 45 (c)); at 1 second, 1000 seconds, 7000 seconds and 30000 seconds for test 4 (Fig. 46 (a)); at 100 seconds, 5000 seconds, 120000 seconds and 300000 seconds (Fig. 46 (b)). Then, time increment for FEM have been adapted and interpolation between time increment for SGP have been done. It result that:

88/103 (i)

For tests 1 to 3, the RMSE between the SGP and experimental results, as well as the SGP and FEM present small error levels (see Table 10): ranges from 0.097% to 1.02% and from 0.0126% to 0.108% for the test 1; from 0.0422% to 0.093% and from 0.00179% to 0.1233% for the test 2; from 0.022% to 0.0346% and from 0.0387% to 0.0902% for the test 3, respectively. (ii)

For tests 4 and 5, the RSME between the SGP and F.E. results and the MBE between the SGP and experimental results present small error levels (see Table 10): it ranges from 0.00127% to 0.027% and from 0.1226% to 1.33% for the test 4; from 0.164% to 0.362% and from 0.0207% to 2.45% for test 5, respectively. Considering the above results, it can be concluded that the SGP algorithm produce good results for the simulation of one-dimensional impregnation.

The CPU time for these five tests are reported in Table 11 for SGP and F.E.M methods, considering a 3.07 GHz, x64-based processor with 16 GB of RAM. The computational time of the SGP algorithm is significantly lower than that of the F.E.M. The efficiency of the SGP algorithm is highest for the test 1 (by a factor higher than 2) and test 4 (by a factor of 4). This variation in the ratio between the CPU time of both methods can be explained by the difference in the evolution of timestep between both methods.

First consider the F.E.M for the (1D) Richards' equation and put = Ψ

. Its discretization is then given by [START_REF] Wu | A finite-element algorithm for modeling variably saturated flows[END_REF]:

where , , and are time, time increment, and iteration level, respectively; + , = ( + , ).

We see that the value of time increment is chosen in respect of Richards' equation [START_REF] Szymkiewicz | Modelling Water Flow in Unsaturated Porous Media, GeoPlanet: Earth and Planetary Sciences[END_REF]. The Eq. ( 103) is the system of the discrete nonlinear equations which must be transferred to discrete linear equations by using an iteration method. More precisely, the "exact" left-hand side term of Eq. ( 103) is replaced by an approximate term by means of assumptions. The difference between these two terms is called "the mass-balance error". This error is demanded to be "small enough" to ensure the accuracy of the numerical solutions. As a consequence, the time increment is a small value in respect of . Herein, is the critical time increment [START_REF] Vermeer | An accuracy condition for consolidation by finite elements[END_REF] and is the maximum value of time increment which is deduced by the maximum wetting liquid pore pressure change allowed in an increment (Anon, 2011). It means that the time increments are restricted to guarantee that this value is not exceed at any node.

Second, consider the SGP algorithm. It is to notice that the evolution of the solution is the evolution of the standard deviation, consistent with the analytic solution at any time. Hence, time increment is independent of the classical limitations for the numerical integration of the Richards' equation. That is why the variation in the ratio between the CPU time of both methods is not always the same.

Anh Khoa NGUYEN METHODE DE PERCOLATION EN GRADIENT AUTO -ORGANISEE POUR UNE SIMULATION NUMÉRIQUE D'IMPREGNATION DANS LE MILIEUX POREUX

Résumé :

La modélisation numérique du phénomène d'imprégnation dans les milieux poreux a une très large gamme d'applications. Le développement d'un modèle doit tenir compte de la complexité multiphysique du problème afin de rester fidèle à la physique du phénomène. Dans le cadre de cette thèse, nous proposons un modèle numérique pour la prédiction de l'évolution de l'imprégnation non-réactive quasi-unidimensionnelle et en particulier de sa courbe de succion capillaire. Ce modèle est basé sur une approche purement probabiliste, avec des paramètres liés aux lois physiques régissant le procédé d'imprégnation, passant outre la résolution des équations aux dérivées partielles qui sont à l'origine des difficultés numériques actuelles. Les enjeux de cette méthode, appelée Self-Organized Gradient Percolation, sont doubles : réduction significative des temps de calcul et résultats fiables exempts de toutes oscillations parasites. Les points clés de la méthode résident dans la simulation de l'évolution du front d'imprégnation selon une fonction de densité de probabilité et dans l'identification de la relation entre les paramètres du modèle et les propriétés matériaux. Les premiers résultats, comparés à des résultats numériques provenant d'autres méthodes et à des résultats expérimentaux, sont prometteurs. Une évolution de ce modèle dans la deuxième dimension est proposée dans les perspectives.

Mots clés : imprégnation, succion capillaire, percolation en gradient, simulation.

SELF-ORGANIZED GRADIENT PERCOLATION METHOD FOR NUMERICAL SIMULATION OF IMPREGNATION IN POROUS MEDIA

Abstract:

The modeling of the impregnation of porous media with liquid is a research topic with a wide scope of applications. The numerical models that are faithful to the laws of physics require a multiphysics approach taking into account the complexity of this multiphysics coupling. The aim of the thesis is to propose a numerical model which is able to reproduce the evolution of the capillary pressure profile for the (quasi) one-dimensional non-reactive unsaturated impregnation. This model is based on a pure probabilistic approach, where its parameters are related to the laws of the physics governing the impregnation process and avoiding the resolution of the partial differential classical equation which is the arising origin of the current numerical locks. This method, called Self-Organised Gradient Percolation, has two main objectives: drastically reducing CPU time and ensuring reliable results to be free from spurious oscillations. The key points of the method are: the assumption that the capillary pressure profile is considered as a probability density function and the strong coupling between the model parameters and the material properties. The first results, compared to numerical results from finite element method and experimental results, are promising. An extension of the model to the two-dimensional case is proposed in the outlook. Keywords: impregnation, capillary suction, gradient percolation, simulation.

Laboratoire de Mécanique

Gabriel Lamé (LaMé) Univ. Orléans, LaMé (EA7494, Univ. Orléans, Univ. Tours, INSA CVL), France