An efficient method to estimate sorption isotherm curve coefficients

ABSTRACT

This paper aims at estimating the sorption isotherm coefficients of a wood fiber material using experimental data. First, the mathematical model, based on convective transport of moisture, the Optimal Experiment Design (OED) and the experimental set-up are presented. Then, measurements of relative humidity within the material are carried out, after searching the OED using the computation of the sensitivity functions and a priori values of the unknown parameters. It enables to plan the experimental conditions in terms of sensor positioning and boundary conditions out of 20 possible designs, ensuring the best accuracy for the identification method and, thus, for the estimated parameter. After the measurements, the parameter estimation problem is solved. The determined sorption isotherm coefficients calibrate the numerical model to fit better the experimental data. However, some discrepancies still appear since the hysteresis effects on the sorption capacity are not included in the model. Therefore, the latter is improved proposing an additional differential equation for the sorption capacity to consider the hysteresis effects. The OED approach is developed for the estimation of five of the coefficients involved in the hysteresis model. To conclude, the prediction of the model with hysteresis have better reliability when compared to the experimental observations.

KEYWORDS:

• Optimal Experiment Design (OED)
• parameter estimation problem
• convective moisture transfer
• sensitivity functions
• sorption moisture coefficients
• hysteresis

MSC SUBJECT CLASSIFICATIONS:

• 65N21
• 80A23
• 35K55

1. Introduction

Moisture in buildings has been a subject of major concern since 1980s. It may affect energy consumption and demand so we can mention at least four International Energy Agency projects conducted in the last 30 years to promote global research on this subject (Annexes 14, 24, 41 and 55) [1]. Furthermore, moisture can also have a dramatic impact on occupants' health and on material deterioration. Several tools have been developed to simulate the moisture transport in constructions as described in [2], which can be used to predict conduction loads associated to porous elements and mould growth risk in building enclosures. Nevertheless, those tools require input parameters containing temperature- and moisture-dependent hygrothermal properties.

1.1. Moisture transport in constructions

The following system of partial differential equations established by Luikov [3] represents the physical phenomenon of heat and mass transfer through capillary porous materials: (1) $\frac{\mathrm{\partial }U}{\mathrm{\partial }t}=\mathbf{\nabla }\mathbf{\cdot }\left(a_m\mathbf{\nabla }U+\delta a_m\mathbf{\nabla }T\right),$(1) (2) $c_b{\rho }_0\frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\mathbf{\nabla }\mathbf{\cdot }\left(\lambda \mathbf{\nabla }T\right)+r_{12}\mathbf{\nabla }\mathbf{\cdot }\left(a_{m1}{\rho }_0\left(\mathbf{\nabla }U+{\delta }_1\mathbf{\nabla }T\right)\right),$(2) where U is the the specific mass content in the porous body, corresponding to the ratio of mass of water (in vapour and liquid phases) to the dry-basis mass of porous material. The other quantities are T the temperature, $a_m$ the mass transfer coefficient for vapor (denoted with the subscript 1) and liquid inside the body, δ the thermal-gradient coefficient, ${\rho }_0$ the specific mass of the dry body, $c_b$ the specific heat of the body and $r_{12}$ the latent heat of vapourization.

In building physics, those equations represent the physics that occurs in the porous building envelope and indoor porous elements such as furniture, textiles, etc. Regarding the envelope, the phenomenon is investigated to analyse the influence of moisture transfer on the total heat flux passing through the wall, with the objective of estimating heat losses. They are also studied to analyse the durability of walls and to avoid disorders due to the presence of moisture as, for instance, mould growth, shrinking or interstitial condensation. This aspect is of major importance for wall configurations involving several materials with different properties, where moisture can be accumulated at the interface between two materials. Durability problems may also appear when considering important moisture sources as wind-driven rain or rising damp problems. These analyses are performed by computing solutions to the partial differential equations similar to (1(1) $\frac{\mathrm{\partial }U}{\mathrm{\partial }t}=\mathbf{\nabla }\mathbf{\cdot }\left(a_m\mathbf{\nabla }U+\delta a_m\mathbf{\nabla }T\right),$(1) ). For this, numerical methods are used due to the nonlinearity of the material properties, depending on moisture content and temperature, and the non-stationary boundary conditions, defined as Robin-type and varying according to climatic data. Most of the numerical approaches consider standard discretization techniques. For the time discretisation, the Euler implicit [4–6] or explicit [7] schemes are adopted. Regarding the spatial discretisation, works in [8] are based on finite-differences methods, in [5,6,9,10] on finite-volume methods and in [11,12] on finite-element methods. It is important to note that the numerical solution of these partial differential equations requires, after discretisation, the calculation of large systems of nonlinear algebraic equations (an order of ${10}^6$ for 3-D problems). Furthermore, the problem deals with different time scales. The diffusive phenomena and the boundary conditions evolve on the time scale of minutes while the building performance usual analysis is done for a time interval of one year or even longer when dealing with durability or mould growth issues. Thus, the computation of heat and moisture transfer in porous material in building physics has a non-negligible computational cost. Recently, innovative and efficient methods of numerical simulation have been proposed. Some improved explicit schemes, enabling to overcome the stability restrictions of standard Euler explicit schemes, have been proposed in [13,14]. An accurate and fast numerical scheme based on the Scharfetter–Gummel idea has been proposed in [15] to solve the advective–diffusive moisture differential equation. Some attempts based on model reduction methods have been also proposed with an overview in [8].

1.2. Inverse problems in building physics

While some research focuses on numerical methods to compute the solution of the so-called direct problem to analyse the physical phenomena, some studies aim at solving inverse problems of heat and mass transfer in porous materials. In this case, the focus is the estimation of material properties $({\rho }_0^{\circ },c^{\circ },a_m^{\circ },{\delta }^{\circ })$ using experimental data denoted as $(T_{\exp },U_{\exp })$ by minimizing a thoroughly chosen cost function $\text{J}$: $\begin{array}{rl}\left({\rho }_0^{\circ },c^{\circ },a_m^{\circ },{\delta }^{\circ }\right)& =\text{arg}min\text{J},\\ \mathrm{with}\text{J}& =\left|\left|\left({\mathit{T}}_{\exp }-\mathit{T},{\mathit{U}}_{\exp }-\mathit{U}\right)\right|\right|,\end{array}$ where $||(\bullet ,\bullet )||$ is certain vector-norm in time.

Here, the inverse problem is an inverse medium problem, as it aims at estimating the coefficients of the main equation [16,17]. Two contexts can be distinguished. First, when dealing with existing buildings to be retrofitted, samples cannot be extracted from the walls to determine their material properties. Therefore, some in situ measurements are carried out according to a non-destructive design. The experimental data can be gathered by temperature, relative humidity, flux sensors and infrared thermography, among others. In most cases, measurements are made at the boundary of the domain. From the obtained data, parameter estimation enables to determine the material thermophysical properties. As mentioned before, the properties are moisture and temperature-dependent. Therefore, the parameter identification problem needs to estimate functions that can be parameterized. Moreover, in such investigations, there is generally a few a priori information on the material properties. In [18], the thermal conductivity of an old historic building wall composed of different materials is presented. In [19], the thermophysical properties of materials composing a wall are estimated. In [20], the heat capacity and the thermal conductivity of a heterogeneous wall are determined. Once the needed parameters have been estimated, efficient simulations using the direct model can be performed to predict the wall conduction loads and, at the end, choose adequate retrofitting options.

Another issue arises while comparing the numerical model results and experimental data. Some discrepancies were observed as reported in several studies [15,21] and illustrated in Figure 1. A material, with an initial moisture content $U_0$, is submitted to adsorption and desorption cycles. The moisture content raises up to $U_{\max }$ at the end of the adsorption phase. Then, during the desorption phase, the moisture content decreases until a value $U_1>U_0$ due to the hysteresis effects. When comparing the simulation results to the experimental data, it is observed that predictions of moisture content commonly underestimate the experimental observations of moisture adsorption processes. In other words, the numerical predictions are lower than the experimental values obtained during the adsorption phase. On the other hand, in the course of the desorption phase, the numerical predictions often overestimate the experimental observations, i.e. the simulation values are greater than the experimental ones. At the end of the desorption process, the experimental values are greater again, when compared to both the predicted values and the initial moisture content since the hysteresis phenomenon significantly affects the material moisture sorption capacity. It means the experimental moisture front rushes faster than the numerical predictions. To answer this issue, models can be calibrated using in situ measurements for adapting the material properties to reduce the discrepancies between model predictions and real observations. In [22], moisture- and temperature-dependent diffusivity and thermophysical properties are estimated using only temperature measurements under a drying process. In [23], the moisture permeability and advective coefficients are estimated using relative humidity measurements in a wood fibre material. In these cases, a priori information on the material properties is known thanks to complementary measurements based on well-established standards.

Figure 1. Illustration of the discrepancies observed when comparing experimental data to results from numerical models of moisture transfer in porous material.

In terms of methodology for the estimation of parameters, several approaches can be distinguished in literature. Descent algorithms, based on the Levenberg–Marquardt nonlinear least-square method, are used in [20]. Stochastic approaches, using Bayesian inferences and the Markov chain Monte Carlo algorithms, are applied in [18,24–26]. Approaches based on genetic algorithms based approaches are adopted to minimize the cost function in [27,28]. Model reduction techniques, based on Proper Orthogonal Decomposition , are employed in [29].

1.3. Problem statement

This article presents an efficient ethod for the estimation of moisture sorption isotherm coefficients of a wood fibre material, represented by three parameters, using an experimental facility, aiming at reducing the discrepancies between model predictions and real performance. The estimation of the unknown parameters, based on observed data and identification methods, strongly depends on the experimental protocol. In particular, the boundary conditions imposed to the material and the location of sensors play a major role. In [30], the concept of searching the Optimal Experiment Design (OED) was used to determine the best experimental conditions in terms of imposed flux, quantity and location of sensors, aiming at estimating the thermophysical and hygrothermal properties of a material. The OED provides the best accuracy of the identification method and, thus, the estimated parameters.

Therefore, the main issue in this paper is to use the method to determine the OED considering the experimental set-up before starting the data acquisition process. First, the optimal boundary conditions and location of sensors are defined. Then, the experimental campaign is carried out, respecting the OED. From the experimental data, the parameters are estimated using an interior-point algorithm with constraints on the unknowns. To improve the fidelity of the physical phenomenon, the OED approach is extended to a model that takes into account the hysteresis effects on the sorption curves. Lastly, a comparison between the numerical predictions and the experimental observations reveals smaller discrepancies and a satisfactory representation of the physical phenomena.

This article is organized as follows. Section 2 presents the physical problem with its mathematical formulation and the OED methodology. In Section 3, the existing experimental set-up is described. The OED search providing the different possibilities for the estimation of one or several parameters with the experimental set-up is presented in Section 4. The parameters are estimated in Section 5 and then the main conclusions are finally outlined.

2. Methodology

2.1. Physical problem

In this section, the physical problem of moisture transport in porous material is described. The system of Equation (1(1) $\frac{\mathrm{\partial }U}{\mathrm{\partial }t}=\mathbf{\nabla }\mathbf{\cdot }\left(a_m\mathbf{\nabla }U+\delta a_m\mathbf{\nabla }T\right),$(1) ) proposed by Luikov has been established for non-isothermal conditions. In our study, the temperature variations in the material are assumed negligible. Therefore, in Section 2.1.1, the derivation of the unidimensional advective–diffusive moisture transport equation is detailed. Then, the assumption on the material properties is defined. In the last Section 2.1.3, the boundary conditions of the problem are specified.

2.1.1. Moisture transport

The term moisture includes the water vapour, denoted by index 1 and the liquid water, denoted by index 2, both migrating through the porous matrix of the material. The differential equation describing the mass conservation for each species can be formulated as [3]: (3) $\frac{\mathrm{\partial }{w}_i}{\mathrm{\partial }t}=-\mathbf{\nabla }\mathbf{\cdot }{\mathit{j}}_i+I_i,i\in \left\{1,2\right\},$(3) where $w_i$ is the volumetric concentration of species i, ${\mathit{j}}_i$, the total flux of the species i and $I_i$, the volumetric source. It is assumed that solid water is not present in the porous structure and $I_1+I_2=0$. Thus, by summing Equation (2(3) $\frac{\mathrm{\partial }{w}_i}{\mathrm{\partial }t}=-\mathbf{\nabla }\mathbf{\cdot }{\mathit{j}}_i+I_i,i\in \left\{1,2\right\},$(3) ) for $i\in \{1,2\}$, we obtain: (4) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(w_1+w_2\right)=-\mathbf{\nabla }\mathbf{\cdot }\left({\mathit{j}}_1+{\mathit{j}}_2\right).$(4) The total volumetric concentration of moisture is denoted by $w\stackrel{\mathrm{def}}{:=}{w}_1+w_2.$ It can be noted that we have the following relation $w={\rho }_{0}U$ between the potentials w and U from Equation (1b(2) $c_b{\rho }_0\frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\mathbf{\nabla }\mathbf{\cdot }\left(\lambda \mathbf{\nabla }T\right)+r_{12}\mathbf{\nabla }\mathbf{\cdot }\left(a_{m1}{\rho }_0\left(\mathbf{\nabla }U+{\delta }_1\mathbf{\nabla }T\right)\right),$(2) ) derived by Luikov. Since the temperature variations in the material are assumed negligible, we can write: $\frac{\mathrm{\partial }w}{\mathrm{\partial }t}=\frac{\mathrm{\partial }w}{\mathrm{\partial }\phi }\cdot \frac{\mathrm{\partial }\phi }{\mathrm{\partial }{P}_v}\cdot \frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }t},$

