Ainagul Jumabekova
email: ainagul.jumabekova@univ-smb.fr

Julien Berger
email: berger.julien@univ-smb.fr

Denys Dutykh
email: denys.dutykh@univ-smb.fr

Hervé Le Meur
email: herve.lemeur@math.u-psud.fr

Aurélie Foucquier
email: aurelie.foucquier@gmail.com

Mickaël Pailha

Christophe Ménézo
email: christophe.menezo@univ-smb.fr

Mickael Pailha
email: mickael.pailha@univ-smb.frch

An efficient numerical model for liquid water uptake in porous material and its parameter estimation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Nowadays efficient energy consumption of buildings is one of the most challenging tasks. In case of the historical buildings, this problem has to deal not only with the energy performance but also with the building protection itself [START_REF] Falchi | Rising damp in historical buildings: A Venetian perspective[END_REF]32]. The study of moisture's amount in the building walls takes the largest part for energy reduction since it can raise the heat loss through the walls. The increase of moisture quantity may occur due to different reasons as wind-driven rain, rising damp or simply vapour diffusion [4].

Among all the potential sources of moisture, rising damp is a major issue since it is an important source of liquid water. It may strongly impact the building structure by modifying the mechanical behaviour of the wall. Moreover, the water from the ground brings dissolved salts. During the evaporation process, these salts precipitate and crystals appear that deteriorate the material [27,[START_REF] Scherer | Stress from crystallization of salt[END_REF]. Some examples of such deterioration can be found in [22] or [26].

Thus, one understands the importance of studying the physical phenomena occurring in rising damp. As mentioned by Franzoni in [START_REF] Franzoni | Rising damp removal from historical masonries: A still open challenge[END_REF], the main challenges related to the rising damp problem is to find accurate, fast and cheap methods to characterize the moisture content during the water uptake of a material or a wall. Several methods to characterize the amount of moisture already exist and can be divided in two main groups: (i) invasive (traditional like weighting-drying, chromatography) and (ii) non-invasive (dielectric, microwave) [17]. But usually, these approaches are very complex procedures requiring extensive technical support. Moreover, the obtained results do not describe an exact situation of the whole building. Thus, as underlined in [START_REF] Falchi | Rising damp in historical buildings: A Venetian perspective[END_REF], it is worth investigating other approaches by proposing numerical models to perform simulations and obtain quantitative results.

Several models were proposed in the literature for the simulation of liquid water uptake in a porous material. A first approach is to propose analytical solutions of the problem as for instance in [START_REF] Guimarães | Mathematical analysis of the evaporative process of a new technological treatment of rising damp in historic buildings[END_REF]18]. However, the main drawback is that these solutions work only for idealized conditions. As an example, in [18] a step-function is used for the diffusivity coefficient, neglecting its nonlinear variation with the fields. This assumption is not valid for the practical case study. Other approaches are based on standard numerical methods in [14, [START_REF] Häupl | Coupled heat air and moisture transfer in building structures[END_REF][START_REF] Peszynksa | Numerical methods for unsaturated flow with dynamic capillary pressure in heterogeneous porous media[END_REF][START_REF] Torres | The influence of the thickness of the walls and their properties on the treatment of rising damp in historic buildings[END_REF]. However, these approaches also have an important drawback. When using explicit Euler approach, the standard stability conditions must be respected implying very small values of discretization time grids due to the high nonlinearity of the problem. Thus these numerical models have a high computational cost. These works promote the use of implicit schemes for their unconditional stability properties. However, they require a large number of sub-iterations during the computation to treat nonlinearities, which also lead to an increase in the computational cost of the numerical models. It is of major importance to propose efficient numerical models to represent the physical phenomena of capillary adsorption in a porous material.

When proposing efficient numerical models, one important point is its reliability to represent the physical phenomena. For this, the numerical predictions need to be confronted to experimental observations. In the case of rising damp, there are several material properties involved in the definition of the physical model. To have an accurate numerical prediction, the moisture diffusivity and the hydraulic conductivity must be known precisely. As mentioned in [14], there is a lack of data in the literature for these material properties. Thus, when comparing the numerical predictions to experimental observation, the material properties may be estimated to calibrate the numerical model. This procedure requires solving the parameter estimation problem. The solution of such type of the inverse problem requires the computation of the so-called direct problem several times. Thus, since current numerical models have high computational cost, the estimation of material properties becomes an obstacle.

Therefore, the objective of this article is twofold. First, an innovative numerical model is proposed for capillary adsorption phenomena in a porous material. It is based on the Scharfetter-Gummel numerical scheme. This approach is particularly efficient for socalled advection-diffusion equations as highlighted from a mathematical point of view in [START_REF] Gosse | Viscous equations treated with L-splines and Steklov-Poincare operator in two dimensions[END_REF] and illustrated in [2, 3] for the case of heat and moisture transfer in building porous materials. In our work, the proposed numerical model is compared to the standard methods in the context of capillary adsorption phenomena. The issue is to prove its reliability, accuracy and smaller computational cost. Then, the numerical model is compared to experimental data to evaluate its accuracy to represent the physical phenomena. To calibrate the model, a parameter estimation problem is solved in a good time to determine the material properties of a brick using observed data of water uptake.

The paper is organized in the following way. In Section 2, the physical model is described. Following Section 3 presents the Scharfetter-Gummel scheme. Its accuracy and advantages are compared to classical methods on numerical case studies. Section 4 describes the experimental facility used to generate observation of water uptake. Section 5 aims at comparing the numerical predictions and the experimental observations, as well as estimating the material properties of the brick to calibrate the model.

Physical model

Liquid transfer

Let us present the physical model of liquid water uptake process in a porous material. The physical model is inspired by the one proposed by Philipp and De Vries in 1957 [25]. Since the problem deals with high moisture content, the exchange between vapor and liquid water is supposed negligible. Moreover, it is assumed that no chemical reaction occurs. Thus, the liquid water conservation equation is given by:

∂θ ∂t = -∇ • j ,
where θ m 3 /m 3 is the volume basis of liquid water content in the material and -→ j kg/(m 2 .s) is the liquid flow through the capillaries. The liquid flow is driven by diffusion, advection due to air pressure difference and gravity forces [START_REF] Hall | Water Transport in Brick, Stone and Concrete[END_REF]:

j = j d + j a + j g .
The diffusion liquid transfer flux is expressed as:

j d = D θ ∇θ + D T ∇T ,
where D θ m 2 /s is the liquid transport coefficient under a water content gradient and D T m 2 /(s.K) is the liquid transport coefficient associated to a temperature gradient. Both depend on water content θ . Since the experiments are conducted under isothermal conditions, the liquid diffusion under temperature gradient is assumed negligible. Thus, we have:

j d = D θ ∇θ .
The transfer is also driven by gravity forces described by the following expression:

j g = K i ,
where K m/s is the hydraulic conductivity depending on θ , and i is taken as pointing upward as illustrated in Figure 1. To take into account the motion due to the gradient pressure, an advection flux is introduced. It represents liquid water transfer induced by the filtration of dry air:

j a = Π θ v ,
where v m/s is the mass average velocity and Π the porosity of the material. Since the filtration of dry air only influences the water front at x = H , the mass average velocity is expressed as:

v = 1 - H L v 0 ,
where H is the water height in the brick, calculated as the total moisture content at the current instant t :

H ( t ) = ˆL 0 θ ( x, t ) θ sat dx ,
where θ sat -saturation moisture content.

Let us define coefficient a 0 and approximate it using Darcy's law as:

a 0 = Π v 0 ≈ Π k a µ ∆P L ,
where k a m 2 is the intrinsic permeability of the material, and µ Pa • s is the dynamic viscosity of the fluid. Finally, the 1-dimensional liquid transfer equation is expressed as:

∂θ ∂t = ∂ ∂x D θ ∂θ ∂x -a 0 ( 1 -H ) θ i -K i . (2.1)
In the physical model, Dirichlet boundary conditions are set since the surface convective transfer coefficients are unknown. At the boundary where liquid uptakes, the water content is set to saturation. At the top of the brick, the water content remains equal to the initial condition where no liquid is present in the brick. Thus, the boundary and initial conditions are defined as:

θ ( x = 0 , t ) =    0 , t = 0 θ sat , t > 0 , θ ( x = L , t ) = 0 , θ ( x , t = 0) = 0 .

Dimensionless formulation

For numerical analysis, it is very useful to work with the unitless formulation of the model. The first reason is the computational accuracy. Scaling the problem to the appropriate order, many terms will have the same order, so the effects of numerical errors are minimized when computing the residual. The second reason is that many useful relationships exist between dimensionless numbers, which show how specific things influence the whole system.

The following dimensionless variables for the moisture content, the time and space domains are defined:

u = θ θ sat , x ⋆ = x L , t ⋆ = t t ref ,
where θ sat -saturation moisture content, L -length of the brick, t ref -set to one hour.

Then all the material properties are modified considering a reference value:

d ⋆ ( u ) = D θ ( θ ) D ref , a ⋆ 0 = a 0 a ref , k ⋆ ( u ) = K( θ ) k ref .
The Fourier number measures the relative importance of the heat and mass transfers through the material:

Fo = t ref D ref L 2 .
The Péclet number demonstrates the relative importance of the advection against the diffusion transfer:

Pe = a ref L D ref .
The importance of the gravity forces with the respect to the diffusion is quantified by the Bond number:

Bo = k ref L θ sat D ref .
Finally, the following dimensionless equation of liquid transfer is obtained:

∂u ∂t ⋆ = Fo ∂ ∂x ⋆ d ⋆ ( u ) ∂u ∂x ⋆ -Pe a ⋆ 0 1 -H ⋆ u -Bo k ⋆ ( u ) , (2.2)
where the dimensionless water front height is given by:

H ⋆ = ˆ1 0 u ( x ⋆ , t ⋆ ) d x ⋆ . (2.3)
Further in article, for the sake of clarity the upper script ⋆ is omitted, and Eq. (2.2) transforms into:

∂u ∂t = Fo • ∂ ∂x d ( u ) • ∂u ∂x -Pe • a ( u ) • u -Bo • k ( u ) , (2.4)
with the following boundary and initial conditions:

u( x = 0 , t) =    0 , t = 0 1 , t > 0 , u( x = 1 , t) = 0 , u( x , t = 0) = 0 .
Next section will introduce how to solve the retrieved equation effectively.

Numerical method

The issue is now to propose an efficient numerical model to compute an accurate solution with a reduced computational cost. For this, the Scharfetter-Gummel numerical scheme for an advection-diffusion equation with gravity flux is introduced. Its mathematical properties are discussed. For the sake of simplicity, the following differential equation is considered:

∂u ∂t = ∂ ∂x d • ∂u ∂x -a • u -k (u) , t > 0 , x ∈ 0 , 1 , (3.1)
where a , d are constants and coefficient k (u) depends on u . Dirichlet boundary conditions are taken as

u ( 0, t ) = u L ( t ) , u ( 1, t ) = u R ( t ) .
A brief demonstration of the uniqueness of the solution of Eq. (3.1) is provided for the linearized equation in Appendix A. It can be remarked that the uniqueness of advectiondiffusion equation has been demonstrated in [5] considering a nonlinear advective term and a linear diffusive one with Dirichlet-type boundary conditions. Some demonstrations have been made for multi-dimensional systems of advection-diffusion equations also with Dirichlet-type boundary conditions. In [START_REF] Evans | Partial Differential Equations[END_REF], the problem considered is linear. In [START_REF] Kusano | On the first boundary problem for quasilinear systems of parabolic differential equations in non-cylindrical domains[END_REF], it is quasi-linear since the diffusion coefficients depend on space and time whereas the advective and source terms are nonlinear. Some proofs have been also produced for a system of diffusion equations with diffusion coefficients depending on the fields and Neumann-type boundary conditions in [START_REF] Araújo | Cross-Diffusion Systems for Image Processing: II. The Nonlinear Case[END_REF].

Let us discretize uniformly the space and time intervals, with the parameters ∆x and ∆t , respectively. The discrete values of function u (x, t) are defined as u n j ≡ u ( x j , t n ) , where j ∈ { 1, . . . , N } and n ∈ { 1, . . . , N t } .

The Scharfetter-Gummel numerical scheme

The flux J is denoted by the following expression:

J = d • ∂u ∂x -a • u -k (u) . For the cell C = x ∈ x j- 1 2 , x j+ 1 2
illustrated in Figure 2, the semi-discretization of the equation (3.1) results as:

d u j d t = 1 ∆x J n j+ 1 2 -J n j-1 2 . (3.2)
Within the Scharfetter-Gummel approach, the assumption is that the flux J is constant on the dual cell C ⋆ = x j , x j+1 . The following boundary value problem can be written:

J n j+ 1 2 = d • ∂u ∂x -a • u -k n j+ 1 2 , ∀x ∈ [ x j , x j+1 ] , (3.3a) u = u n j , x = x j , (3.3b) u = u n j+1 , x = x j+1 . (3.3c)
where k n j+ 1 2 is approximated by:

k n j+ 1 2 = k 1 2 u n j + u n j+1 .
It is important to note that problem (3.3) has two unknowns u ( x , t ) and J n j+ 1 2 with two constrains at x = x j and x = x j+1 . Therefore, one can obtain the exact solution of Eq. (3.3) as:

J n j+ 1 2 = -a • u n j+1 -u n j e a ∆x d 1 -e a ∆x d -k n j+ 1 2 . (3.4)
Applying expression (3.4) to Eq. (3.2), and using an Euler explicit approach, we obtain the expression to compute u n+1 j :

u n+1 j = u n j + ∆t ∆x • -a • u n j+1 + u n j-1 e a ∆x d 1 -e a ∆x d + a • 1 + e a ∆x d 1 -e a ∆x d • u n j -k n j+ 1 2 + k n j-1 2 .
In this work, several approaches are used for the temporal discretization. As presented above Euler explicit approach can be applied. Additionaly, the solution is computed with adaptive time step the Adams-Bashforth-Moulton algorithm, using Matlab™ function ODE113 [START_REF] Shampine | The MATLAB ODE Suite[END_REF].

Properties of Scharfetter-Gummel scheme

The main advantage of the Scharfetter-Gummel numerical scheme is its explicit formulation of the solution. Thus, the nonlinearities are handled without costly sub-iterations at each time steps as it is the case with implicit approaches. Another significant point, the scheme is well-balanced as well as asymptotically preserving [START_REF] Gosse | Viscous equations treated with L-splines and Steklov-Poincare operator in two dimensions[END_REF]. The Scharfetter-Gummel scheme has first order accuracy over the time and space O (∆x + ∆t) [11]. Moreover, the flux is also approximated to the first order O ( ∆x ) .

For Eq. (3.1), the Courant-Friedrichs-Lewy (CFL) stability condition for the Scharfetter-Gummel numerical scheme combined with an Euler explicit approach, is calculated with the expressions from [START_REF] Gosse | Viscous equations treated with L-splines and Steklov-Poincare operator in two dimensions[END_REF]:

∆t • max j a + dk du u x j+ 1 2 • tanh a + dk du u x j+ 1 2 • ∆x 2 d -1
∆x .

(3.5)

Figure 2. Stencil of the Scharfetter-Gummel numerical scheme

Through this condition, a nonlinear relation stands between ∆t and ∆x . But if large space discretization ∆x is taken, time and space grid become proportional to each other ∆t ≃ ∆x . Thus, the stability condition is relaxed compared to the classical central finite-difference scheme combined with Euler explicit approach ∆t ≃ ∆x 2 . As mentioned before, Eq (3.3) has two unknowns u ( x , t ) and J n j+ 1

2

. The solution J n j+ 1 2 is given by Eq. (3.3). The exact expression of solution u can be computed for x ∈ C ⋆ :

u n ( x ) = u n j e a ∆x d -u n j+1 e a ∆x d -1 + u n j+1 -u n j e a ∆x d -1 • e a • (x -x j ) d ,
x ∈ x j , x j+1 .

Extension for the nonlinear case

In case where the coefficients a, d and k of Eq. (3.1) all depend on u , we approximate them by using the frozen coefficients assumption on the dual cell x j , x j+1 .

Thus, the solution of the boundary-value problem (3.3) is written as:

J n j+ 1 2 = -a n j+ 1 2 • u n j+1 -u n j e θ n j+ 1 2 1 -e θ n j+ 1 2 -k n j+ 1 2
, where

θ n j+ 1 2 = a n j+ 1 2 • ∆x d n j+ 1 2
, and

a n j+ 1 2 = a 1 2 u n j + u n j+1 , d n j+ 1 2 = d 1 2 u n j + u n j+1 , k n j+ 1 2 = k 1 2 u n j + u n j+1 . (3.6)
For the nonlinear case the CFL stability condition for the Scharfetter-Gummel numerical scheme is is calculated with the expressions from [START_REF] Gosse | Viscous equations treated with L-splines and Steklov-Poincare operator in two dimensions[END_REF]:

∆t • max j d j max j p j+ 1 2 d j+ 1 2 tanh p j+ 1 2 • ∆x 2 d j+ 1 2 -1 ∆x ,
where

p j+ 1 2 = a j+ 1 2 + u j+ 1 2 • da du u x j+ 1 2 + dk du u x j+ 1 2
.

Comparing numerical models

In order to validate the numerical model, the error between the solution u num (x, t) and a reference one u ref (x, t) is evaluated as a function of x according to the formula:

ε 2 (x) ≡ 1 N t Nt j=1 u num x , t j -u ref x, t j 2 ,
where N t is the number of temporal steps. The local uniform error is defined as the maximum value of ε 2 (x) :

ε ∞ ≡ max x ∈ [ 0, L x ] ε 2 (x) .
The significant correct digits (scd) of the solution is calculated according to [START_REF] Söderlind | Adaptive time-stepping and computational stability[END_REF]:

scd(u) ≡ -log 10 u( x, t end ) -u ref ( x, t end ) u ref ( x, t end ) ∞ .
The reference function u ref x, t is given by a numerical solution of the differential equation based on the Chebyshev polynomial and adaptive spectral methods and obtained using the function pde23t from the Matlab™ open source package Chebfun [6].

Numerical validation

In this Section, the purpose is to perform numerical computations to validate the proposed model. Thus, only dimensionless quantities are considered with no necessary physical meanings. x 0 = 0.5

x 0 = 0.99

x 0 = 0.1 SG Ref (b)
Figure 3. Variation of the field u as a function of (a) space and (b) and time t .

Linear case

First, we study a case with constant functions d ( u ), a ( u ) and a linear function of k ( u ):

d ( u ) = d 0 = 0.05 , a ( u ) = a 0 = 0.02 , k ( u ) = k 1 u ,
with k 1 = 0.5 . It should be noted that the dimensionless numbers are set to unity in Eq. (2.4). The following boundary and initial conditions are set:

u ( 0 , t ) = 0.2 1 -cos ( π t ) 2 , u ( 1 , t ) = 0.3 sin( π t ) 2 , u ( x , 0 ) = 0 .
The time domain is defined as t ∈ 0 , 3 . According to Equation (3.5), the Scharfetter-Gummel scheme CFL condition is given in this case by:

∆t tanh a 0 + k 1 • ∆x 2 d 0 a 0 + k 1 • ∆x .
First, the problem is solved by implementing an adaptive in time approach, using Matlab™ function ODE113 with absolute and relative tolerances set to 10 -4 . The space discretization parameter is set ∆x = 10 -2 . Given the numerical value of the parameter, the CFL condition is ∆t 10 -3 . Figure 3 shows the variations of the solution u ( x , t ) as a function of space and time. It demonstrates a satisfactory agreement between the Scharfetter-Gummel numerical solution and the reference one. Figure 4 highlights the L 2 error value ε 2 ( x ) at the order of 10 -4 , which validates the numerical model for this case. This result is consistent with the value of the discretization parameter and the tolerances set in the ODE solver. An analysis of the numerical error is performed by varying one of the discretization parameters ∆x and ∆t values, while the other parameter remains constant. Different numerical schemes are compared for the computation of the solution of problem (3.1). The Scharfetter-Gummel scheme is implemented with (i ) the adaptive time step approach using ODE113 solver and a tolerance set to 10 -4 and (ii ) the Euler explicit approach. For comparison, a central finite-difference scheme using the Euler explicit method is also used.

Figure 5(a) shows the variation of the error ε 2 as a function of ∆x with a fixed time discretization parameter ∆t = 10 -4 . The CFL condition ∆t 2.5 • 10 -4 is respected until ∆x = 5•10 -3 . Beyond this limit, the numerical model cannot provide a bounded solution. These results also confirm a higher accuracy of the numerical model using the Scharfetter-Gummel scheme. Moreover, the Scharfetter-Gummel numerical scheme seems to have a second-order accuracy in space O( ∆x 2 ) which is better than the theoretical results mentioned in Section 3.2. The variations of the computational time with the requested level of accuracy is shown in Figure 5(b). For the accuracy ε ∞ = O ( 10 -3 ) , the finitedifference is ten-time slower than the Scharfetter-Gummel approach combined with the Euler explicit method. It can be noted that the Scharfetter-Gummel scheme with the adaptive time step is the fastest one thanks to the nonlinearity of the stability condition and adaptive time step. The computational time can be reduced by a hundred times. Figure 5(c) displays the variation of the significant correct digits with the computational time. In order to reach scd = O( 1.5 ) the finite-difference scheme requires ten times larger CPU time than the Scharfetter-Gummel scheme with Euler explicit approach. The latter is ten times slower than the one with an adaptive approach. Figure 6 presents another error study conducted by varying ∆t while the space discretization is fixed to ∆x = 10 -2 . For the given parameters, the CFL condition is ∆t 10 -3 . As noticed in Figure 6, a bounded solution can be computed until this condition is respected. It can also be noted that the Scharfetter-Gummel scheme with Euler explicit approach has a greater degree of accuracy. Moreover, the scheme is first-order accurate in time O ( t ) as mentioned by the theoretical results.

Nonlinear case

The previous validation is performed for a linear case in order to verify the theoretical results and advantages of the Scharfetter-Gummel numerical scheme. Since, the parameter estimation problem to be solved in Section 5 considers nonlinear coefficients, this Section aims at validating the numerical model for such cases. We consider the problem (2.4) with the following nonlinear coefficients: The dimensionless numbers are set to unity. The initial and boundary conditions are:

a ( u ) = 0.02 + 0.01 u 2 , d ( u ) = 0.05 + 0.03 u 2 , k ( u ) = 0.05 u 2 .
u ( x , 0 ) = 0 , u ( 0 , t ) = 0.8 sin π t / 3 + 0.2 sin π t / 5 , u ( 1 , t ) = 0.5 sin π t / 4 + 0.3 sin π t / 7 .
The solution of this problem is computed for the time horizon t ∈ 0 , 15 . The Scharfetter-Gummel numerical scheme is used with a space discretization step ∆ x = 0.01 and an adaptive time step with error tolerances set to 10 -4 . As shown in Figure 7, the profiles of the solution shows a very satisfactory agreement between the Scharfetter-Gummel solution and the reference one. It validates the proposed numerical model. The computational time of the numerical model is compared to three other ones: (i) the Scharfetter-Gummel scheme with Euler explicit approach, (ii) the central finite-differences approach with Euler explicit and (iii) the central finite-differences with adapting time step. As noticed in Table 1, the Scharfetter-Gummel with the adaptive time step scheme is the faster than Scharfetter-Gummel with Euler explicit approach. The computational time is of the same order between the two adaptive in time approaches (Scharfetter-Gummel and central finitedifferences). However, the Scharfetter-Gummel with adaptive time step approach is ten times more accurate than central finite-difference with adaptive time step.

As mentioned in Section 3.3, the numerical model requires the computation of the coefficients d

(a) (b) (c) (d) (f) (e) (g) ε ∞ 2 • 10 -3 1.9 • 10 -3 2 • 10 -3 1.9 • 10 -3 2 • 10 -3 1.9 • 10 -3 1.
k j+ 1 2 = k u j + u j+1 2 , (3.7a) k j+ 1 2 = 1 2 k u j + k u j+1 , (3.7b)
k j+ 1 2 = 2 k u j • k u j+1 k u j + k u j+1 , (3.7c)
k j+ 1 2 = 1 3 k u j + k u j • k u j+1 + k u j+1 , (3.7d)
k j+ 1 2 = k u j • k u j+1 , (3.7e)
k j+ 1 2 = 1 2 k u j 3 + k u j+1 3 1 3 , (3.7f) k j+ 1 2 = 1 2 k u j 2 + k u j+1 2 . (3.7g)
For each interpolation expression, the CPU time to compute the solution and the error with the reference solution are analyzed. The results are synthesized in Table 2 and Figure 8. All interpolation expressions have a similar order of accuracy. Nevertheless, it highlights that the approach using the arithmetic mean between u j and u j+1 to compute the coefficients d is slightly faster for this particular case.

This section examined and analyzed the Scharfetter-Gummel scheme applying to the diffusion-advection equation with linear and nonlinear coefficients. These results demonstrate a significant reduction of the computational cost without losing the accuracy while using the Scharfetter-Gummel scheme with an adaptive time step approach. Therefore, in further sections, the solutions are computed with the Scharfetter-Gummel combined with an adaptive time step approach. Since the numerical model is validated, the next section introduces the experimental facility employed to obtain the observations used in the parameter estimation problem.

Experimental facility

This section presents the facility to produce the experimental observations that will be compared to the numerical predictions of the model to compare its reliability to represent the physical phenomena.

Experimental observations

To perform the liquid water uptake tests, a refractory brick measuring 11 × 22 × 5 cm is used. The brick is initially sanded allowing a better visualization of the rise of the water. The vertical faces are then taped to be airtight and watertight. As illustrated in Figure 9, a tray is filled with distilled water and maintained at a constant level thanks to a bottle returned with a pierced cap. The brick is placed vertically on a support at the surface of the water. The experiment begins at the moment when the underside of the brick is in contact with the water and lasts at least 7 h . The rising damp is observed on one of the faces 11 × 22 cm. For this, a camera is placed at 1 m to take a picture of the brick every 3 min . The experiment occurs in a box opaque to the light to control the luminosity as well as to minimize the convective and radiative exchanges. Some LED lamps provide a constant lighting in the box.

At the top of the brick, two different boundary conditions are imposed. For the first set of experiments, the top of the brick is in contact with the open air at atmospheric pressure. For the second set, a relative pressure of -50 Pa is applied thanks to a chamber with a fan. An illustration of the design is shown in Figure 10(a). With this device, the air velocity in the brick is assumed as constant in time. The experiment is taken at the constant room temperature T = 293 K.

The height of the water front is taken in the middle of the brick from the pictures of each experiment every 30 min . As shown in Figure 10(b), graph paper is attached to the brick in order to measure the height front. With this setting, the uncertainty measurement of the height scales with σ h = 0.5 cm .

For the case when the brick is exposed to ambient air pressure, the experiment takes t max = 15 h . For the other experiment, with a pressure difference, it lasts t max = 7 h . Figure 11 shows the obtained experimental data for both cases.

To verify that the liquid front is homogeneous inside the brick and not occurring only at the interfaces where the pictures are taken, a complementary experiment is performed. The brick is preliminary cut in the middle and then fulfilled with silicone to avoid the liquid uptake in this area. The liquid uptake experiment is performed. At the end of the experiments, the cut enables to rapidly break the brick without perturbing the liquid front. As illustrated in Figure 12, it is confirmed that the liquid front is homogeneous inside the brick. The white silicone that fulfills the cut in the middle of the brick can be noticed.

Defining the material properties

The physical model representing the liquid water uptake test is described in Section 2. Now, it is important to define the material properties involved in the definition of the diffusion, gravity and advection fluxes for the brick used in the experiments. These properties are uncertain [14] and will be determined in the next Section using the experimental data obtained with the facility.

Data from the literature are used to define the dependency of these properties on the moisture content. The data is projected on the polynomial functions, using the cftool from the Matlab™ environment, and polynomials with the highest root mean square error Rsquare (the square of the correlation between the response values and the predicted response values) are selected. These measurements provide a priori values of the coefficients of the polynomials, which act as initial guess when solving the parameter estimation problem. The polynomial functions have been chosen for the material properties to fit the data from the literature. Another set of functions can be applied, and the same steps of parameter estimation are utilized. It should be noticed that data fitting is performed only for the liquid water content.

The reference parameters used for computations are With given reference values, dimensionless representation d ( u ) of the liquid transport coefficient is:

L = 0.22 m , t ref = 1 h , D ref = 10 -6 m 2 /s , k ref = 10
d ( u ) = d 4 u 4 + d 1 u + d 0 ,
where d 4 = 0.6 , d 1 = -0.04 and d 0 = 0.0067 . The dimensionless function a ( u ) of the advection part:

a ( u ) = a 0 1 -H ,
where a 0 = 0.7 , and H = ´1 0 u ( x , t ) dx . Last, the dimensionless function k ( u ) of the liquid conductivity part:

k ( u ) = k 3 u 3 ,
where k 3 = 0.8 . Let us discuss which coefficient of d 0 , d 1 or d 4 contributes the most in the context of liquid uptake or specifically u = O ( 1 ) . For this purpose, the normalized derivative

d i • ∂d ( u ) ∂d i
relative to each coefficient is taken and compared with each other. The results are presented in Figure 13. One may conclude that coefficient d 4 has more influence on the value of function d ( u ), when liquid uptake occurs u = O( 1 ) . Therefore, in this article only coefficient d 4 is estimated.

d i • ∂d ( u ) ∂d i = { d 4 u 4 ; d 1 u ; d 0 } . u(-) 0

Comparison of the numerical predictions with experimental data

The efficiency of the proposed numerical model is demonstrated in Section 3.5 by comparison to reference solution. Now, the purpose is to assess the reliability of the numerical model to predict the physical phenomena. For this, the experimental observations described in previous section are used for comparison with the numerical predictions. With the given parametrization in Section 4.2, Eq.(2.4) writes as:

∂u ∂t = Fo • ∂ ∂x d ( u ) • ∂u ∂x -Pe • a ( u ) • u -Bo • k ( u ) , (5.1)
where

d ( u ) = d 4 u 4 + d 1 u + d 0 , a ( u ) = a 0 1 - ˆ1 0 u ( x , t ) dx , k ( u ) = k 3 u 3 .
Since literature lacks of consistent data [14], the material properties, represented by functions k ( u ) , d ( u ) and a ( u ) , are uncertain. Thus, parameter estimation problem are solved to determine the uncertain parameters and calibrate the numerical model with the experimental data. Within the procedure of parameter estimation, first we demonstrate the structural and practical identifiability of the three unknown parameters a 0 , k 3 , d 4 .

Then, the parameter estimation problem is solved and the reliability of the numerical model is discussed.

Structural Identifiability

The aim of this section is to demonstrate the formal identifiability of the unknown parameters. The set of unknown parameters is defined as:

P ≡ { P i } = { a 0 , k 3 , d 4 } .
We assume here, u ( x , t ) is the only observable field.

A parameter P i ∈ P is Structurally Globally Identifiable (SGI), if the following condition is satisfied [START_REF] Walter | Global approaches to identifiability testing for linear and nonlinear state space models[END_REF]:

∀t , u ( P ) = u ( P ⋆ ) ⇒ P i = P ⋆ i .
In

Fo • ∂ ∂x d( u ) • ∂u ∂x -Pe • a( u ) • u -Bo • k( u ) ≡ (5.2) Fo • ∂ ∂x d ⋆ ( u ) • ∂u ∂x -Pe • a ⋆ ( u ) • u ⋆ -Bo • k ⋆ ( u ) .
In addition, straightforward replacements in the Eq. (5.1) shows us the following:

∂u ∂t = Fo • 4 d 4 u 3 + d 1 • ∂u ∂x + d 4 u 4 + d 1 u + d 0 • ∂ 2 u ∂x 2 -Pe • a 0 • 1 - ˆ1 0 u dx • ∂u ∂x -3 Bo • k 3 • u 2 • ∂u ∂x ,
which shows that Equality (5.2) becomes:

Fo • 4 d 4 u 3 + d 1 • ∂u ∂x + d 4 u 4 + d 1 u + d 0 • ∂ 2 u ∂x 2 -Pe • a 0 • 1 - ˆ1 0 u dx • ∂u ∂x -3 Bo • k 3 • u 2 • ∂u ∂x ≡ Fo • 4 d ,⋆ 4 u 3 + d 1 • ∂u ∂x + d ⋆ 4 u 4 + d 1 u + d 0 • ∂ 2 u ∂x 2 -Pe • a ⋆ 0 1 - ˆ1 0 u dx • ∂u ∂x -3 Bo • k ⋆ 3 • u 2 • ∂u ∂x ,
and the final expression is

d 4 -d ,⋆ 4 • 4 u 3 + u 4 ∂ 2 u ∂x 2 -Pe • a 0 -a ,⋆ 0 • 1 - ˆ1 0 u dx • ∂u ∂x -3 Bo • k 3 -k ,⋆ 3 • u 2 • ∂u ∂x ≡ 0 Since that 4 u 3 + u 4 ∂ 2 u ∂x 2 , 1 - ˆ1 0 u dx • ∂u ∂x , u 2 • ∂u ∂x

Practical Identifiability

As demonstrated in previous Section, the three unknown parameters d 4 , k 3 and a 0 are globally identifiable. This result remains theoretical. Before solving the parameter estimation problem, it is of major importance to study the sensitivity coefficients of the parameters to demonstrate the practical identifiability. If the sensitivity coefficients are either small or correlated, the estimation problem is difficult and very sensitive to measurement errors. The sensitivity coefficient is defined as the first derivative of the (numerical) observations with respect to an unknown parameter [9,35]:

Y P i ( t ) = σ p σ H ∂H ∂P i = σ p σ H ˆ1 0 X P i dx , X P i = ∂u ∂P i ,
where σ H and σ p are scaling factors. In this case, σ H corresponds to the uncertainty on the experimental observations. The quantity σ p is set to unity since all unknown parameter have the same uncertainty on the a priori values.

To compute the sensitivity coefficients, Equation (5.1) is differentiated with respect to each unknown parameter. Three differential equations are obtained enabling to compute the sensitivity coefficients. For the parameter k 3 , we define X k 3 def := ∂u ∂k 3 obtained by solving the following equation:

∂X k 3 ∂t = Fo • ∂ ∂x d ( u ) • ∂X k 3 ∂x + d ( u ) • ∂u ∂x -Pe • a ( u ) -Bo • k ( u ) • X k -Bo u 3 + Pe a 0 ˆ1 0 X k 3 dx ,
where

d ( u ) = 4 d 4 u 3 + d 1 , k ( u ) = 3 k 3 u 2 .
For the parameter d 4 , we define X d 4 def := ∂u ∂d 4 obtained by solving the following equation:

∂X d 4 ∂t = Fo • ∂ ∂x d ( u ) • ∂X d 4 ∂x + d ( u ) • ∂u ∂x -Pe • a ( u ) -Bo • k ( u ) • X d + u 4 ∂u ∂x + Pe a 0 ˆ1 0 X d 4 dx .
And last, for the parameter a 0 , we define X a 0 def := ∂u ∂a 0 obtained by solving the following equation:

∂X a 0 ∂t = Fo • ∂ ∂x d ( u ) • ∂X a 0 ∂x + d ( u ) • ∂u ∂x -Pe • a ( u ) -Bo • k ( u ) • X a -Pe • 1 - ˆ1 0 u( x , t ) dx -a 0 ˆ1 0 X a 0 dx • u .
The boundary and initial conditions are the same for each coefficient X P i :

X P i ( x = 0 , t ) = 0 , X P i ( x = 1 , t ) = 0 , X P i ( x , t = 0 ) = 0 .
The three differential equations are solved with Scharfetter-Gummel numerical scheme.

For the investigation of the practical identifiability, the a priori parameters values are used. It is not correlated with other coefficients, demonstrating that it is possible to estimate it. However, its magnitudes are very small showing that its estimation with satisfactory accuracy is a difficult task. One can argue that these observations are due to the a priori values for the coefficient a 0 . However, as noticed in Figure 14(b), even multiplying by 15 or 100 times the a priori value of a 0 , the magnitude of the sensitivity coefficient almost doesn't change.

Estimation and comparison with the experimental observations

Since the formal and practical identifiability of the three unknown parameters d 4 , a 0 and k 3 have been demonstrated, the aim is now to estimate them using the experimental observations obtained with the facility. First, the methodology to solve the inverse problem is briefly described. Then, the results are presented and the reliability of the numerical model predictions are discussed.

Methodology to solve the parameter estimation problem

The parameter estimation problem is solved by minimizing the following cost function by the least squares method:

J d 4 , k 3 , a 0 : = H d 4 , k 3 , a 0 -H exp 2 .
The value of H results from the solution of the direct problem (5.1) for a given set of parameters d 4 , k 3 , a 0 . The value of H exp is given by the measurements from the experimental facility and interpolated on the time grid of the numerical scheme.

The cost function J is minimized using function fmincon from the Matlab™ environment. This method uses the interior-point algorithm with constraints on the unknown parameters. Here, the upper and lower boundary constraints for the parameters are defined based on the preliminary calculations:

d 4 ∈ 0, 2 , k 3 ∈ 0, 2 , a 0 ∈ 0, 2 .
To estimate the quality of the solution of parameter estimation problem, the normalized Fisher matrix is defined according to [19,33]:

F = F ij , ∀ i, j ∈ { 1, ..., N p } , F ij = 1 σ 2 H ˆtmax 0 Y P i Y P j dt ,
where Y P i is the sensitivity coefficient of the solution related to the parameter P i , σ H the measurement uncertainty and N p the number of parameters. The matrix F measures the total sensitivity of the system for the measurements to variations of the entire set of parameters P. Under some assumptions detailed in [35], the inverse of the Fisher matrix is the matrix of variance of the parameters considered as random variables of the given observable fields. In other words, it summarizes the quality of the information obtained in the parameter identification process. Thus, the inverse matrix of F is used to assess the estimation uncertainty by computing an error estimator for the parameter P i :

η i = F -1 i i .
High values of η i indicate a possible high error during the parameter identification process.

Results and discussion

The estimation process is performed with the following discretization parameters ∆ x ⋆ = 0.05 , ∆ t ⋆ = 10 -2 , using the Scharfetter-Gummel numerical scheme and the Adams-Bashforth-Moulton time adaptive algorithm. The solution of the optimization problem requires the computation of the solution u of Eq. (2.4) for each set of the parameters ( d 4 , k 3 , a 0 ). During the optimization process for the case without pressure, the direct problem Eq. (2.4) is solved 114 times. Each iteration requires approximately 100 s to solve the direct problem. It results in a overall performance of 3 h . Considering the results from previous section and particularly Table 1, the same process using the standard finite-difference scheme and Euler approach would last for about 7.5 h . For the experimental data with pressure, the calculation of solution u of Eq. (2.4) is computed almost 270 times. The estimation takes about 7.5 h . This value is comparatively small to the 17.5 h required with the model based on standard approach. The results highlights that the Scharfetter-Gummel numerical scheme, applying to parameter estimation problem, save important computational efforts compared to the standard approach with finite-difference scheme.

Figure 15 shows the time variation of the water height for both experiments with pressure and without pressure. It compares the solution computed using the estimated values of the parameters and the experimental data. The results demonstrate a good agreement with the experimental observations, while the solution with the a priori values cannot approximate the experimental data well. One may say that the calibrated numerical model has a As a result of the optimization process, the estimated values, the residuals of the cost function and the error estimators are reported in the Table 3. The error estimator is small for parameters k 3 and d 4 , proving that the accuracy of estimation is decent. However, the high values of the error estimator of parameter a 0 indicate an unsatisfactory estimation with a very high uncertainty. These results are consistent with ones obtained when analyzing the sensitivity coefficient with a priori values in Section 5.2. In addition, as noticed in Table 3, there is no significant difference between the a priori and estimated In order to demonstrate the convergence of the optimization algorithm, the following procedure is implemented. The parameter a 0 cannot be estimated accurately so its value is fixed for this analysis. Then several values of initial guesses for parameters d 4 and k 3 are taken. The estimation procedure is performed for each case. Figure 16 shows the path from each starting point until the estimated parameters values. As one may conclude the final points are in a similar range for all the tested initial guesses. It validates the convergence of the optimization method. It should be noted that the figure does not show the intermediary computations of the algorithm procedure.

Figure 17 displays the volume of liquid water relative to space and time when the brick is exposed to the open air or to a pressure difference of -50 Pa . The propagation of the water front through the brick seems faster in the case of the pressure difference. As noted in Figure 15, the experimental data with pressure show that the velocity of water uptake is faster since the same height of the water is achieved twice faster than with the experiment without pressure. One may conclude that the pressure influences the velocity of water uptake. However, as reported in indicate that material properties may be different between the two bricks used during the experiments. Moreover, the error estimators of the advection coefficients are really poor indicating that for the given value of pressure difference, it has a minor influence on the water uptake. It should also be noted that the sensitivity coefficients of parameter a 0 has really low magnitudes of variations as shown in Figure 14(b).

In Figure 15, it can be noticed that the full length of the brick is not reached, validating the chosen physical model with Dirichlet boundary conditions presented in Section 2.

One may argue that other types of boundary conditions could have been chosen at the top of the brick x = 1 . The Robin-type boundary condition is avoided. Indeed the boundary flow is proportional to the ambient and surface conditions:

d ( u ) • ∂u ∂x -Pe • a 0 • 1 - ˆ1 0 u ( x , t ) dx • u -Bo • k ( u ) = Bi • u -u ∞ , x = 1 ,
where Bi is a dimensionless surface transfer coefficient and u ∞ are the ambient conditions in the facility. It implies to estimate this coefficient Bi which is rather difficult experimentally. It can be estimated by adding the coefficient Bi in the inverse problem, at the price of additional computational costs. In addition, the expected values of the surface coefficient are slow and the boundary conditions can be modified into a homogeneous Neumann one:

d ( u ) • ∂u ∂x -Pe • a 0 • 1 - ˆ1 0 u ( x , t ) dx • u -Bo • k ( u ) = 0 , x = 1 , (5.3)
To investigate the importance of modifying the definition of the boundary conditions at x = 1 , additional computation are performed considering the problem Eq. (5.1) with the estimated parameters ( d 4 , k 3 , a 0 ) and the boundary condition Eq. (5.3). The boundary condition at x = 0 and the initial condition are unchanged: Another open question is the definition of the liquid height. Previously, the height is defined by the integral of u over the space domain. An alternative definition is to the height of the water front in the brick using a threshold:

u( x = 0 , t) =    0 , t = 0 1 , t > 0 , u( x , t = 0) = 0 .
H = max x : u ( x , t ) ≥ û , (5.4)
where û is a chosen value of the amount of liquid. Figure 19 compares the numerical solution according to the definition of H with Eq. (2.3) or (5.4). The numerical predictions are varying with the definition. On one hand, the integral definition is a more general approach but it assumes that the water vapor is negligible in the material. On the other hand, the definition (5.4) requires to determine the value of the threshold. For these reasons, the definition of the liquid height in the brick according remains an open question.

Conclusion

The moisture in the walls greatly affects the overall building performance. Such a process as rising damp has a strong impact on the energy consumption of the building. In this work, the physical model represents the water uptake process in a single brick, based on diffusion, advection and gravity fluxes. To investigate these physical phenomena it is important to have efficient numerical models in terms of computational cost and accuracy. These models require also to have a good reliability to represent the physical phenomena. To answer this issue, an innovative approach is proposed based on the Scharfetter-Gummel numerical scheme. The first part of the paper studies the properties of the numerical model such as accuracy, stability conditions and computational time with a reference solution for both linear and nonlinear cases in Section 3.5. The investigations show that the Scharfetter-Gummel numerical model is more accurate and faster than the standard approach based on finite-differences and Euler explicit scheme. The Scharfetter-Gummel numerical model has an explicit formulation avoiding costly sub-iterations at each time step. Since its stability condition is nonlinear and scales with ∆t ≃ ∆x for large space discretization, the Scharfetter-Gummel numerical scheme combined with an adaptive time step approach is particularly efficient to save computational efforts. Along with numerical methods, the reliability of the proposed numerical model is analyzed by comparing the numerical predictions to experimental observations. The experimental facility is presented in Section 4. Empirical data are gathered through a set of two experiments. First, the height of the rising front is obtained for a brick under normal conditions. Then, the height front is measured in another brick submitted to a pressure of -50 Pa . Then, the comparison between numerical simulations and the experimental observations of the height of the rising front is carried out in Section 5. Since the literature lacks consistent data for the material properties [14], the uncertain parameters are determined by solving a parameter estimation problems. Before solving the parameter estimation problem, the structural and practical identifiability of the three unknown parameters are demonstrated in Sections 5.1 and 5.2, respectively. The results of the parameter estimation problem show that the diffusivity and liquid conductivity are estimated with satisfactory accuracy. The advection coefficient cannot be estimated with accuracy even with a pressure difference of -50 Pa at the top of the brick. The sensitivity of the numerical solution with respect to this parameter is much smaller compared to the diffusivity and liquid conductivity. With the estimated parameters, there is a satisfactory agreement between the numerical predictions and the experimental observations, highlighting a good reliability of the model. Important computational efforts are saved thanks to the efficiency of the numerical model. The Scharfetter-Gummel numerical scheme enables to save by 50% the computational cost compared to the standard approach.

Even if this work enhances the efficiency of the Scharfetter-Gummel numerical scheme for the solution of advection-diffusion equation with gravity flux, the rising damp problem require simulation of the physical phenomena in 2-dimensions. Thus, further research should be conducted to extend the numerical model in this way.

In order to study the uniqueness of the solution, we use the so-called energy method [START_REF] Evans | Partial Differential Equations[END_REF]. For that purpose, we assume that two solutions u 1 ( x , t ) and u 2 ( x , t ) satisfy Eqs. (A.1) and (A.2). We define w( x , t ) as: w ( x , t ) The Cauchy-Schwartz inequality writes:

ˆΩ f g ˆΩ f 2 dx 1 2 ˆΩ g 2 dx 1 2
, for f and g square-integrable realx-value functions. Thus, the inequality (A.5) yields to:

Figure 1 .

 1 Figure 1. Illustration of the problem of water uptake in the brick.

Figure 4 .

 4 Figure 4. Error ε 2 as a function of space.

Figure 5 .

 5 Figure 5. (a) Variation of the error ε ∞ as a function of ∆x . (b) Variation of the error ε ∞ as a function of the CPU time of the numerical model. (c) Variation of the accuracy digits value scd as a function of the CPU time of the numerical model.

Figure 6 .

 6 Figure 6. Variation of the error ε ∞ as a function of ∆t .

  of the dual cell x j , x j+1 . According to Eq. (3.6), these coefficients are interpolated using the mean values of u j and u j+1 . The accuracy of this interpolation is tested by using different expressions in the computation of the solution u . The coefficients d are calculated by one the following Table2. Variation of the error of the solution and of the CPU time of the numerical models according to the interpolation formulas Eq. (3.7) of the coefficients d ( u ) , a ( u ) and k ( u ) .

Formula

Figure 7 .Figure 8 .

 78 Figure 7. Variation of the field u as a function of x (a) and t (b) with nonlinear coefficients

Figure 9 .

 9 Figure 9. Illustrations of the experimental facility. It should be noted that the opaque box has been removed to take the pictures and that the light brick is not the one considered for the present studies.

Figure 10 .Figure 11 .Figure 12 .

 101112 Figure 10. (a) Illustrations of the brick submitted to a relative pressure of -50 Pa linking the top of the brick to a fan. (b) Measurement of the water front height using the graph paper.

  -10 m/s and a ref = 10 -9 m/s . The saturated volumetric moisture content θ sat , introduced in [20], equals θ sat = 0.3065 m 3 /m 3 . It induces the following dimensionless numbers: Fo = 0.074 , Pe = 2.2 • 10 -4 , Bo = 0.7 .

Figure 13 .

 13 Figure 13. The normalized derivative of d ( u ) relative to each coefficient d 0 , d 1 and d 4

  are independent, then d 4 , k 3 and a 0 are SGI [23].

The time values of Y k 3

 3 Figure 14(a). One may conclude that the estimation of the parameters k 3 and d 4 will likely have a good result since their sensitivity coefficients have large magnitudes and they are not correlated with each other. On the other hand, the time variation of the third coefficient Y a 0 is given in Figure 14(b).

Figure 14 .

 14 Figure 14. Time variation of the sensitivity coefficients (a) Y k 3 , Y d 4 and (b) Y a 0 .

Figure 15 .

 15 Figure 15. Water uptake level in the brick respective to the time for case without pressure (a) and (b) case with pressure. Figure displays experimental data with uncertainty σ h = 0.5 cm and its interpolation, a priori data and results with the estimated parameters with residual ε 2 = 0.05 for (a) case and ε 2 = 0.0087 for (b) case.

Figure 16 .

 16 Figure 16. Illustration of optimization procedure convergence.

Figure 18

 18 Figure18displays the difference between the solutions computed with Dirichlet and Neumann boundary conditions. The error ε 2 is lower than 10 -4 and one may conclude that the boundary conditions do not change significantly the solution of the numerical model. These results are consistent when analyzing Figures17(c) and 17(a). It can be observed that the flow at x = 1 remains null.Another open question is the definition of the liquid height. Previously, the height is defined by the integral of u over the space domain. An alternative definition is to the height of the water front in the brick using a threshold:

Figure 18 .

 18 Figure 18. Error ε 2 between solutions with Dirichlet and homogeneous Neumann boundary condition at x = 1 .

Figure 19 .

 19 Figure 19. Influence of the definition of the water front height of the brick on the numerical predictions.

def:= u 1 (

 1 x , t ) -u 2 ( x , t ) , solution of the equation: and initial conditions:w ( 1 , t ) = 0 , w ( 0 , t ) = 0 , w ( x , 0 ) = 0 . (A.4) If the energy E ( t ) def := ˆΩ w 2 ( x , t) dx is decreasing and E ( 0 ) = 0 , then the solution of Eqs (A.1) and (A.2) is unique :dE dt 0 .To study the validity of this condition, both sides of Eq. (A.3) are multiplied by w and integrated over Ω : dx , By performing an integration by parts, we obtain:

Table 1 .

 1 Computational time required to compute the solution of the nonlinear case with the different numerical models.

	Numerical model	CPU time [ s ] CPU time [ % ]	ε ∞
	SG with adaptive time step	105	42	1.9 • 10 -4
	SG with Euler explicit	249	100	5.8 • 10 -5
	Central FD with adaptive time step	101	40	1.1 • 10 -3
	Central FD with Euler explicit	230	92	2.1 • 10 -2

Table 3 .

 3 A priori and estimated values of dimensionless coefficients

	Parameters	A priori values Estimated value Error η i Estimated value Error η i No pressure With pressure
	a ⋆ 0	0.7	0.0052	± 14	0.55	± 206
	k ⋆ 3	0.8	0.8257	± 0.01	1.3	± 0.21
	d ⋆ 4	0.6	1.0	± 0.002	1.72	± 0.02
	Residual ε 2	1	0.05		0.0087
	satisfactory reliability to represent the physical phenomena of liquid water uptake. Similar
	conclusions can be noted for the experiments with pressure.

Table 3 ,

 3 results from the optimization problem gives different values for the set of parameters ( d 4 , k 3 , a 0 ) with and without pressure. Particularly, values for parameters d 4 and k 3 vary from one case to another. It may Figure 17. Variation of the dimensionless amount of liquid water u in case with and without pressure. Figures (a) and (b) displays case without pressure, relative to the height with t ∈ { 2 , 4 , 6 , 8 , 10 , 12 } , and time with x ∈ { 0.2 , 0.4 , 0.6 , 0.8 } respectively. Figures (c) and (d) displays case with pressure, relative to the height with time t ∈ { 2 , 4 , 6 , 7 } , and time with x ∈ { 0.2 , 0.4 , 0.6 , 0.8 } respectively.

		1						1
		0.8						0.8
	( -)	0.6	t > 0				( -)	0.6				x > 0
	u	0.4					u	0.4
		0.2	t =					0.2	x = 0.2
			2
		0						0
		0	0.2	0.4	0.6	0.8	1	0	2	4	6	8	10	12	14	16
				x ( -)						t ( -)
				(a)							(b)
		1						1
		0.8						0.8
	( -)	0.6	t > 0			( -)	0.6				x > 0
	u	0.4					u	0.4
		0.2	t =					0.2	x = 0.2
			2
		0						0
		0	0.2	0.4	0.6	0.8	1	0		2		4		6
				x ( -)						t ( -)
				(c)							(d)

Acknowledgements

This work was partly funded by the "Conseil Savoie Mont Blanc" (CSMB), the French Atomic and Alternative Energy Center (CEA) and the French Environment and Energy Management Agency (ADEME) (through the research program CAPVENT). The authors also acknowledge the Junior Chair Research program "Building performance assessment, evaluation and enhancement" from the University of Savoie Mont Blanc in collaboration with the French Atomic and Alternative Energy Center (CEA) and Scientific and Technical Center for Buildings (CSTB). The authors also would like to acknowledge Dr. B. Rysbaiuly for his valuable discussions.

A. Demonstration of uniqueness of the solution

This section proposes to demonstrate the uniqueness of the solution of liquid uptake in porous media. Since u ( x , t ) > 0 at t > 0 , equation (3.1) can be transformed as

then regularized into a general advection-diffusion equation:

where ã ( u

. Further in the Section, we assume the linearized representation ã of ã ( u ) around u 0 :

where u 0 is a given function of ( x , t ) .

The following boundary and initial conditions are considered:

For u ∈ W 1 , p 0 , the Poincaré inequality states:

, where C Ω is a constant depending on p and Ω only. Thus, for a one-dimensional problem and for p = 2 , the inequality (A.6) becomes:

which can be re-written as:

Therefore, we have uniqueness of the solution of Eqs (A.1) and (A.2) if k -ã C Ω 0 . In other terms, the generalized advection coefficient has to be small compared to the diffusion one to have unique solution.