where φ is the relative humidity. By its definition $\phi =\frac{P_v}{P_s},$ thus, we have: $\frac{\mathrm{\partial }\phi }{\mathrm{\partial }{P}_v}=\frac{1}{P_s},$ with $P_s$ the saturation pressure. Considering the moisture sorption curve describing the material property between the moisture content w and the relative humidity φ, denoted as $w=\text{f}\left(\phi \right)$, it can be written: (5) $\frac{\mathrm{\partial }w}{\mathrm{\partial }t}=\frac{\mathrm{\partial }\text{f}\left(\phi \right)}{\mathrm{\partial }\phi }\frac{1}{P_s}\frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }t}.$(5) We define the moisture storage coefficient as (6) $c\stackrel{\mathrm{def}}{:=}\frac{\mathrm{\partial }\text{f}\left(\phi \right)}{\mathrm{\partial }\phi }\frac{1}{P_s}.$(6) Moreover, the moisture transfer occurs due to capillary migration, moisture diffusion and advection of the vapour phase. Here, advective transfer represents the movement of species under the presence of airflow through the porous matrix. The convection process designates both diffusion and advection transfer. Thus, the fluxes can be expressed as [15,31,32]: (7) ${\mathit{j}}_1+{\mathit{j}}_2=-d\mathbf{\nabla }{P}_v+\mathit{a}{P}_v,$(7) where $P_v$ is the vapour pressure, d is the global moisture transport coefficient and a the global advection coefficient defined as: $a\stackrel{\mathrm{def}}{:=}\frac{\mathsf{v}}{R_{v}T},$ where T is the fixed temperature, $\mathsf{v}$ the constant mass average velocity and $R_v$ the water vapour gas constant.

Therefore, using Equation (3(4) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(w_1+w_2\right)=-\mathbf{\nabla }\mathbf{\cdot }\left({\mathit{j}}_1+{\mathit{j}}_2\right).$(4) ) and the results from Equations (4(5) $\frac{\mathrm{\partial }w}{\mathrm{\partial }t}=\frac{\mathrm{\partial }\text{f}\left(\phi \right)}{\mathrm{\partial }\phi }\frac{1}{P_s}\frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }t}.$(5) ) and (6(7) ${\mathit{j}}_1+{\mathit{j}}_2=-d\mathbf{\nabla }{P}_v+\mathit{a}{P}_v,$(7) ), the physical problem of unidimensional advective–diffusive moisture transport through a porous material can be mathematically described as: (8) $c\frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }t}=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left[d\frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }x}\right]-a\frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }x}.$(8)

2.1.2. Material properties

The moisture capacity c is assumed to be a second-degree polynomial of the relative humidity, while the moisture permeability d is considered as a first-degree polynomial of the relative humidity: (9) $c=c_0+c_1\phi +c_2{\phi }^2,$(9) (10) $d=d_0+d_1\phi .$(10)

These assumptions correspond to a third-order polynomial function of the relative humidity for the moisture sorption curve. It should be noted that other functions can be used to describe the moisture sorption curve c and/or the moisture permeability d. However, for the material under investigation, these properties have been determined using the traditional cup method and the gravimetric approach, presented in preliminary studies [33,34] and expressed using these functions. These properties will be used as a priori ones in the algorithm when searching the OED in Section 4. It is important to stress out that the hysteresis of the sorption curve is not considered in the physical model.

2.1.3. Boundary conditions

At $x=0$, the surface is in contact with the ambient air at temperature $T^{\mathrm{\infty }}$ and relative humidity ${\phi }^{\mathrm{\infty }}$. Thus, the boundary condition is expressed as: (11) $d\frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }x}-a{P}_v=h(P_v-P_s\left(T^{\mathrm{\infty }}\right){\phi }^{\mathrm{\infty }}\left(t\right)),$(11) where h is the convective vapour transfer coefficient, considered as constant. At $x=L$, the surface is impermeable. Thus, the total moisture flow vanishes at this boundary where the velocity $\mathsf{v}$ and the diffusion flow $d\left(\mathrm{\partial }{P}_v/\mathrm{\partial }x\right)$ are null: (12) $d\frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }x}-a{P}_v=0.$(12) At t=0, the vapour pressure is supposed to be uniform within the material (13) $P_v=P_v^i.$(13) It is of great concern in the construction of the numerical model that the boundary conditions satisfies the initial condition.

2.2. Dimensionless formulation

For building porous materials such as concrete, insulation and brick, the coefficients scale as ${10}^2\text{kg}/({\text{m}}^3.\text{Pa})$ for the sorption curve c and ${10}^{-11}\text{s}$ for the moisture permeability d and the advection coefficients a. Therefore, while performing a mathematical and numerical analysis of a given practical problem, it is of capital importance to obtain a unitless formulation of governing equations. For this, the vapour pressure is transformed to a dimensionless quantity: $u=\frac{P_v}{P_v^{\mathrm{ref}}},u_{\mathrm{i}}=\frac{P_v^i}{P_v^{\mathrm{ref}}},u_{\mathrm{\infty }}=\frac{P_s\left(T^{\mathrm{\infty }}\right){\phi }^{\mathrm{\infty }}}{P_v^{\mathrm{ref}}}.$

The time and space domains are also modified: $x^{\star }=\frac{x}{L},t^{\star }=\frac{t}{t^{\mathrm{ref}}},$ where L is the thickness of the material and $t^{\mathrm{ref}}$ a characteristic time. The material properties are changed considering a reference value for each parameter: $\begin{array}{r}{c}^{\star }=\frac{c}{c_0},d^{\star }=\frac{d}{d_0}.\end{array}$ In this way, dimensionless numbers are highlighted: $\begin{array}{r}\text{P'e}=\frac{a\cdot L}{d_0},\text{Bi}=\frac{h\cdot L}{d_0},\text{Fo}=\frac{t^{\mathrm{ref}}\cdot d_0}{L^2\cdot c_0}.\end{array}$

The dimensionless moisture Biot number $\text{Bi}$ quantifies the importance of the moisture transfer at the bounding surface of the material. The transient transfer mechanism is characterized by the Fourier number $\text{Fo}$, whereas the Péclet number $\text{P'e}$ measures only the importance of moisture advection. The quantities $c^{\star }\left(u\right)$ and $d^{\star }\left(u\right)$ give the variation of storage and permeability coefficients from the reference state of the material. The dimensionless governing equations are finally written as: (14) $\begin{array}{rlrl}{c}^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}& =\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}^{\star }}\left(d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}-\text{P'e}u\right),& t^{\star }& >0,x^{\star }\in \left[0,1\right],\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=\text{Bi}\cdot \left(u-u_{\mathrm{\infty }}\right),& t^{\star }& >0,x^{\star }=0,\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=0,& t^{\star }& >0,x^{\star }=1,\\ u& =u_{\mathrm{i}},& t^{\star }& =0,x^{\star }\in \left[0,1\right],\end{array}$(14) where functions $c^{\star }\left(u\right)$ and $d^{\star }\left(u\right)$ are given by: $\begin{array}{rl}{c}^{\star }\left(u\right)& =1+c_1^{\star }u+c_2^{\star }{u}^2,\\ d^{\star }\left(u\right)& =1+d_1^{\star }u.\end{array}$

In this study, the reference parameters correspond to the conditions at $T={23}^{\circ }C$ and $\phi =0.5$, which gives a vapour pressure of $P_v^{\mathrm{ref}}=1404\text{Pa}$. Since the temperature is assumed as constant, a dimensionless value $u=2$ corresponds to $\phi =1$. The field varies within the interval $u\in [0,2]$.

The direct problem, defined by Equation (12(14) $\begin{array}{rlrl}{c}^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}& =\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}^{\star }}\left(d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}-\text{P'e}u\right),& t^{\star }& >0,x^{\star }\in \left[0,1\right],\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=\text{Bi}\cdot \left(u-u_{\mathrm{\infty }}\right),& t^{\star }& >0,x^{\star }=0,\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=0,& t^{\star }& >0,x^{\star }=1,\\ u& =u_{\mathrm{i}},& t^{\star }& =0,x^{\star }\in \left[0,1\right],\end{array}$(14) ), is solved using a finite-difference standard discretisation method. An embedded adaptive in time Runge–Kutta scheme combined with a Scharfetter–Gummel spatial discretisation approach is used [15]. It is adaptive and embedded to estimate local error in time with low extra cost. The algorithm was implemented in the Matlab^TM environment. For the sake of notation compactness, the upper-script ⋆ standing for dimensionless parameters is no longer used.

2.3. The Optimal Experiment Design

Efficient computational algorithms for recovering parameters $\mathbf{P}$ given an observation $u_{\exp }$ of the field $u(x,t)$ have already been proposed. Readers may refer to [35] for a primary overview of different methods. They are based on the minimization of the cost function $\text{J}\left[\mathbf{P}\right]$. For this purpose, it is required to equate to zero the derivatives of $\text{J}\left[\mathbf{P}\right]$, with respect to each of the unknown parameters $p_m$ to find critical points. Associated to this necessary condition for the minimization of $\text{J}\left[\mathbf{P}\right]$, the scaled dimensionless local sensitivity function [36,37] is introduced: (15) ${\mathrm{\Theta }}_m(x,t)=\frac{{\sigma }_p}{{\sigma }_u}\frac{\mathrm{\partial }u}{\mathrm{\partial }{p}_m},\mathrm{\forall }m\in \left\{1,\dots ,M\right\},$(15) where ${\sigma }_u$ is the variance of the error measuring $u_{\exp }$. The parameter scaling factor ${\sigma }_p$ equals to 1 as we consider that prior information on parameter $p_m$ has low accuracy. It is important to note that all algorithms have been developed considering the dimensionless problem in order to compare only the order of variation of parameters and observation, avoiding the effects of units and scales.

The function ${\mathrm{\Theta }}_m$ measures the sensitivity of the estimated field u with respect to changes in the parameter $p_m$ [35,38,39]. A small magnitude of ${\mathrm{\Theta }}_m$ indicates that large changes in $p_m$ induce small changes in u. The estimation of parameter $p_m$ is therefore difficult, in this case. When the sensitivity coefficient ${\mathrm{\Theta }}_m$ is small, the inverse problem is necessarily ill-conditioned. If the sensitivity coefficients are linearly dependent, the inverse problem is also ill-posed. Therefore, to get an optimal evaluation of parameters $\mathbf{P}$, it is desirable to have linearly independent sensitivity functions ${\mathrm{\Theta }}_m$ with large magnitudes for all parameters $p_m$. These requirements ensure the best conditioning of the computational algorithm to solve the inverse problem and thus the better accuracy of the estimated parameter.

It is possible to define the experimental design in order to meet these requirements. The issue is to find the optimal sensor location $X^{\circ }$ and the optimal amplitude ${\phi }^{\mathrm{\infty },\circ }$ of the relative humidity of the ambient air at the material bounding surface, $x=0$. To search this OED, we introduce the following measurement plan: (16) $\pi \stackrel{\mathrm{def}}{:=}\left\{X,{\phi }^{\mathrm{\infty }}\right\}.$(16)

In the analysis of optimal experiments for estimating the unknown parameter(s) $\mathbf{P}$, a quality index describing the recovering accuracy is the $D-$optimum criterion [40–45]: (17) $\mathrm{\Psi }=det\left[F\left(\pi \right)\right],$(17) where $F\left(\pi \right)$ is the normalized Fisher information matrix [46,47] defined as: (18) $F\left(\pi \right)=\left[{\mathrm{\Phi }}_{ij}\right],\mathrm{\forall }(i,j)\in {\left\{1,\dots ,M\right\}}^2,$(18) (19) ${\mathrm{\Phi }}_{ij}=\sum _{n=1}^N{\int }_0^{\tau }{\mathrm{\Theta }}_i(x_n,t){\mathrm{\Theta }}_j(x_n,t)\text{d}t.$(19)

The matrix $F\left(\pi \right)$ characterizes the total sensitivity of the system as a function of the measurement plan π (Equation (14(16) $\pi \stackrel{\mathrm{def}}{:=}\left\{X,{\phi }^{\mathrm{\infty }}\right\}.$(16) )). The OED search aims at finding a measurement plan ${\pi }^{\star }$ for which the objective function (Equation (15(17) $\mathrm{\Psi }=det\left[F\left(\pi \right)\right],$(17) )) reaches the maximum value: (20) ${\pi }^{\circ }=\left\{X^{\circ },{\phi }^{\mathrm{\infty },\circ }\right\}=\arg \underset{\pi }{max}\mathrm{\Psi }.$(20)

To solve Equation (17(20) ${\pi }^{\circ }=\left\{X^{\circ },{\phi }^{\mathrm{\infty },\circ }\right\}=\arg \underset{\pi }{max}\mathrm{\Psi }.$(20) ), a domain of variation ${\mathrm{\Omega }}_{\pi }$ is considered for the sensor position X and the amplitude ${\phi }^{\mathrm{\infty }}$ of the boundary conditions. Then, the following steps are carried out for each value of the measurement plan $\pi =\{X,{\phi }^{\mathrm{\infty }}\}$ in the domain ${\mathrm{\Omega }}_{\pi }$. First, the direct problem, defined by Equations (7(8) $c\frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }t}=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left[d\frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }x}\right]-a\frac{\mathrm{\partial }{P}_v}{\mathrm{\partial }x}.$(8) )–(11(13) $P_v=P_v^i.$(13) ), is solved. Then, given the solution $P_v$ for a fixed value of the measurement plan, the next step consists of computing the sensitivity coefficients ${\mathrm{\Theta }}_m=\left(\mathrm{\partial }u/\mathrm{\partial }{p}_m\right)$, using also an embedded adaptive time Runge–Kutta scheme combined with central spatial discretisation. Then, with the sensitivity coefficients, the Fisher matrix (16(18) $F\left(\pi \right)=\left[{\mathrm{\Phi }}_{ij}\right],\mathrm{\forall }(i,j)\in {\left\{1,\dots ,M\right\}}^2,$(18) )(a,b) and the $D-$optimum criterion (15(17) $\mathrm{\Psi }=det\left[F\left(\pi \right)\right],$(17) ) are computed. The solution of the direct and sensitivity problems is obtained for a given a priori parameter $\mathbf{P}$ and, in this case, the validity of the OED depends on this knowledge. If there is no prior information, the methodology of the OED can be done using an outer loop on the parameter $\mathbf{P}$ sampled using, for instance, Latin hypercube or Halton or Sobol quasi-random samplings. Interested readers may refer to [30] for further details on the computation of sensitivity coefficients.

An interesting remark with this approach is that the probability distribution of the unknown parameter $p_m$ can be estimated from the distribution of the measurements of the field u and from the sensitivity ${\mathrm{\Theta }}_m$. The probability $\mathcal{P}$ of u is given by: $F\left(\overline{u}\right)=\mathcal{P}\left\{u(x,t,p_m)⩽\overline{u}\right\}.$

Using the sensitivity function ${\mathrm{\Theta }}_m$, the probability can be approximated by: $\begin{array}{rl}F\left(\overline{u}\right)& \simeq \mathcal{P}\left\{u(x,t,p_m^{\circ })+{\mathrm{\Theta }}_m\cdot \left(p_m-p_m^{\circ }\right)⩽\overline{u}\right\}.\end{array}$ Assuming ${\mathrm{\Theta }}_m>0$, we get: $F\left(\overline{u}\right)=\mathcal{P}\left\{p_m⩽p_m^{\circ }+\frac{\overline{u}-u(x,t,p_m^{\circ })}{{\mathrm{\Theta }}_m}\right\}.$

Therefore, using a change of variable, the cumulative derivative function of the probability of the unknown parameter $p_m$ is estimated by: $\begin{array}{rl}F\left({\overline{p}}_m\right)& =\mathcal{P}\left\{p_m⩽{\overline{p}}_m\right\}\\ & =F\left(u+{\mathrm{\Theta }}_m\cdot \left({\overline{p}}_m-p_m^{\circ }\right)\right).\end{array}$ When ${\mathrm{\Theta }}_m<0$, the cumulative derivative function of the probability is given by: $F\left({\overline{p}}_m\right)=1-F\left(u+{\mathrm{\Theta }}_m\cdot \left({\overline{p}}_m-p_m^{\circ }\right)\right).$

It gives a local approximation of the probability distribution of the unknown parameter $p_m$, at a reduced computational cost. Moreover, the approximation is reversible. Thus, if one has the distribution of the unknown parameter, it is possible to get the one of field u.

3. Experimental facility

The test facility used to carry out the experiment is illustrated in Figure 2. It is composed of two connected climatic chambers. The temperature of each chamber is controlled independently with a thermostatically controlled water bath allowing water to recirculate in a heat exchanger. The relative humidity is kept fixed using saturated salt solutions of ${\text{MgCl}}_2$ and $\text{NaCl}$. Relative humidity values in chambers 1 and 2 are fixed to ${\phi }_1=33\mathrm{\%}$ and ${\phi }_2=75\mathrm{\%}$, respectively. Two door locks, at each side, allow the operator to insert or remove samples to minimize system disturbances. They enable easy and instantaneous change in humidity boundary conditions for the samples while passing from one chamber to another. Another climatic chamber is also available used to initially condition materials at ${\phi }_0=10\mathrm{\%}$.

Figure 2. Illustration of the RH-box experimental facility.

The temperature and relative humidity fields are measured within the samples with wireless sensors from the HygroPuce range (Waranet industry). The accuracy is $\pm 2\mathrm{\%}$ for the relative humidity and $\pm 0.5{}^{\circ }\mathrm{C}$ for the temperature and the dimensions are $0.6\text{cm}$ thickness and $1.6\text{cm}$ diameter, as illustrated in Figure 3(a). The sensors are placed within the material by cutting the samples. The total uncertainty on the measurement of relative humidity can be evaluated considering the propagation of the uncertainty due to sensor measurement and due to their location. In [23], the total uncertainty on the measurement has been evaluated to $\mathrm{\Delta }\phi =2\mathrm{\%}$.

Figure 3. Sensors of relative humidity and temperature (a) and wood fibre samples (b) with white acrylic seal and with aluminium tape.

The material investigated is the wood fibre, which properties have been determined in [23,33] and are shown in its dimensionless form in Table 1. The reference parameter used to compute the unitless parameters is $t^{\mathrm{ref}}=3600\text{s}$, $d_0=5.17\text{s}$, $c_0=2.85\text{kg}/({\text{m}}^3.\text{Pa})$ and $L=0.08\text{m}$. It constitutes a priori information on the unknown parameters $\text{Fo}$, $c_1$ and $c_2$. The samples are cylindrical, with a $10\text{cm}$ diameter and $8\text{cm}$ thickness in order to avoid border effects and to minimize perturbations by sensors placed within the sample. Moreover, to ensure unidimensional moisture transfer and a null flux condition at x=1, the samples are covered with aluminium tape and glued on a white acrylic seal, as illustrated in Figure 3(b). The convective moisture transport coefficient at x=0 has been estimated experimentally in [23,48]. The corresponding Biot number is reported in Table 1.

Table 1. A priori dimensionless material properties of wood fibre from [23,33].

Storage coefficient c	$c\left(u\right)=c_0+c_{1}u+c_2{u}^2$
	$c_0=1$, $c_1=-9.79\cdot {10}^{-1}$, $c_2=1.06$
Diffusion coefficient d	$d\left(u\right)=d_0+d_{1}u$
	$d_0=1$, $d_1=0.29$
PÉCLET number	$\text{P'e}=1.1\cdot {10}^{-2}$
FOURIER number	$\text{Fo}=4\cdot {10}^{-3}$
BIOT number	$\text{Bi}=13.7$

Finally, the experimental facility is used to submit the samples to a single or multiple steps of relative humidity. For a single step, the boundary conditions are defined as: $\begin{array}{r}{u}_{\mathrm{\infty }}\left(t\right)=\left\{\begin{array}{ll}{u}_{\mathrm{i}},& t=0,\\ u_{\mathrm{c}},& t>0.\end{array}\right.\end{array}$ For the case of multiple steps, we set: $\begin{array}{r}{u}_{\mathrm{\infty }}\left(t\right)=\left\{\begin{array}{ll}{u}_{\mathrm{i}},& t=0,\\ u_{\mathrm{c}}^1,& t\in \left(0,\tau \right],\\ u_{\mathrm{c}}^2,& t\in \left(\tau ,2\tau \right],\\ u_{\mathrm{c}}^3,& t\in \left(2\tau ,3\tau \right],\end{array}\right.\end{array}$ where the initial condition belongs to $u_{\mathrm{i}}\in \{0.2,0.6,1.5\}$, the climatic chamber boundary condition $(u_{\mathrm{c}},u_{\mathrm{c}}^1,u_{\mathrm{c}}^2,u_{\mathrm{c}}^3)\in \{0.6,1.5\}$ and the duration of the step $\tau \in \{24,48,72,96,120,144,168,192\}$. A total of 20 designs are possible for providing measurements to estimate the unknown parameters $\text{Fo}$, $c_1$ and $c_2$. A synthesis of the possible designs is provided in Table 2. It should be noted that, according to the reference parameters, unitless values u=0.2, u=0.66, u=1.5 correspond to $\phi =0.1$, $\phi =0.33$ and $\phi =0.75$, respectively.

Table 2. Possible designs according to the experimental facility.

			Multiple step
Single step			Case 1				Case 2
Design	$u_{\mathrm{i}}$	$u_{\mathrm{c}}$	Design	τ	Design	τ	Design	τ	Design	τ
1	0.2	0.66	5	1	9	5	13	1	17	5
2	0.2	1.5	6	2	10	6	14	2	18	6
3	0.66	1.5	7	3	11	7	15	3	19	7
4	1.5	0.66	8	4	12	8	16	4	20	8

Case 1: $u_{\mathrm{i}}=0.2$, $u_{\mathrm{c}}^1=0.66$, $u_{\mathrm{c}}^2=1.5$, $u_{\mathrm{c}}^3=0.66$.

Case 2: $u_{\mathrm{i}}=0.2$, $u_{\mathrm{c}}^1=1.5$, $u_{\mathrm{c}}^2=0.66$, $u_{\mathrm{c}}^3=1.5$.

4. Searching the OED

4.1. Estimation of one parameter

This section focuses on the estimation of one parameter within $\text{Fo}$, $c_1$ or $c_2$ with experiments coming from single- or multiple-step designs. It should be noted that by estimating parameter $\text{Fo}$, the complete sorption isotherm curve is defined, according to the dimensionless quantities defined in Section 2.2. The equations to compute the sensitivity functions are given in Appendix 1. In addition, the demonstration of structural identifiability of the three parameters is provided in Appendix 2.

4.1.1. Single step

Figure 4(a) gives the variation of the criterion Ψ for the four single-step designs. For the estimation of parameters $\text{Fo}$ or $c_1$, the criterion reaches its maximal value for design 2, corresponding to a step from $u_{\mathrm{i}}=0.1$ to $u_{\mathrm{c}}=1.5$. For parameter $c_2$, the design 4 is the optimal one. It can be noted that, for parameter $c_1$, the relative criterion Ψ attains $80\mathrm{\%}$ for the design 4. It could be an interesting alternative to estimate this parameter. The variation of the criterion is related to the sensitivity function of each parameter. As noticed in Figure 5(a,b), functions Θ have higher magnitudes of variation for the OED. The variation of Ψ as a function of the sensor location X is shown in Figure 4(b) for the OED. The optimal sensor location is at the boundary of the material opposite from the perturbations. If required for practical purpose, the sensor can be placed in the interval $X\in [0.9,1]$, ensuring to reach $95\mathrm{\%}$ of the criterion Ψ. Results have similar tendencies for the three parameters.

Figure 4. Variation of the criterion Ψ for the four possible single-step designs (a) and as a function of the sensor position X for the OED (b), in the case of estimating one parameter.

Figure 5. Sensitivity coefficients Θ for parameters $\text{Fo}$, $c_1$ and $c_2$ for the OED (a) and for design 1 (b) ($X=X^{\circ }$).

4.1.2. Multiple steps

Figure 9(a) shows the variation of the relative criterion Ψ for the designs considering multiple steps of relative humidity. It increases with the duration τ of the steps. Thus, for the group of designs 5–12 and the groups 13 –20, the criteria reach their maximum for designs 12 and 20, respectively, corresponding to the step duration $\tau =8$. The groups 5–12 correspond to a multiple step $u_{\mathrm{c}}^1=0.66$, $u_{\mathrm{c}}^2=1.5$, $u_{\mathrm{c}}^3=0.66$. For them, the criterion does not attain $80\mathrm{\%}$ of the maximal criteria. Therefore, it is preferable to choose among designs 18–20, with a multiple step $u_{\mathrm{c}}^1=1.5$, $u_{\mathrm{c}}^2=0.66$, $u_{\mathrm{c}}^3=1.5$, and a duration $\tau ⩾6$ to estimate the parameters. Figure 7(a–d) compares the sensitivity function of each parameter for three different designs. The quantity Θ has higher magnitude of variation for the OED than for the others. Moreover, for the design 5, the duration of the step is so short that there is almost no variation on the sensitivity when occurring the first step for $t\in [0,24]$. As for the previous case, the optimal sensor position is $X\in [0.9,1]$.

Figure 6. Variation of the criterion Ψ for the 16 possible designs (a) and as a function of the sensor position X for the OED (b), in the case of estimating one parameter.

Figure 7. Sensitivity coefficients Θ for parameters $\text{Fo}$, $c_1$ and $c_2$ for the OED (design 20) (a), design 12 (b), design 5 (c) and design 13 (d) ($X=X^{\circ }$).

Figure 8. Variation of the criterion Ψ for the four possible designs of single step of relative humidity (a) and as a function of the sensor position X for the OED (b), in the case of estimating the couple of parameters $(\text{Fo},c_2)$.

Figure 9. Variation of the criterion Ψ for the 16 possible designs of multiple steps of relative humidity (a) and as a function of the sensor position X for the OED (b), in the case of estimating the couple of parameters $(\text{Fo},c_2)$.

4.2. Estimation of several parameters

The issue is now to estimate two or three parameters defining the moisture capacity $\text{Fo}$, $c_1$ and $c_2$. First of all, it is important to notice in Figures 5 and 7 that the sensitivity function Θ of the parameters has a strong correlation. The interval of variation of the correlation coefficients for all the designs is: $\begin{array}{rl}\mathrm{Cor}\left(\text{Fo},c_1\right)& \in \left[0.94,0.99\right],\\ \mathrm{Cor}\left(c_1,c_2\right)& \in \left[0.92,0.99\right],\\ \mathrm{Cor}\left(\text{Fo},c_2\right)& \in \left[0.71,0.95\right].\end{array}$ Therefore, the estimation of the three parameters at the same time using only one experiment might be a difficult task. In addition, over all the possible designs, the couple of parameters $(\text{Fo},c_2)$ is the one with the lower correlation. Therefore, the OED search will only consider their estimation.

4.2.1. Single step

Figure 8(a) gives the variation of the criterion Ψ for the four possible designs, considering a single step of relative humidity. The OED is reached for design 4. However, the design 2 represents an interesting alternative as more than $95\mathrm{\%}$ of the maximum criterion is reached. The sensor should be placed between $X\in [0.9,1]$.

4.2.2. Multiple steps

The variation of the criterion Ψ for the 16 designs is shown in Figure 9. It increases with the duration of the steps τ. The OED is reached for the design considering a multiple step $u_{\mathrm{i}}=0.2$, $u_{\mathrm{c}}^1=1.5$, $u_{\mathrm{c}}^2=0.66$ and $u_{\mathrm{c}}^3=1.5$, with a duration $\tau =8$. As for the previous case, the OED is defined for a sensor placed near the boundary of the material $x=1$.

5. Estimation of the unknown parameters

5.1. Methodology

As mentioned before, the sensitivity functions of parameters $\text{Fo}$, $c_1$ and $c_2$ are strongly correlated and the estimation of the three parameters using one single experiment might be a laborious task. To answer this issue, a single step, referenced as experiment A, will be used for the estimation of the parameter $c_1$ and a multiple step, referenced as experiment B, for the parameters $(\text{Fo},c_2)$, which sensitivity functions are less correlated. According to the results of the OED, the sensor is placed near the border $x=1$. For the boundary conditions, the single step will be operated from $u_{\mathrm{i}}=0.2$ to $u_{\mathrm{c}}=1.5$ (design 2 from Table 2). The OED multiple step is defined as $u_{\mathrm{i}}=0.2$, $u_{\mathrm{c}}^1=1.5$, $u_{\mathrm{c}}^2=0.66$, $u_{\mathrm{c}}^3=1.5$ and a duration of each step $\tau =8$ (design 20 from Table 2).

To estimate the unknown parameters, the following cost function is defined by minimizing the residual between the experimental data and the numerical results of the direct model: ${\text{J}}_i^n={\left|\left|u-u_{\exp ,n}\right|\right|}_{{\mathcal{L}}_i},n\in \left\{1,2,3\right\},i\in \left\{2,\mathrm{\infty }\right\}.$

Several expressions of the cost function are tested. The subscript i denotes the norm used for the cost function $\text{J}$: i=2 stands for the standard discrete ${\mathcal{L}}_2$ norm while $i=\mathrm{\infty }$ for the ${\mathcal{L}}_{\mathrm{\infty }}$ (uniform) norm. The upper-script $n=1$ implies that both experiments are considered separately. The single step is used for estimating parameter $c_1$ and the multiple-step experiment for the parameters $(\text{Fo},c_2)$. In such case, there is a cost function according to each experiment: $\begin{array}{rl}{\text{J}}_i^{1A}\left(\text{Fo},c_2\right)& ={\left|\left|u-u_{\exp ,A}\right|\right|}_{{\mathcal{L}}_i},\mathrm{for}\mathrm{the}\mathrm{single}-\mathrm{step}\mathrm{experiment},\\ {\text{J}}_i^{1B}\left(c_1\right)& ={\left|\left|u-u_{\exp ,B}\right|\right|}_{{\mathcal{L}}_i},\mathrm{for}\mathrm{the}\mathrm{multiple}-\mathrm{step}\mathrm{experiment}.\end{array}$

The estimation of the unknown parameters proceeds in an iterative approach as described in the Algorithm 1. In this case, a tolerance $\eta ={10}^{-6}$ has been chosen.

When $n\in \{2,3\}$, both experiments of single and multiple steps are considered at the same time, without any distinction. For $n=2$, parameters $(\text{Fo},c_1,c_2)$ are estimated at once. An additional test, for $n=3$ is carried by considering both experiments to estimate all the material properties parameters $(\text{Fo},c_1,c_2,d_1,\text{P'e})$. Thus, in these cases, the cost functions are defined as: $\begin{array}{rl}{\text{J}}_2^2\left(\text{Fo},c_1,c_2\right)& ={\left|\left|u-u_{\exp ,A}\right|\right|}_{{\mathcal{L}}_2}+{\left|\left|u-u_{\exp ,B}\right|\right|}_{{\mathcal{L}}_2},\\ {\text{J}}_{\mathrm{\infty },\max }^2\left(\text{Fo},c_1,c_2\right)& =max\left\{{\left|\left|u-u_{\exp ,A}\right|\right|}_{{\mathcal{L}}_{\mathrm{\infty }}},{\left|\left|u-u_{\exp ,B}\right|\right|}_{{\mathcal{L}}_{\mathrm{\infty }}}\right\},\\ {\text{J}}_{\mathrm{\infty },+}^2\left(\text{Fo},c_1,c_2\right)& ={\left|\left|u-u_{\exp ,A}\right|\right|}_{{\mathcal{L}}_{\mathrm{\infty }}}+{\left|\left|u-u_{\exp ,B}\right|\right|}_{{\mathcal{L}}_{\mathrm{\infty }}},\end{array}$ for $n=2$, and for $n=3$: $\begin{array}{rl}{\text{J}}_2^3\left(\text{Fo},c_1,c_2,d_1,\text{P'e}\right)& ={\left|\left|u-u_{\exp ,A}\right|\right|}_{{\mathcal{L}}_2}+{\left|\left|u-u_{\exp ,B}\right|\right|}_{{\mathcal{L}}_2},\\ {\text{J}}_{\mathrm{\infty },\max }^3\left(\text{Fo},c_1,c_2,d_1,\text{P'e}\right)& =max\left\{{\left|\left|u-u_{\exp ,A}\right|\right|}_{{\mathcal{L}}_{\mathrm{\infty }}},{\left|\left|u-u_{\exp ,B}\right|\right|}_{{\mathcal{L}}_{\mathrm{\infty }}}\right\},\\ {\text{J}}_{\mathrm{\infty },+}^3\left(\text{Fo},c_1,c_2,d_1,\text{P'e}\right)& ={\left|\left|u-u_{\exp ,A}\right|\right|}_{{\mathcal{L}}_{\mathrm{\infty }}}+{\left|\left|u-u_{\exp ,B}\right|\right|}_{{\mathcal{L}}_{\mathrm{\infty }}}.\end{array}$

Table 3 synthesizes all tests performed according to the definition of the cost function $\text{J}$. The cost function $\text{J}$ is minimized using function fmincon in the Matlab^TM environment, providing an efficient interior-point algorithm with constraint on the unknown parameters [49]. Here, the box-type constraints are defined with upper and lower bound for the parameters: $p^{\circ }\in \left[0.8,1.5\right]\cdot p^{\mathrm{apr}},$ where the upper-scripts ° and $\mathrm{apr}$ denote the estimated and a priori values of the parameters, respectively. The bounds have been defined by performing previous tests and analysing the parameter impact on the calibration.

Table 3. Synthesis of the tests carried out with the expression of the cost function.

	Tests
Definition of the cost function $\text{J}$	1	2	3	4	5	6	7	8	9
Experiments considered separately	x	x	x
Experiments at the same time (3 parameters to estimate)				x	x	x
Experiments at the same time (5 parameters to estimate)							x	x	x
Euclidean norm	x			x			x
Infinite norm		x			x			x
Sum of the infinite norms			x			x			x

In order to quantify the quality of measured data, we estimate the noise inherent to any real physical measurement. By assuming that the noise $\xi \left(\omega \right)$ is centred Gaussian (i.e. $\xi \left(\omega \right)\sim \mathcal{N}(0,{\sigma }^2)$), linear and additive, its variance ${\sigma }^2$ can be thus estimated. Moreover, we assume that the underlying signal is smooth. In order to extract the noise component, the signal is approximated locally by a low-order polynomial representing the trend. Then, the trend is removed by using a special filter, leaving us with the pure noise content, which can be further analysed using the standard statistical techniques. For the considered data, the variance equals: $\begin{array}{rl}& \sigma \simeq 0.01,\mathrm{for}\mathrm{the}\mathrm{single}-\mathrm{step}\mathrm{experiments},\\ & \sigma \simeq 0.008,\mathrm{for}\mathrm{the}\mathrm{multiple}-\mathrm{step}\mathrm{experiments}.\end{array}$ The noise variance does not necessarily correspond to the measurement accuracy. This measure provides a lower bound of this error, i.e. the accuracy cannot be lower than the noise present in the measurement.

5.2. Results

Figure 10(a) shows the variation of the residual between the measured data and the numerical results for different tests performed. The residual is minimized for tests 1, 4 and 7, corresponding to the involving considering the Euclidean norm for the computation of the cost functions. The tendencies are similar for both experiments. It can be noted that the residual is lower when estimating only three parameters $(\text{Fo},c_1,c_2)$ and not all the parameters of the material properties. Figure 10(b) provides the number of computations of the direct problem. As expected, the tests 1–3, considering both experiments separately, require a few more computations of the direct problem, due to the iterative procedure. Globally, the algorithm requires less than 100 computations, which is extremely low compared to stochastic approaches. For instance in [18], ${10}^4$ computations are necessary to estimate the thermal conductivity of two materials by solving an inverse heat transfer problem.

Figure 10. Residual between the measured data and numerical results for both experiments (a) and number of computations of the direct problem (b) for the different definition of the cost function $\text{J}$.

The comparison of the measured data and numerical results is illustrated in Figure 11(a) for the one-step experiment. Figure 11(b) shows it for the multiple-step procedure, for the case 20. Results of the numerical model are provided for the a priori and estimated three parameters. As mentioned in the Introduction, the numerical model with a priori parameters predicts values lower than those obtained from experiments during the moisture adsorption and greater than them along the desorption processes. With the calibrated model, i.e. with the estimated parameters, there is a better agreement between the numerical results and the experimental data. Figure 11(c,d) shows the residual. It is uncorrelated, highlighting a satisfactory estimation of the parameters. Nevertheless, it can be noted that some discrepancies remain between the experimental data and the numerical results. This can be specifically observed at $t=200$, in Figure 11(b,d), for which some explanations are possible. First, the mathematical model may fail in representing the physical phenomenon. Some assumptions such as isothermal conditions, unidimensional transport and constant velocity might contribute to the differences observed between experimental and numerical results. Although the experiment has been conceived to be isothermal, slight variations in the temperature field occur due to mechanisms of phase change that may affect the profile of vapour pressure, which is highly temperature-dependent. An important assumption, very often considered in literature, is the disregard of the hysteresis phenomenon, which may considerably affect the results. This issue can be noted in Figure 11(b) for $t\in [0,300]$. A good agreement is noted during the adsorption process $t\in [0,200]$, since the coefficients have been estimated for a similar experiment. However, during the desorption process, the field u is underestimated by the numerical predictions. On the other hand, despite the low relative humidity uncertainty, other uncertainties appear such as the interference of sensors on the moisture transfer through the sample, contact resistances and no perfect impermeabilization. Another possible explanation is associated to the parametrization of the material properties that can be improved. An interesting alternative could be to search for time varying parameters by adding a regularization term in the cost function $\text{J}$. The convergence of the parameter estimation is shown in Figure 12(a). On the contrary to parameters $\text{Fo}$ and $c_1$, the a priori values of $c_2$ are not far from the estimated one. After one iteration, the algorithm almost estimates the parameters. The number of computations of the direct model for test 1 is given in Figure 12(b). Only two global iterations are required to compute the solution of the inverse problem. More computations are required at the first iteration as the unknown parameters are more distant from the estimated optimal ones.

Figure 11. Comparison of the numerical results with the experimental data (a–b) and their $98\mathrm{\%}$ confidence interval for the single-step (a) and the multiple-step (b) experiments (test 1). Comparison of the residual for the single-step (c) and multiple-step (d) experiments.

Figure 12. Convergence of the parameter estimation, for the test 1 as a function of the number of iterations (a) and number of computations of the direct model (b).

The computation of the sensitivity functions of the parameters to be estimated enables to approximate their probability density functions. The error is assumed as a normal distribution $\mathcal{N}(0,{\sigma }_{\mathrm{eff}}^2)$ with its standard deviation ${\sigma }_{\mathrm{eff}}$ computed by adding the ones due to uncertainty and to the noise: ${\sigma }_{\mathrm{eff}}={\sigma }_{\mathrm{noise}}+{\sigma }_{\mathrm{error}}.$

Thus, the probability distribution is computed for different times as illustrated in Figure 13. As reported in Figure 7(a), the sensitivity function of parameters $\text{Fo}$ and $c_1$ is maximum and minimum at t=207 and $t=4$, respectively. It explains why the probability function is maximum at t=207 for these parameters. Similar results can be observed when comparing the sensitivity function of parameter $c_2$ in Figure 5(a) and the probability function in Figure 13(c).

Figure 13. Probability density function approximated for the estimated parameters, computed using the sensitivity of single-step (c) and a multiple-step (a–b) experiments.

6. Accounting for hysteresis

As mentioned in the description of the physical model, the hysteresis effects are not considered in the sorption capacity of the material. This assumption impacts the comparison between numerical and experimental observations, particularly for the cases of multiple steps of relative humidity (Figure 11(b)). In this section, the physical model is changed to revise this assumption, which is commonly considered.

6.1. Model for hysteresis

The hysteresis impacts the sorption curve $f\left(\phi \right)$ and consequently on the moisture capacity coefficient c defined in Equation (5(6) $c\stackrel{\mathrm{def}}{:=}\frac{\mathrm{\partial }\text{f}\left(\phi \right)}{\mathrm{\partial }\phi }\frac{1}{P_s}.$(6) ). To account for hysteresis, several model are available in literature as proposed in [50] or in [51]. Recently, in [52], a new model of hysteresis has been proposed and validated using experimental data. This model is referred as a smoothed bang–bang model in control literature. It enables intermediary sorption states between the main adsorption $c_{\mathrm{ads}}^{\star }$ and desorption $c_{\mathrm{des}}^{\star }$ curves illustrated in Figure 18(b). The coefficient c is computed using the following dimensionless differential equation: (21) $\frac{\mathrm{\partial }{c}^{\star }}{\mathrm{\partial }{t}^{\star }}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}\right)\left(c^{\star }-c_{\mathrm{ads}}^{\star }\right)\left(c^{\star }-c_{\mathrm{des}}^{\star }\right).$(21) As for the previous case, these two properties are assumed to be second-degree polynomials of the relative humidity: $\begin{array}{rl}{c}_{\mathrm{ads}}^{\star }& =1+c_{a,1}^{\star }u+c_{a,2}^{\star }{u}^2,\\ c_{\mathrm{des}}^{\star }& =c_{d,0}^{\star }+c_{d,1}^{\star }u+c_{d,2}^{\star }{u}^2.\end{array}$ It should be noted here that the dimensionless coefficient $c^{\star }$ is defined as $c^{\star }=c/c_{a,0}^{\star }$. The coefficient β is a numerical parameter, controlling the transition velocity between the two curves. The function $\mathrm{sign}:\mathbb{R}\to \{-1,0,1\}$ is defined as: $\begin{array}{r}\mathrm{sign}\left(x\right)=\left\{\begin{array}{ll}1,& x>0,\\ 0,& x=0,\\ -1,& x<0.\end{array}\right.\end{array}$ The mathematical model of moisture transfer is now defined by a system of one partial differential equation coupled with an ordinary differential equation: (22) $c^{\star }\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}=\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}^{\star }}\left(d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}-\text{P'e}u\right),$(22) (23) $\frac{\mathrm{\partial }{c}^{\star }}{\mathrm{\partial }{t}^{\star }}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}\right)\cdot \left(c^{\star }-c_{\mathrm{ads}}^{\star }\right)\cdot \left(c^{\star }-c_{\mathrm{des}}^{\star }\right),$(23) associated with initial and boundary conditions.

The issue is now to estimate the five coefficients of the adsorption and desorption curves $(c_{a,1}^{\star },c_{a,2}^{\star },c_{d,0}^{\star },c_{d,1}^{\star },c_{d,2}^{\star })$. It should be noted that the demonstration of structural identifiability of the five parameters is provided in Appendix 2. When searching the OED, one need to compute the partial derivative of u according to each coefficient $c_{a,d,m}^{\star }$. Since the model is now composed of two differential equations, for each coefficients, two sensitivity functions have to be computed. For instance, for the coefficient $c_{d,0}^{\star }$, it is required to compute: ${\mathrm{\Theta }}_u=\frac{\mathrm{\partial }u}{\mathrm{\partial }{c}_{d,0}^{\star }}\mathrm{and}{\mathrm{\Theta }}_c=\frac{\mathrm{\partial }c}{\mathrm{\partial }{c}_{d,0}^{\star }}.$

From a mathematical point of view, the computation of ${\mathrm{\Theta }}_c$ is not possible since the function $\mathrm{sign}\left(x\right)$ is not differentiable in the classical sense for x=0. It is differentiable in the sense of distributions. Therefore, a regularized version of the hysteresis model (18(21) $\frac{\mathrm{\partial }{c}^{\star }}{\mathrm{\partial }{t}^{\star }}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}\right)\left(c^{\star }-c_{\mathrm{ads}}^{\star }\right)\left(c^{\star }-c_{\mathrm{des}}^{\star }\right).$(21) ): $\frac{\mathrm{\partial }{c}^{\star }}{\mathrm{\partial }{t}^{\star }}=\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}\right)\cdot \left(c^{\star }-c_{\mathrm{ads}}^{\star }\right)\cdot \left(c^{\star }-c_{\mathrm{des}}^{\star }\right),$ where the function $\mathrm{Rsign}:\mathbb{R}\to \mathbb{R}$ is defined as: $\mathrm{Rsign}\left(x\right)=\tanh \left(\delta x\right),$ with δ a sufficiently large real parameter. In this study, the following value is taken $\delta ={10}^8$. A numerical validation of the regularized hysteresis model is provided in Appendix 3.

6.2. Searching the OED

The OED is now searched among the experimental designs described in Section 3 for both single and multiple steps. The issue is to estimate one or several parameters among the five coefficients $(c_{a,1}^{\star },c_{a,2}^{\star },c_{d,0}^{\star },c_{d,1}^{\star },c_{d,2}^{\star })$. For this, the sensitivity functions ${\mathrm{\Theta }}_{u,a,d,m}=\mathrm{\partial }u/\mathrm{\partial }{c}_{a,d,m}^{\star }$ and ${\mathrm{\Theta }}_{c,a,d,m}=\mathrm{\partial }c/\mathrm{\partial }{c}_{a,d,m}^{\star },m\in \{0,1,2\}$ are computed using central finite-difference and a Runge–Kutta  approach provided by ode45 solver in Matlab^TM [53]. The equations to compute the sensitivity functions are given in Appendix 1.

First, the OED is searched using the adsorption coefficients estimated in Section 5 and the desorption coefficients obtained from literature [33]: $c_{a,1}^{\star }=-0.99,c_{a,2}^{\star }=1.003,c_{d,0}^{\star }=0.75,c_{d,1}^{\star }=-1.24,c_{d,2}^{\star }=\mathrm{0.89.}$

The Fourier number corresponds to the estimated one in Section 5 and equals $\text{Fo}=6.1\times {10}^{-3}$. The investigations are performed for both the single and multiple steps of relative humidity, identified in Table 2. For the sake of compactness, the results are illustrated only for the multiple-step experiments. By improving the model with hysteresis effects, one has to compute the sensitivity coefficients for both equations (19a(22) $c^{\star }\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}=\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}^{\star }}\left(d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}-\text{P'e}u\right),$(22) ) and (19b(23) $\frac{\mathrm{\partial }{c}^{\star }}{\mathrm{\partial }{t}^{\star }}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}\right)\cdot \left(c^{\star }-c_{\mathrm{ads}}^{\star }\right)\cdot \left(c^{\star }-c_{\mathrm{des}}^{\star }\right),$(23) ). The sensitivity coefficients ${\mathrm{\Theta }}_c$ and ${\mathrm{\Theta }}_u$ are shown in Figures 15(a,c) for the adsorption coefficients $c_{a,1}^{\star }$ and $c_{a,2}^{\star }$. The variation of the sensitivity coefficients follows the changes in the boundary conditions $u_{\mathrm{\infty }}$. The design 12 corresponds to the OED for each parameter to be estimated, as shown in Figure 14. On the contrary to the results obtained in Section 4, the design 20 does not allow to estimate the unknown parameters with accuracy. The influence of including the hysteresis effects in the mathematical model can be remarked in the determination of the OED. Indeed, by comparing Figures 15 and 16, it can be noticed that the sensitivity coefficients ${\mathrm{\Theta }}_c$ and ${\mathrm{\Theta }}_u$ have larger magnitudes for the OED than for the design 20. For multiple-step experiments, a strong correlation is observed between the five coefficients $(c_{a,1}^{\star },c_{a,2}^{\star },c_{d,0}^{\star },c_{d,1}^{\star },c_{d,2}^{\star })$: $\begin{array}{r}\begin{array}{cccccc}& c_{a,1}^{\star }& c_{a,2}^{\star }& c_{d,0}^{\star }& c_{d,1}^{\star }& c_{d,2}^{\star }\\ c_{a,1}^{\star }& 1& 0.99669& 0.99974& 0.99473& 0.97735\\ c_{a,2}^{\star }& & 1& 0.99663& 0.999& 0.99\\ c_{d,0}^{\star }& & & 1& 0.99567& 0.97941\\ c_{d,1}^{\star }& & & & 1& 0.99388\\ c_{d,2}^{\star }& & & & & 1\end{array}\end{array}$ indicating that it is not possible to estimate more than two parameters of the model for such experiments. For single case experiments, the correlation is lower: $\begin{array}{r}\begin{array}{cccccc}& c_{a,1}^{\star }& c_{a,2}^{\star }& c_{d,0}^{\star }& c_{d,1}^{\star }& c_{d,2}^{\star }\\ c_{a,1}^{\star }& 1& 0.98981& 0.97193& 0.89302& 0.78741\\ {}_{a,2}^{\star }& & 1& 0.98349& 0.93466& 0.85096\\ c_{d,0}^{\star }& & & 1& 0.9713& 0.90551\\ c_{d,1}^{\star }& & & & 1& 0.98031\\ c_{d,2}^{\star }& & & & & 1\end{array}\end{array}$

Thus, for this case, it is possible to estimate a couple of two parameters among $(c_{a,1}^{\star },c_{d,2}^{\star })$ or $(c_{a,2}^{\star },c_{d,2}^{\star })$. The OED results are synthesized in Table 4. For all experiments, the sensor should be located at the bottom of the material, near the impermeable face. The designs 2 ($u_{\mathrm{i}}=0.2$ to $u_{\mathrm{c}}=1.5$) and 12 ($u_{\mathrm{i}}=0.2$, $u_{\mathrm{c}}^1=0.66$, $u_{\mathrm{c}}^2=1.5$, $u_{\mathrm{c}}^3=0.66$ and a duration of each step $\tau =8$) provide the highest accuracy to estimate the parameters. Due to high correlation between the coefficients, it is important to note that it is not possible to estimate all five coefficients considering the possible designs described in Table 2. If one aims at estimating the five coefficients, additional designs have to be planned with the experimental facility, including for instance other levels of relative humidity.

Figure 14. Variation of the criterion Ψ for the 16 possible designs for the adsorption coefficients (a) and the desorption coefficients (b).

Figure 15. Sensitivity coefficients ${\mathrm{\Theta }}_u$ (a,b) and ${\mathrm{\Theta }}_c$ (c,d) for the OED (design 12, $X=X^{\circ }=1$).

Figure 16. Sensitivity coefficients ${\mathrm{\Theta }}_u$ (a,b) and ${\mathrm{\Theta }}_c$ (c,d) for the design 20, ($X=X^{\circ }=1$).

Table 4. Synthesis of OED for the estimation of one or several parameters when the model account for hysteresis.

		Optimal Experimental Design
	Parameter(s) to be estimated	Single-step	Multiple-step	Sensor position
One parameter	$c_{a,1}^{\star }$, $c_{a,2}^{\star }$, $c_{d,0}^{\star }$, $c_{d,1}^{\star }$ or $c_{d,2}^{\star }$	design 2	design 12	$X^{\circ }\in [0.85,1]$
Multiple parameter	$(c_{a,1}^{\star },c_{d,2}^{\star })$ or $(c_{a,2}^{\star },c_{d,2}^{\star })$	design 2	–	$X^{\circ }\in [0.85,1]$

6.3. Comparing the numerical results with the experimental observations

Previous section aimed at illustrating the possibility of searching the OED with an improved model that includes the hysteresis effects. Since it is rather difficult to estimate the five coefficients due to strong correlations of the sensitivity functions, the purpose is now to show the influence of taking into account the hysteresis effects in the model. It should be noted that the regularized model is not needed so that the normal hysteresis model given by Equation (19b(23) $\frac{\mathrm{\partial }{c}^{\star }}{\mathrm{\partial }{t}^{\star }}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}\right)\cdot \left(c^{\star }-c_{\mathrm{ads}}^{\star }\right)\cdot \left(c^{\star }-c_{\mathrm{des}}^{\star }\right),$(23) ) is used for the present case study. To avoid stability restrictions, an implicit–explicit numerical scheme, detailed in Appendix 4, is employed to compute numerically the solution.

The numerical predictions are compared with the experimental observations for the multiple-step experiment (design 20). The adsorption coefficients, estimated in Section 5, are used together with the desorption coefficients obtained from literature [33]. The Fourier number equals $\mathrm{Fo}=6.1\times {10}^{-3}$ and the other parameters have the same numerical value as the ones mentioned in Section 5. Figure 17(a) shows the comparison between the numerical predictions and the experimental data. The results from the physical model including the hysteresis effects are more satisfactory. Indeed, the residual is lower for this model than for the model without hysteresis, as shown in Figure 17(b). Particularly, the importance of the hysteresis can be noted for the desorption phase $t\in [200,300]$ and the second adsorption phase $t\in [400,600]$. The solution of Equation (19b(23) $\frac{\mathrm{\partial }{c}^{\star }}{\mathrm{\partial }{t}^{\star }}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}\right)\cdot \left(c^{\star }-c_{\mathrm{ads}}^{\star }\right)\cdot \left(c^{\star }-c_{\mathrm{des}}^{\star }\right),$(23) ) enables to compute the time evolution of the sorption coefficient c as illustrated in Figure 18(a). It is compared with the function $(c\circ u)\left(t\right)$ where u is computed using the model without hysteresis given by Equation (12(14) $\begin{array}{rlrl}{c}^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}& =\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}^{\star }}\left(d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}-\text{P'e}u\right),& t^{\star }& >0,x^{\star }\in \left[0,1\right],\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=\text{Bi}\cdot \left(u-u_{\mathrm{\infty }}\right),& t^{\star }& >0,x^{\star }=0,\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=0,& t^{\star }& >0,x^{\star }=1,\\ u& =u_{\mathrm{i}},& t^{\star }& =0,x^{\star }\in \left[0,1\right],\end{array}$(14) ). The time variation of the coefficients is similar for the first adsorption part corresponding to $t\in [0,200]$. Indeed, for this period, the computed sorption coefficient equals the adsorption curve as noticed in Figure 18(b). After this period, the coefficient computed with the hysteresis model decreases comparing to the one without hysteresis and oscillates between the adsorption and desorption curves as shown in Figure 18(b).

Figure 17. Comparison of the numerical results with the experimental data (a) and their residual (b).

Figure 18. Variation of the sorption coefficient c according to time (a) and according to the computed field u (b).

7. Final remarks

7.1. Conclusion

In the context of building physics, inverse problems are encountered to estimate moisture-dependent hygrothermal properties of porous materials, using measurements associated to heat and moisture transport. Two applications are distinguished. In the first case, concerning the diagnosis of existing building walls, there is a few a priori estimation of material properties. Moreover, measurements must be non-intrusive and non-destructive. In the second case, measurements are performed in the laboratory to calibrate the numerical model with the experimental data. This article is encompassed in these conditions, focused on the estimation of moisture sorption isotherm coefficients of a wood fibre material.

First, the OED methodology has been described and used for searching the OED, ensuring to provide the best accuracy of the identification method for the parameter estimation. The approach is based on the sensitivity functions of the unknown parameters, enabling to determine sensor location within the material and boundary conditions, according to an existing facility among 20 possible designs. It has been carried out considering a priori values of the unknown parameters. The facility allows to submit material to relative humidity steps on one surface, being all others moisture impermeable. Results have enhanced two designs: (i) single step of relative humidity from $10\mathrm{\%}$ to $75\mathrm{\%}$ and (ii) multiple steps of relative humidity $10-75-33-75\mathrm{\%}$, with a duration period of 8 days for each step. For each design, the sensor has to be placed as close as possible to the impermeable boundary.

Then, experimental data have been provided according to the OED results for the two selected designs. The parameter estimation has been conducted by minimizing a cost function between the experimental data and the numerical results. The estimation has been accomplished using an interior point algorithm. Nine tests have been performed for the definition of the cost function $\text{J}$. The ${\mathcal{L}}_2$ and ${\mathcal{L}}_{\mathrm{\infty }}$ have been evaluated. Two series of tests aimed at estimating the three parameters using both experiments at the same time or separately with an iterative algorithm. The third series intended to estimate all five parameters of the material properties. Results have shown that the ${\mathcal{L}}_2$ norm provided better results of the parameter estimation problem. Moreover, it was better to consider both experiments separately to estimate only three parameters of the problem. Within this approach, the algorithm requires only two iterations to compute the solution with less than 100 computations of the direct model. This approach has a really low computational cost compared to stochastic approaches, needing an order of ${10}^4$ computations for similar problems. Another advantage of this approach is to use the sensitivity functions, computed during the search of the OED, to provide an approximation of the probability distribution function of the estimated parameters at a lower computational cost.

As highlighted in Section 5, with the estimated parameters, a better agreement between the numerical model and the experimental data is observed. However, the importance of hysteresis effects were highlighted. Particularly, when cycles of desorption–adsorption processes take place, some discrepancies remain between experimental data and numerical predictions. Therefore, a new mathematical model has been proposed to take into account the hysteresis effects on the sorption coefficients. A second differential equation has been added and increased to five the number of coefficients to be estimated. Two coefficients correspond to the adsorption curve and three to the desorption one. A regularized version of the hysteresis model was proposed to have continuous differentiable functions and, therefore, to be able to compute the sensitivity coefficients. Then, the OED was explored by computing the sensitivity coefficients of the five parameters of a family of scanning curves of adsorption and desorption processes. This clearly highlighted the possibility of including the hysteresis effects in the OED approach. The results draw attention to two designs: (i) a single step of relative humidity from $10\mathrm{\%}$ to $75\mathrm{\%}$ and (ii) multiple steps of relative humidity $10-33-75-33\mathrm{\%}$, with a 8-day time period for each step. The sensors have to be placed near the impermeable boundary.

Finally, the numerical predictions, considering the hysteresis phenomenon, have been compared with the experimental observations of a multiple-step design. An efficient implicit–explicit numerical scheme was proposed to compute the solution of the hysteretic model. The parameters of this model correspond to the estimated ones in Section 5 for the adsorption curve and to the a priori ones provided in the literature. The comparison has shown that the discrepancies are reduced, fitting better experimental data. During the simulation, the computed sorption coefficients oscillated between the ones from the main adsorption and desorption curves.

7.2. Outlooks and open-problems

An interesting perspective of improvement concerns the assumptions related to the moisture sorption isotherm coefficients $c\left(u\right)$. A parametrization was previously defined $c\left(u\right)=1+c_{1}u+c_2{u}^2$ and the parameter estimation problem aimed at identifying coefficients $c_1$, $c_2$ (and $\text{Fo}$). An ambitious outlook could aim at estimating directly the function $c\left(u\right)$, with inspiration from the following studies [54,55].

Nomenclature

Table

Latin letters
a	moisture advection coefficient ($\mathrm{s}/\mathrm{m}$)
$a_m$	mass transfer coefficient (${\mathrm{m}}^2/\mathrm{s}$)
c	moisture storage capacity ($\mathrm{kg}/({\mathrm{m}}^3.\mathrm{Pa})$)
$c_b$	specific heat $\left(\frac{\mathrm{J}}{\left(\mathrm{kg}.\mathrm{K}\right)}\right)$
d	moisture permeability (s)
h	convective vapour transfer coefficient ($\mathrm{s}/\mathrm{m}$)
L	length (m)
$P_s$	saturation pressure (Pa)
$P_v$	vapour pressure (Pa)
$r_{12}$	latent heat of evaporation (J/kg)
$R_v$	water gas constant (J/(kg.K))
T	temperature (K)
t	time coordinate (s)
U	relative moisture concentration (%)
x	space coordinate (m)
$\mathsf{v}$	mass average velocity (m/s)
w	specific moisture content(${\mathrm{kg}/\mathrm{m}}^3$)

Table

Greek letters
λ	thermal conductivity (W/(m.K))
δ	thermal-gradient coefficient (${\mathrm{K}}^{-1}$)
φ	relative humidity (-)
ρ	specific mass (${\mathrm{kg}/\mathrm{m}}^3$)

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was partly funded by the French Environment and Energy Management Agency (Agence de lEnvironnement et de la Maîtrise de lEnergie, ADEME), the “Assemblée des Pays de Savoie” (APS) and the French National Research Agency (Agence Nationale de la Recherche, ANR) through its Sustainable Cities and Buildings programme (MOBAIR project ANR-12-VBDU-0009). The authors acknowledge the Junior Chair Research programme “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) and the support of Centre National de la Recherche Scientifique (CNRS/INSIS) (INSIS, Cellule Energie, Projets Exploratoires, Cellule énergie) under the programme “Projets Exploratoires – 2017”. The authors also acknowledge the Brazilian Agency CNPQ of the Ministry of Science, Technology and Innovation, for the financial support.

Appendices

A.1. Appendix 1. Equations for the computation of the sensitivity coefficients

This section provides the equations derived analytically from the mathematical model (Equation (12(14) $\begin{array}{rlrl}{c}^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}& =\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}^{\star }}\left(d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}-\text{P'e}u\right),& t^{\star }& >0,x^{\star }\in \left[0,1\right],\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=\text{Bi}\cdot \left(u-u_{\mathrm{\infty }}\right),& t^{\star }& >0,x^{\star }=0,\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=0,& t^{\star }& >0,x^{\star }=1,\\ u& =u_{\mathrm{i}},& t^{\star }& =0,x^{\star }\in \left[0,1\right],\end{array}$(14) )), to compute the sensitivity coefficients.

A.1.1 Model without hysteresis

The sensitivity coefficients are denoted as follows: ${\mathrm{\Theta }}_0\stackrel{\mathrm{def}}{:=}\frac{\mathrm{\partial }u}{\mathrm{\partial }\text{Fo}},{\mathrm{\Theta }}_1\stackrel{\mathrm{def}}{:=}\frac{\mathrm{\partial }u}{\mathrm{\partial }{c}_1},{\mathrm{\Theta }}_2\stackrel{\mathrm{def}}{:=}\frac{\mathrm{\partial }u}{\mathrm{\partial }{c}_2}.$

For the sake of clarity, the superscript ${}^{\star }$ is omitted. The sensitivity coefficients verified the following equations. For ${\mathrm{\Theta }}_0$: $\begin{array}{rl}c\left(u\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_0}{\mathrm{\partial }t}=& \text{Fo}d\left(u\right)\frac{\mathrm{\partial }{}^2{\mathrm{\Theta }}_0}{\mathrm{\partial }{x}^2}+\text{Fo}\left(2d_1\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_0}{\mathrm{\partial }x}+\text{Fo}{d}_1{\mathrm{\Theta }}_0\frac{\mathrm{\partial }{}^{2}u}{\mathrm{\partial }{x}^2}\\ & +d_1{\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }x}\right)}^2+d\left(u\right)\frac{\mathrm{\partial }{}^{2}u}{\mathrm{\partial }{x}^2}-\text{P'e}\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\left(c_1{\mathrm{\Theta }}_0+2c_{2}u{\mathrm{\Theta }}_0\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }t}.\end{array}$ For ${\mathrm{\Theta }}_1$: $\begin{array}{rl}c\left(u\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_1}{\mathrm{\partial }t}=& \text{Fo}d\left(u\right)\frac{\mathrm{\partial }{}^2{\mathrm{\Theta }}_1}{\mathrm{\partial }{x}^2}+2\text{Fo}{d}_1\frac{\mathrm{\partial }{\mathrm{\Theta }}_1}{\mathrm{\partial }x}\frac{\mathrm{\partial }u}{\mathrm{\partial }x}+\text{Fo}{d}_1{\mathrm{\Theta }}_1\frac{\mathrm{\partial }{}^{2}u}{\mathrm{\partial }{x}^2}\\ & -\text{Fo}\text{P'e}\frac{\mathrm{\partial }{\mathrm{\Theta }}_1}{\mathrm{\partial }x}-\left(u+c_1{\mathrm{\Theta }}_1+2c_{2}u{\mathrm{\Theta }}_1\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }t}.\end{array}$ and for ${\mathrm{\Theta }}_2$: $\begin{array}{rl}c\left(u\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_2}{\mathrm{\partial }t}=& \text{Fo}d\left(u\right)\frac{\mathrm{\partial }{}^2{\mathrm{\Theta }}_2}{\mathrm{\partial }{x}^2}+2\text{Fo}{d}_1\frac{\mathrm{\partial }{\mathrm{\Theta }}_2}{\mathrm{\partial }x}\frac{\mathrm{\partial }u}{\mathrm{\partial }x}+\text{Fo}{d}_1{\mathrm{\Theta }}_2\frac{\mathrm{\partial }{}^{2}u}{\mathrm{\partial }{x}^2}\\ & -\text{Fo}\text{P'e}\frac{\mathrm{\partial }{\mathrm{\Theta }}_2}{\mathrm{\partial }x}-\left(c_1{\mathrm{\Theta }}_2+2c_{2}u{\mathrm{\Theta }}_2+u^2\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }t},\end{array}$

A.1.2 Model with hysteresis

The model with hysteresis includes two differential equations, recalled here: $\begin{array}{rl}c\frac{\mathrm{\partial }u}{\mathrm{\partial }t}& =\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(d\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}u\right),\\ \frac{\mathrm{\partial }c}{\mathrm{\partial }t}& =\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c-c_{\mathrm{ads}}\left(u\right)\right)\cdot \left(c-c_{\mathrm{des}}\left(u\right)\right).\end{array}$ Therefore, two sensitivity coefficients have to be computed. The differential equations for ${\mathrm{\Theta }}_{u,a,1}=\mathrm{\partial }u/\mathrm{\partial }{c}_{a,1}$ and ${\mathrm{\Theta }}_{c,a,1}=\mathrm{\partial }c/\mathrm{\partial }{c}_{a,1}$ are: $\begin{array}{rl}\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,a,1}}{\mathrm{\partial }t}=& \text{Fo}d\left(u\right)\frac{\mathrm{\partial }{}^2{\mathrm{\Theta }}_{u,a,1}}{\mathrm{\partial }{x}^2}+\text{Fo}\left(2d_1\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,a,1}}{\mathrm{\partial }x}+\text{Fo}{d}_1\frac{\mathrm{\partial }{}^{2}u}{\mathrm{\partial }{x}^2}{\mathrm{\Theta }}_{u,a,1}\\ & -{\mathrm{\Theta }}_{c,a,1}\frac{\mathrm{\partial }u}{\mathrm{\partial }t},\\ \frac{\mathrm{\partial }{\mathrm{\Theta }}_{c,a,1}}{\mathrm{\partial }t}=& \beta {\mathrm{Rsign}}^{\mathrm{\prime }}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,a,1}}{\mathrm{\partial }t}\left(c-c_{\mathrm{ads}}\right)\cdot \left(c-c_{\mathrm{des}}\right)\\ & +\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left({\mathrm{\Theta }}_{c,a,1}-\left(c_{1,a}{\mathrm{\Theta }}_{u,a,1}+u+2c_{2,a}u{\mathrm{\Theta }}_{u,a,1}\right)\right)\cdot \left(c-c_{\mathrm{des}}\right)\\ & +\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c-c_{\mathrm{ads}}\right)\cdot \left({\mathrm{\Theta }}_{c,a,1}-\left(c_{1,d}{\mathrm{\Theta }}_{u,a,1}+2c_{2,d}u{\mathrm{\Theta }}_{u,a,1}\right)\right).\end{array}$ For ${\mathrm{\Theta }}_{u,a,2}=\mathrm{\partial }u/\mathrm{\partial }{c}_{a,2}$ and ${\mathrm{\Theta }}_{c,a,2}=\mathrm{\partial }c/\mathrm{\partial }{c}_{a,2}$: $\begin{array}{rl}\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,a,2}}{\mathrm{\partial }t}=& \text{Fo}d\left(u\right)\frac{\mathrm{\partial }{}^2{\mathrm{\Theta }}_{u,a,2}}{\mathrm{\partial }{x}^2}+\text{Fo}\left(2d_1\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,a,2}}{\mathrm{\partial }x}+\text{Fo}{d}_1\frac{\mathrm{\partial }{}^{2}u}{\mathrm{\partial }{x}^2}{\mathrm{\Theta }}_{u,a,2}\\ & -{\mathrm{\Theta }}_{c,a,2}\frac{\mathrm{\partial }u}{\mathrm{\partial }t},\\ \frac{\mathrm{\partial }{\mathrm{\Theta }}_{c,a,2}}{\mathrm{\partial }t}=& \beta {\mathrm{Rsign}}^{\mathrm{\prime }}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,a,2}}{\mathrm{\partial }t}\left(c-c_{\mathrm{ads}}\right)\cdot \left(c-c_{\mathrm{des}}\right)\\ & +\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left({\mathrm{\Theta }}_{c,a,2}-\left(c_{1,a}{\mathrm{\Theta }}_{u,a,2}+u^2+2c_{2,a}u{\mathrm{\Theta }}_{u,a,2}\right)\right)\cdot \left(c-c_{\mathrm{des}}\right)\\ & +\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c-c_{\mathrm{ads}}\right)\cdot \left({\mathrm{\Theta }}_{c,a,2}-\left(c_{1,d}{\mathrm{\Theta }}_{u,a,2}+2c_{2,d}u{\mathrm{\Theta }}_{u,a,2}\right)\right).\end{array}$ For ${\mathrm{\Theta }}_{u,d,0}=\mathrm{\partial }u/\mathrm{\partial }{c}_{d,0}$ and ${\mathrm{\Theta }}_{c,d,0}=\mathrm{\partial }c/\mathrm{\partial }{c}_{d,0}$: $\begin{array}{rl}\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,d,0}}{\mathrm{\partial }t}=& \text{Fo}d\left(u\right)\frac{\mathrm{\partial }{}^2{\mathrm{\Theta }}_{u,d,0}}{\mathrm{\partial }{x}^2}+\text{Fo}\left(2d_1\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,d,0}}{\mathrm{\partial }x}+\text{Fo}{d}_1\frac{\mathrm{\partial }{}^{2}u}{\mathrm{\partial }{x}^2}{\mathrm{\Theta }}_{u,d,0}\\ & -{\mathrm{\Theta }}_{c,d,0}\frac{\mathrm{\partial }u}{\mathrm{\partial }t},\\ \frac{\mathrm{\partial }{\mathrm{\Theta }}_{c,d,0}}{\mathrm{\partial }t}=& \beta {\mathrm{Rsign}}^{\mathrm{\prime }}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,d,0}}{\mathrm{\partial }t}\left(c-c_{\mathrm{ads}}\right)\cdot \left(c-c_{\mathrm{des}}\right)\\ & +\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left({\mathrm{\Theta }}_{c,d,0}-\left(c_{1,a}{\mathrm{\Theta }}_{u,d,0}+2c_{2,a}u{\mathrm{\Theta }}_{u,d,0}\right)\right)\cdot \left(c-c_{\mathrm{des}}\right)\\ & +\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c-c_{\mathrm{ads}}\right)\cdot \left({\mathrm{\Theta }}_{c,d,0}-\left(1+c_{1,d}{\mathrm{\Theta }}_{u,d,0}+2c_{2,d}u{\mathrm{\Theta }}_{u,d,0}\right)\right).\end{array}$ For ${\mathrm{\Theta }}_{u,d,1}=\mathrm{\partial }u/\mathrm{\partial }{c}_{d,1}$ and ${\mathrm{\Theta }}_{c,d,1}=\mathrm{\partial }c/\mathrm{\partial }{c}_{d,1}$: $\begin{array}{rl}\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,d,1}}{\mathrm{\partial }t}=& \text{Fo}d\left(u\right)\frac{\mathrm{\partial }{}^2{\mathrm{\Theta }}_{u,d,1}}{\mathrm{\partial }{x}^2}+\text{Fo}\left(2d_1\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,d,1}}{\mathrm{\partial }x}+\text{Fo}{d}_1\frac{\mathrm{\partial }{}^{2}u}{\mathrm{\partial }{x}^2}{\mathrm{\Theta }}_{u,d,1}\\ & -{\mathrm{\Theta }}_{c,d,1}\frac{\mathrm{\partial }u}{\mathrm{\partial }t},\\ \frac{\mathrm{\partial }{\mathrm{\Theta }}_{c,d,1}}{\mathrm{\partial }t}=& \beta {\mathrm{Rsign}}^{\mathrm{\prime }}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,d,1}}{\mathrm{\partial }t}\left(c-c_{\mathrm{ads}}\right)\cdot \left(c-c_{\mathrm{des}}\right)\\ & +\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left({\mathrm{\Theta }}_{c,d,1}-\left(c_{1,a}{\mathrm{\Theta }}_{u,d,1}+2c_{2,a}u{\mathrm{\Theta }}_{u,d,1}\right)\right)\cdot \left(c-c_{\mathrm{des}}\right)\\ & +\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c-c_{\mathrm{ads}}\right)\cdot \left({\mathrm{\Theta }}_{c,d,1}-\left(c_{1,d}{\mathrm{\Theta }}_{u,d,1}+u+2c_{2,d}u{\mathrm{\Theta }}_{u,d,1}\right)\right).\end{array}$ For ${\mathrm{\Theta }}_{u,d,2}=\mathrm{\partial }u/\mathrm{\partial }{c}_{d,2}$ and ${\mathrm{\Theta }}_{c,d,2}=\mathrm{\partial }c/\mathrm{\partial }{c}_{d,2}$: $\begin{array}{rl}\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,d,2}}{\mathrm{\partial }t}=& \text{Fo}d\left(u\right)\frac{\mathrm{\partial }{}^2{\mathrm{\Theta }}_{u,d,2}}{\mathrm{\partial }{x}^2}+\text{Fo}\left(2d_1\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,d,2}}{\mathrm{\partial }x}+\text{Fo}{d}_1\frac{\mathrm{\partial }{}^{2}u}{\mathrm{\partial }{x}^2}{\mathrm{\Theta }}_{u,d,2}\\ & -{\mathrm{\Theta }}_{c,d,2}\frac{\mathrm{\partial }u}{\mathrm{\partial }t},\\ \frac{\mathrm{\partial }{\mathrm{\Theta }}_{c,d,2}}{\mathrm{\partial }t}=& \beta {\mathrm{Rsign}}^{\mathrm{\prime }}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\frac{\mathrm{\partial }{\mathrm{\Theta }}_{u,d,2}}{\mathrm{\partial }t}\left(c-c_{\mathrm{ads}}\right)\left(c-c_{\mathrm{des}}\right)\\ & +\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left({\mathrm{\Theta }}_{c,d,2}-\left(c_{1,a}{\mathrm{\Theta }}_{u,d,2}+2c_{2,a}u{\mathrm{\Theta }}_{u,d,2}\right)\right)\cdot \left(c-c_{\mathrm{des}}\right)\\ & +\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c-c_{\mathrm{ads}}\right)\cdot \left({\mathrm{\Theta }}_{c,d,2}-\left(c_{1,d}{\mathrm{\Theta }}_{u,d,2}+u^2+2c_{2,d}u{\mathrm{\Theta }}_{u,d,2}\right)\right).\end{array}$ The function ${\mathrm{Rsign}}^{\mathrm{\prime }}:\mathbb{R}\to \{0,1\}$ is the regularized Dirac function: $\begin{array}{r}{\mathrm{Rsign}}^{\mathrm{\prime }}\left(x\right)=\left\{\begin{array}{ll}1,& x=0,\\ 0,& x\ne 0.\end{array}\right.\end{array}$

Figure A1. Time variation of the field u (a) and of the sorption capacity c (b) for the model with hysteresis and the regularised one.

Figure A2. Variation of the sorption capacity c with the field u (a) according to the case study. Time variation of the ${\mathcal{L}}_2$ error (b) for the fields u and c between the model with hysteresis and the regularised one.

A.2 Appendix 2. Proving the structural identifiability of the parameters

This section aims at justifying the identifiability of the unknown parameters for both model with and without hysteresis. We recall that a parameter P is Structurally Globally Identifiable (SGI) if the following condition is satisfied [56]: $y\left(P\right)\equiv y\left(P^{\mathrm{\prime }}\right)⟹P\equiv P^{\mathrm{\prime }},$

where y is the response of the model depending on parameter P.

A.2.2 Model without hysteresis

First, the SGI is demonstrated for the estimation of parameters $(\text{Fo},c_1,c_2)$ in the model (12(14) $\begin{array}{rlrl}{c}^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}& =\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}^{\star }}\left(d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}-\text{P'e}u\right),& t^{\star }& >0,x^{\star }\in \left[0,1\right],\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=\text{Bi}\cdot \left(u-u_{\mathrm{\infty }}\right),& t^{\star }& >0,x^{\star }=0,\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=0,& t^{\star }& >0,x^{\star }=1,\\ u& =u_{\mathrm{i}},& t^{\star }& =0,x^{\star }\in \left[0,1\right],\end{array}$(14) ), recalled here without the superscript ${}^{\star }$: (24) $\left(1+c_{1}u+c_2{u}^2\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }t}-\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\left(1+d_{1}u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}u\right)=0,$(24) It is assumed that u is observable. So as to prove identifiability, it is assumed that another set of parameters, denoted with a superscript ${}^{\mathrm{\prime }}$, holds: (25) $\left(1+c_1^{\mathrm{\prime }}{u}^{\mathrm{\prime }}+c_2^{\mathrm{\prime }}(u^{\mathrm{\prime }}{)}^2\right)\frac{\mathrm{\partial }{u}^{\mathrm{\prime }}}{\mathrm{\partial }t}-{\text{Fo}}^{\mathrm{\prime }}\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\left(1+d_1{u}^{\mathrm{\prime }}\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}{u}^{\mathrm{\prime }}\right)=0,$(25) If $u\equiv u^{\mathrm{\prime }}$ then $\mathrm{\partial }u/\mathrm{\partial }t\equiv \mathrm{\partial }{u}^{\mathrm{\prime }}/\mathrm{\partial }t$ and $\mathrm{\partial }u/\mathrm{\partial }x\equiv \mathrm{\partial }{u}^{\mathrm{\prime }}/\mathrm{\partial }x$. Thus, from Equations (A1(24) $\left(1+c_{1}u+c_2{u}^2\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }t}-\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\left(1+d_{1}u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}u\right)=0,$(24) ) and (A2(25) $\left(1+c_1^{\mathrm{\prime }}{u}^{\mathrm{\prime }}+c_2^{\mathrm{\prime }}(u^{\mathrm{\prime }}{)}^2\right)\frac{\mathrm{\partial }{u}^{\mathrm{\prime }}}{\mathrm{\partial }t}-{\text{Fo}}^{\mathrm{\prime }}\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\left(1+d_1{u}^{\mathrm{\prime }}\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}{u}^{\mathrm{\prime }}\right)=0,$(25) ) we have: $\left(1+c_{1}u+c_2{u}^2\right)\equiv \left(1+c_1^{\mathrm{\prime }}{u}^{\mathrm{\prime }}+c_2^{\mathrm{\prime }}(u^{\mathrm{\prime }}{)}^2\right),$

and $\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\left(1+d_{1}u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}u\right)\equiv {\text{Fo}}^{\mathrm{\prime }}\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\left(1+d_1{u}^{\mathrm{\prime }}\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }x}-\text{P'e}{u}^{\mathrm{\prime }}\right).$ Since $u\equiv u^{\mathrm{\prime }}$ and $u^2\equiv (u^{\mathrm{\prime }}{)}^2$, it follows that $c_1\equiv c_1^{\mathrm{\prime }},c_2\equiv c_2^{\mathrm{\prime }},{\text{Fo}}^{\mathrm{\prime }}\equiv {\text{Fo}}^{\mathrm{\prime }}.$ Therefore, parameters $(\text{Fo},c_1,c_2)$ are SGI for the model without hysteresis.

A.2.3 Model with hysteresis

Now, the SGI for the five parameters $(c_{a,1}^{\star },c_{a,2}^{\star },c_{d,0}^{\star },c_{d,1}^{\star },c_{d,2}^{\star })$ is demonstrated for the model with hysteresis. For this, Equation (18(21) $\frac{\mathrm{\partial }{c}^{\star }}{\mathrm{\partial }{t}^{\star }}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}\right)\left(c^{\star }-c_{\mathrm{ads}}^{\star }\right)\left(c^{\star }-c_{\mathrm{des}}^{\star }\right).$(21) ) is recalled omitting the superscript ${}^{\star }$: $\frac{\mathrm{\partial }c}{\mathrm{\partial }t}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\left(c-\left(1+c_{a,1}u+c_{a,2}{u}^2\right)\right)\left(c-\left(c_{d,0}+c_{d,1}u+c_{d,2}{u}^2\right)\right).$

Similarly, to prove the identifiability, another set of parameters is assumed: $\begin{array}{rl}& \frac{\mathrm{\partial }{c}^{\mathrm{\prime }}}{\mathrm{\partial }t}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }{u}^{\mathrm{\prime }}}{\mathrm{\partial }t}\right)\left(c^{\mathrm{\prime }}-\left(1+c_{a,1}^{\mathrm{\prime }}{u}^{\mathrm{\prime }}+c_{a,2}^{\mathrm{\prime }}(u^{\mathrm{\prime }}{)}^2\right)\right)\\ & \left(c^{\mathrm{\prime }}-\left(c_{d,0}^{\mathrm{\prime }}+c_{d,1}^{\mathrm{\prime }}{u}^{\mathrm{\prime }}+c_{d,2}^{\mathrm{\prime }}(u^{\mathrm{\prime }}{)}^2\right)\right).\end{array}$ It is assumed that $c\equiv c^{\mathrm{\prime }}$, $u\equiv u^{\mathrm{\prime }}$ and then $\mathrm{\partial }c/\mathrm{\partial }t\equiv \mathrm{\partial }{c}^{\mathrm{\prime }}/\mathrm{\partial }t$. Using the symbolic code Maple^TM, by expansion it can be demonstrated that: $c_{a,1}\equiv c_{a,1}^{\mathrm{\prime }},c_{a,2}\equiv c_{a,2}^{\mathrm{\prime }},c_{d,0}\equiv c_{d,0}^{\mathrm{\prime }},c_{d,1}\equiv c_{d,1}^{\mathrm{\prime }},c_{d,,2}\equiv c_{d,2}^{\mathrm{\prime }},$ and that the parameters are SGI.

A.3 Appendix 3. Numerical validation of the regularized hysteresis model

A computation using the regularized hysteresis model is carried out: $\begin{array}{rl}{c}^{\star }\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}& =\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}^{\star }}\left(d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}-\text{P'e}u\right),\\ \frac{\mathrm{\partial }{c}^{\star }}{\mathrm{\partial }t}& =\beta \mathrm{Rsign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c^{\star }-c_{\mathrm{ads}}^{\star }\left(u\right)\right)\cdot \left(c^{\star }-c_{\mathrm{des}}^{\star }\left(u\right)\right),\end{array}$ where the numerical values of the coefficients are: $\begin{array}{rl}& \text{P'e}=1\cdot {10}^{-2},\text{Fo}=2\cdot {10}^{-2},d^{\star }\left(u\right)=0.86+0.25u,\beta ={10}^{-3},\\ & c_{\mathrm{des}}^{\star }\left(u\right)=2.54-4.17u+3.02u^2,c_{\mathrm{ads}}^{\star }\left(u\right)=3.36-6.11u+3.37u^2.\end{array}$ The initial and boundary conditions are defined in Equation (12(14) $\begin{array}{rlrl}{c}^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{t}^{\star }}& =\text{Fo}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}^{\star }}\left(d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}-\text{P'e}u\right),& t^{\star }& >0,x^{\star }\in \left[0,1\right],\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=\text{Bi}\cdot \left(u-u_{\mathrm{\infty }}\right),& t^{\star }& >0,x^{\star }=0,\\ d^{\star }\left(u\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }{x}^{\star }}& -a^{\star }u=0,& t^{\star }& >0,x^{\star }=1,\\ u& =u_{\mathrm{i}},& t^{\star }& =0,x^{\star }\in \left[0,1\right],\end{array}$(14) ) with the following numerical values: $\begin{array}{r}\text{Bi}=13.7,u_{\mathrm{i}}=0.2,u_{\mathrm{\infty }}=\left\{\begin{array}{ll}1.5,& t>0,t<100,\\ 0.5,& t>100,t<200,\\ 1.0,& t>200,t<300,\\ 0.2,& t>300,t<400,\\ 1.8,& t>400,t<500.\end{array}\right.\end{array}$ The case study is similar to the experimental designs investigated in this work. The simulation time horizon is t=700. The time variations of the field u are given in Figure A1(a). The results of the regularized and non-regularized model are compared to the model with no-hysteresis, considering only the adsorption curve. The importance of the hysteresis in the time variations of u can be noted. Figure A1(b) shows the time variation of c by computation of its differential equation. The variations of the sorption coefficient c with the computed values of u are shown in Figure A2(a). In both Figures A1(b) and A2(a), a perfect agreement is observed between the regularized and non-regularized models. Moreover, as noticed in Figure A2(b), the ${\mathcal{L}}_2$ error of the fields u and c, computed between the regularized and non-regularized models, scales with ${10}^{-5}$ and ${10}^{-3}$, respectively. The agreement between the results of the two models is very satisfactory, validating the implementation of the regularized hysteretic model.

A.4 Appendix 4. Implicit–explicit numerical scheme for the hysteresis model

To relax stability restriction, an implicit–explicit numerical scheme is used to compute the solution of the hysteresis model: (26) $\frac{\mathrm{\partial }c}{\mathrm{\partial }t}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c-c_{\mathrm{ads}}\left(u\right)\right)\cdot \left(c-c_{\mathrm{des}}\left(u\right)\right).$(26) For the sake of clarity, the ⋆ superscript have been omitted. A uniform discretisation is considered for time intervals. The discretisation parameter is denoted using $\mathrm{\Delta }t$ for the time. The discrete values of function $c\left(t\right)$ are denoted by $c^n\stackrel{\mathrm{def}}{:=}c\left(t^n\right)$ with $n\in \{1,\dots ,N_t\}$.

When $\mathrm{\partial }u/\mathrm{\partial }t<0$, Equation (A3(26) $\frac{\mathrm{\partial }c}{\mathrm{\partial }t}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c-c_{\mathrm{ads}}\left(u\right)\right)\cdot \left(c-c_{\mathrm{des}}\left(u\right)\right).$(26) ) is discretized according to: $\frac{1}{\mathrm{\Delta }t}\left(c^{n+1}-c^n\right)=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c^{n+1}-c_{\mathrm{ads}}\right)\cdot \left(c^n-c_{\mathrm{des}}\right),$ which gives the explicit expression of $c^{n+1}$: $c^{n+1}=\frac{c^n-\mathrm{\Delta }t\beta \mathrm{sign}\left({\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }t}}\right)\left(c^n-c_{\mathrm{des}}\right)c_{\mathrm{ads}}}{1-\mathrm{\Delta }t\beta \mathrm{sign}\left({\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }t}}\right)\left(c^n-c_{\mathrm{des}}\right)}.$

Similarly, when $\mathrm{\partial }u/\mathrm{\partial }t>0$, Equation (A3(26) $\frac{\mathrm{\partial }c}{\mathrm{\partial }t}=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\cdot \left(c-c_{\mathrm{ads}}\left(u\right)\right)\cdot \left(c-c_{\mathrm{des}}\left(u\right)\right).$(26) ) is discretized according to: $\frac{1}{\mathrm{\Delta }t}\left(c^{n+1}-c^n\right)=\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\left(c^n-c_{\mathrm{ads}}\right)\left(c^{n+1}-c_{\mathrm{des}}\right),$ to obtain the explicit expression of $c^{n+1}$: $c^{n+1}=\frac{c^n-\mathrm{\Delta }t\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\left(c^n-c_{\mathrm{ads}}\right)c_{\mathrm{des}}}{1-\mathrm{\Delta }t\beta \mathrm{sign}\left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\right)\left(c^n-c_{\mathrm{ads}}\right)}.$ This numerical scheme provides robust and stable results as already shown in Figures 17 and 18.