Marouane Talbi

Marc Prat
email: mprat@imft.fr

Percolating and non-percolating liquid phase continuum model of drying in capillary porous media with application to solute transport in the very low Péclet number limit

A three equation continuum model of drying is presented. The model explicitly considers the liquid phase as formed by a percolating liquid phase and a non-percolating liquid phase. The model is tested against pore network simulations. A quite good agreement is obtained between the predictions of the continuum model and data obtained by volume averaging the pore network simulation results. Then, the model is extended to the case where a solute is present in the liquid phase. This leads to the consideration of a five equation continuum model as opposed to the classically considered two equation model. The model is tested when diffusion is the solute dominant transport mechanism. In agreement with the pore network simulations, the five equation continuum model predicts that the solute concentration in the percolating liquid phase is greater than in the non-percolating liquid phase in the considered situation. The work illustrates the key role of the liquid fragmentation process occurring during drying on the solute dynamics. Counter-intuitively, although diffusion is dominant, it is shown that the solution concentration varies over the liquid phase as the result of the liquid phase fragmentation process.

I. INTRODUCTION

Drying in porous media is a topic a great interest with applications in chemical engineering, [START_REF] Mujumdar | Handbook of Industrial Drying[END_REF], process engineering, [START_REF] Ito | Effect of through-plane distribution of polytetrafluoroethylene in carbon paper on in-plane gas permeability[END_REF], food processing [START_REF] Kumar | Intermittent drying of food products: A critical review[END_REF], civil engineering [START_REF] Villani | The Influence of Pore Solutions Properties on Drying in Cementitious Materials[END_REF], soil sciences [START_REF] Or | Advances in soil evaporation physics-A review[END_REF], to name only a few. Its modeling and numerical simulation has a long history as reviewed in [START_REF] Whitaker | Simultaneous heat, mass and momentum transfer in porous media. A theory of drying[END_REF] or [START_REF] Vu | Mass and heat transport models for analysis of the drying process in porous media: a review and numerical implementation[END_REF]. As presented in [START_REF] Whitaker | Simultaneous heat, mass and momentum transfer in porous media. A theory of drying[END_REF] and [START_REF] Vu | Mass and heat transport models for analysis of the drying process in porous media: a review and numerical implementation[END_REF], the drying process in capillary porous media, defined as porous media with pores greater than about 1µm in equivalent diameter, is commonly modelled within the framework of the continuum approach to porous media using a strongly non-linear diffusion equation governing the evolution of the saturation in the medium [START_REF] Pel | Analytic solution for the non-linear drying problem[END_REF]. However, this approach does not explicitly take into account an important feature of the drying process, namely the fact, sketched in Fig. 1, that the gradual replacement of the liquid by the gas phase in the pores resulting from the evaporation process leads to the fragmentation of the liquid phase into liquid clusters [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF]. FIG. 1. Schematic of liquid phase fragmentation process during drying with main cluster (percolating liquid phase) and isolated clusters (non-percolating liquid phase).

The fragmentation mechanism is inherent to the invasion percolation process [START_REF] Wilkinson | Invasion percolation: a new form of percolation theory[END_REF] at the core of the drying process, e.g. [START_REF] Prat | Percolation model of drying under isothermal conditions in porous media[END_REF]. One can refer to [START_REF] Prat | Isothermal drying of non-hygroscopic capillary-porous materials as an invasion percolation process[END_REF] for an illustration of the fragmentation process, i.e. the formation of a new cluster from an existing cluster as the result of a pore invasion. As discussed for instance in [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF], one can distinguish a first period in the drying process, referred to as stage 1, where the liquid phase is distributed between the largest liquid cluster spanning the porous medium and a number of smaller clusters. The largest cluster is referred to as the main cluster whereas the smaller clusters are referred to as the isolated clusters.

As shown in [START_REF] Ahmad | From micro-scale to macroscale modeling of solute transport in drying capillary porous media[END_REF], this structuration of the liquid phase has an impact on the evolution of the spatial and temporal distribution of a dissolved species during drying. In other words, the question arises as to whether the liquid phase fragmentation process can be taken into account in the continuum approach to porous media so as to develop better models of the transport of a solute during drying. This holds true as well for other drying situations where particles are present in the liquid phase, e.g. [START_REF] Ito | Effect of through-plane distribution of polytetrafluoroethylene in carbon paper on in-plane gas permeability[END_REF][START_REF] Keita | MRI evidence for a receding-front effect in drying porous media[END_REF]. To this end, we present in the current paper a three equation continuum model whose main variables are the vapor partial pressure, the main cluster saturation and the isolated clusters saturation. Regarding the liquid phase evolution, the three equation continuum model is based on the theory of biphasic flow in porous media presented in [START_REF] Hilfer | Capillary pressure, hysteresis and residual saturation in porous media[END_REF][START_REF] Hilfer | Macroscopic capillarity without a constitutive capillary pressure function[END_REF][START_REF] Hilfer | Macroscopic capillarity and hysteresis for flow in porous media[END_REF][START_REF] Doster | Numerical solutions of a generalized theory for macroscopic capillarity[END_REF] proposing to treat microscopically percolating fluid regions differently from microscopically non-percolating regions. With our definitions, the microscopically percolating liquid region corresponds to the main cluster whereas the isolated clusters correspond to the microscopically non-percolating regions. As in a series of previous works [START_REF] Ahmad | From micro-scale to macroscale modeling of solute transport in drying capillary porous media[END_REF][START_REF] Moghaddam | Evaporation in capillary porous media at the perfect piston-like invasion limit: Evidence of nonlocal equilibrium effects[END_REF][START_REF] Moghaddam | Kinematics in a slowly drying porous medium: Reconciliation of pore network simulations and continuum modeling[END_REF][START_REF] Ahmad | Non-local equilibrium continuum modeling of partially saturated drying porous media: Comparison with pore network simulations[END_REF], the method to discuss the relevance of the continuum model is to proceed via comparisons with simulations with a pore network model (PNM) of drying. However, since PNMs have been developed and a PNM is used here as a reference, one can wonder why continuum models are still worth of interest. Due to computational issues, PNM simulations are actually limited to small networks. The largest network considered so far in drying PNM simulations is 80 × 80× 80 [START_REF] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[END_REF]. As reported in [START_REF] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[END_REF], the simulation of the full drying for a single realization required approximately 18 h on an Intel Xeon 3.6 GHz workstation. Simulations on larger networks require significantly more time (e.g., 100×100×100 pore networks would require several days per realization [START_REF] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[END_REF]). As an example, consider a mean distance between pores of 10 µm. Then, a 100×100×100 pore network would correspond to 1 mm 3 of porous medium. It is clear that in most applications the computational domain of interest is much larger. In other words, continuum models are usually the only option to deal with a drying situation of interest in the applications because of their much greater computational efficiency. As an illustration of this major difference in term of computational efficiency between continuum models and PNM simulations, the continuum model solutions presented later in the paper are obtained from an Excel sheet in less than one second whereas the corresponding PNM simulations take a few days because the PNM results are averaged over several network realizations. Even when the spatial domain of interest is not very large like in fuel cells for instance, e.g. [START_REF] Medici | Understanding water transport in polymer electrolyte fuel cells using coupled continuum and porenetwork models[END_REF][START_REF] Belgacem | Coupled continuum and condensation-evaporation pore network model of the cathode in polymer-electrolyte fuel cell[END_REF], a full PNM approach is not possible due to the high pore size contrast between the various porous layers in the cell. Relying on continuum models is then again necessary either according to a full continuum approach, e.g. [START_REF] Nandjou | A pseudo-3D model to investigate heat and water transport in large area PEM fuel cells-Part 1: Model development and validation[END_REF], or hybrid approaches combining PNM and continuum modelling, e.g. [START_REF] Medici | Understanding water transport in polymer electrolyte fuel cells using coupled continuum and porenetwork models[END_REF][START_REF] Belgacem | Coupled continuum and condensation-evaporation pore network model of the cathode in polymer-electrolyte fuel cell[END_REF]. In brief, developing better continuum models is more than ever highly desirable due their high computational efficiency and extending their capabilities is also highly desirable to deal with complex drying situations involving dissolved species and / or particles. In the present paper, we consider the archetypical situation sketched in Fig. 2 where only the sample top surface is in contact with the external air. As an additional illustration of the computational efficiency issue between PNM and continuum model simulations, it can be noted that the continuum model formulation for the situation sketched in Fig. 2 is 1D whereas the pore network simulations are 3D. We focus on the situations where the evaporation rate is sufficiently low for the temperature variations to be negligible. This situation is frequently encountered in the laboratory experiments with water at the room temperature and referred to as "isothermal" drying. The corresponding drying process is often described in three main periods [START_REF] Van Brakel | Mass transfer in convective drying[END_REF]: the constant rate period (CRP), the falling rate period (FRP) and the receding front period (RFP) but it is also usual to rather consider drying as a two stage process [START_REF] Brutsaert | Desorption and the two stages of drying of natural tallgrass prairie[END_REF], with stage 1 corresponding to the CRP and stage 2 to the combination of the FPR and RFP. As in [START_REF] Ahmad | From micro-scale to macroscale modeling of solute transport in drying capillary porous media[END_REF], we mainly focus on the period where the main cluster spans the porous sample, i.e. stage 1.

The paper is organized as follows: In Section II, the NLE two equation model is briefly recalled and the three equations continuum model is described. The drying PNM is summarized in Section III. Results of PNM simulations are presented in Section IV. A comparison between PNM results and the three equation continuum model solution is presented in Section V. The model is extended so as to consider the presence of a solute in Section VI where comparisons between the continuum model and PNM simulations are also presented. This is followed by Section VII which proposes a discussion. Section VIII consists of the main conclusions of the study.

II. CONTINUUM MODEL

A. NLE two equation continuum model

The non-local equilibrium (NLE) two-equation model is first recalled. As discussed in some details in [START_REF] Moghaddam | Evaporation in capillary porous media at the perfect piston-like invasion limit: Evidence of nonlocal equilibrium effects[END_REF], the NLE two-equation model is preferred to the conventional approach, i.e. [START_REF] Pel | Analytic solution for the non-linear drying problem[END_REF],

for the following reason. The conventional approach is based on the use of the equilibrium desorption isotherm to relate the vapor partial pressure and the saturation whereas the impact of adsorption phenomena is presumably negligible in the relatively big pores of capillary porous media. This questionable aspect is circumvented by the NLE two equation model whose main variables are the saturation and the vapor partial pressure and which does not rely on the desorption isotherm. As described in [START_REF] Ahmad | Non-local equilibrium continuum modeling of partially saturated drying porous media: Comparison with pore network simulations[END_REF], (

) 2
where ε, t, Deff, S, 𝜌 𝑙 and 𝑃 𝑣 denote the porosity, time, effective vapor diffusivity, liquid saturation, water density and water vapor partial pressure, respectively. 𝑀 𝑣 , R and T represent the molar mass of water, universal gas constant, and temperature;𝑚 ˙ is the liquid-vapor phase change rate. The latter is also referred to as the NLE phase change term. As shown in [START_REF] Ahmad | Non-local equilibrium continuum modeling of partially saturated drying porous media: Comparison with pore network simulations[END_REF], 𝑚 ˙ can be expressed as

𝑚 ˙≈ -𝑎 𝑔𝑙 𝑀 𝑣 𝑅𝑇 𝛽(𝑃 𝑣𝑠 -𝑃 𝑣 ) (3)
where agl is the specific interfacial area and 𝛽 is a coefficient, Pvs is the saturation vapor pressure since adsorption phenomena are not considered.

B. Three equation continuum model

The main novelty with the three equation model is to explicitly consider that the liquid phase can be split in the percolating liquid phase, also referred to as the main cluster, and the non-percolating liquid phase. In the three equation model, the liquid saturation is thus expressed as

𝑆 = 𝑆 1 + 𝑆 2 (4)
where subscript 1 is for the percolating liquid phase, i.e. the main cluster, and subscript 2 for the non-percolating liquid phase, i.e. the isolated clusters. Subscript 3 refers to the vapor (the gas phase is a binary gas formed by air and the vapor of the evaporating species).

Assuming a homogeneous porous medium, mass balance equations for the percolating and nonpercolating liquid phases are expressed as

𝜀𝜌 𝑙 𝜕𝑆 1 𝜕𝑡 + 𝛻. (𝜌 𝑙 𝑈 𝑙1 ) = -𝑚 ˙12 -𝑚 ˙13 (5)
𝜀𝜌 𝑙 𝜕𝑆 2 𝜕𝑡 + 𝛻. (𝜌 𝑙 𝑈 𝑙2 ) = -𝑚 ˙21 -𝑚 ˙23 (6)
where𝜀 is the porous medium porosity, 𝜌 𝑙 is the liquid density, 𝑈 𝑙1 is the filtration velocity in the percolating liquid phase, 𝑈 𝑙2 is the filtration velocity in the non-percolating liquid phase,𝑚 ˙12

is the mass transfer rate between phase 1 and phase 2, 𝑚 ˙13 is the evaporation rate of phase 1 per unit volume of porous medium, 𝑚 ˙21 is the mass transfer rate between phase 2 and phase 1, 𝑚 ˙23 is the evaporation rate of phase 2 per unit volume of porous medium. Since a new isolated cluster actually forms as the result of the fragmentation of the main cluster, we also have,

𝑚 ˙12 = -𝑚 ˙21 (7)
It can be noted that in isothermal drying an isolated cluster cannot reconnect to the main cluster, hence 𝑚 ˙12 > 0.

The gas phase forms a single cluster in the drying process.

where

𝑚 ˙= 𝑚 ˙13 + 𝑚 ˙23 (9)
The phase change rate is expressed as for the NLE two equation continuum model as

𝑚 ˙= 𝑎 𝑙𝑔 𝛽 𝑀 𝑣 𝑅𝑇 (𝑃 𝑣𝑠 -𝑃 𝑣 ) (10)
with

𝑚 ˙13 = 𝑎 𝑙1𝑔 𝛽 𝑀 𝑣 𝑅𝑇 (𝑃 𝑣𝑠 -𝑃 𝑣 ) (11)
𝑚 ˙23 = 𝑎 𝑙2𝑔 𝛽 𝑀 𝑣 𝑅𝑇 (𝑃 𝑣𝑠 -𝑃 𝑣 ) (12)
Where𝑎 𝑙1𝑔 respectively) is the specific interfacial area between phase 1 (phase 2 respectively) and the gas phase. It can be noted that 𝑎 𝑙𝑔 = 𝑎 𝑙1𝑔 + 𝑎 𝑙2𝑔 .

Following [START_REF] Doster | Numerical solutions of a generalized theory for macroscopic capillarity[END_REF], the mass transfer rate 𝑚 ˙12 between the percolating and non-percolating liquid phases is expressed as

-𝑚 ˙12 = 𝜂𝜀𝜌 𝑙 ( 𝑆 2 -𝑆 𝑖𝑟𝑟 𝑆 𝑖𝑟𝑟 -𝑆 ) 𝜕𝑆 𝜕𝑡 ( 13
)
where 𝜂 is a numerical factor, 𝑆 𝑖𝑟𝑟 is the irreducible saturation. Eq.( 13) represents the continuum modeling of the fragmentation process. It expresses that the occurrence of new clusters from the main cluster is proportional to the liquid phase saturation variation rate whereas the factor

( 𝑆 2 -𝑆 𝑖𝑟𝑟 𝑆 𝑖𝑟𝑟 -𝑆
) is qualitatively consistent with the fact that the closer the system to the irreducible saturation, the greater is the probability to form a new cluster from the main cluster since the main cluster structure becomes increasingly ramified as the irreducible saturation is approached [START_REF] Wilkinson | Invasion percolation: a new form of percolation theory[END_REF]. This functional form is qualitatively consistent with the nonlinear increasing variation of the number of clusters during the drying process [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF] (also illustrated in Fig. 4 below from PNM simulations). As we shall see, it also leads to quantitatively consistent results.

The above model is simplified by introducing additional assumptions. The flow in the percolating phase is modelled using the generalized Darcy's law,

𝑈 𝑙1 = -𝑘𝑘 𝑟1 𝜇 𝛻𝑃 𝑙1 (14)
where 𝑈 𝑙1 is the percolating liquid phase filtration velocity, 𝑃 𝑙1 is the pressure in the percolating liquid phase, 𝑘 is the medium permeability, 𝑘 𝑟1 is the percolating phase relative permeability, 𝜇 is the liquid viscosity. By introducing the capillary pressure curve Pc1(S1), where Pc1(S1) is the local pressure difference between the gas phase and the percolating liquid phase, Eq.( 5) can be expressed as

𝜀𝜌 𝑙 𝜕𝑆 1 𝜕𝑡 + 𝛻. (𝜌 𝑙 𝐷 𝑙1 (𝑆 1 )𝛻𝑆 1 ) = -𝑚 ˙12 -𝑚 ˙13 (15)
Where

𝐷 𝑙1 (𝑆 1 ) = -𝑘𝑘 𝑟1 𝜇 𝑑𝑃 𝑐1 𝑑𝑠 1 (16)
The non-percolating liquid phase is assumed immobile. Eq.( 6) is thus simplified as

𝜀𝜌 𝑙 𝜕𝑆 2 𝜕𝑡 = -𝑚 ˙21 -𝑚 ˙23 (17)

C. Boundary conditions

For the 1D case sketched in Fig. 2a, the boundary conditions at the bottom read, -𝜌 𝑙 𝐷 𝑙1 (𝑆 1 )𝛻𝑆 1 . 𝑛 = 0

and

-𝜀(1 -𝑆)𝐷 𝑒𝑓𝑓 𝑀 𝑣 𝑅𝑇 𝛻𝑃 𝑣 . 𝑛 = 0 (19)
where𝑛 is a unit normal vector to the considered surface.

The definition of the boundary conditions at the top surface is still a debated question in the drying theory [START_REF] Ahmad | Non-local equilibrium continuum modeling of partially saturated drying porous media: Comparison with pore network simulations[END_REF][START_REF] Moghaddam | A pore network study of evaporation from the surface of a drying non-hygroscopic porous medium[END_REF][START_REF] Talbi | About Schlünder's model: A numerical study of evaporation from partially wet surfaces[END_REF][START_REF] Talbi | Coupling between internal and external mass transfer during stage 1 evaporation in capillary porous media: interfacial resistance approach, submitted to[END_REF]. In other words, the coupling at the surface between the transport phenomena in the porous medium and in the external air in contact with the porous medium surface is not yet sufficiently well understood and modelled. In the case of the three equation continuum model, the detailed modelling of the mass transfer at the top surface is particularly challenging since one has to consider that the vapor leaves the porous medium surface from three categories of pores: i) by evaporation from liquid filled pores belonging to the main cluster (percolating liquid phase), ii) by evaporation from liquid filled pores belonging to isolated clusters connected to the surface (non-percolating liquid phase), iii) gaseous pores. However, this challenging modelling issue is left for a future work. In the present study, the focus is on the modelling of the liquid phase internal fragmentation. A simplified approach is therefore adopted. As explained in more details later in the paper, this simplified approach consists in assuming that the overall evaporation rate is an input data for the model and not an outcome.

D Macroscopic parameters

As can be seen from the above, the use of the three equation continuum model implies in general to specify several parameters, most of them being non-linear functions of saturation. However, the focus in this paper is on the liquid phase fragmentation process, which is taken into account in the three equation continuum model via Eq.( 13). For this reason, a particular drying regime is considered (as explained in the next section). The consideration of this particular regime allows focusing on 𝑚 ˙12 (Eq.( 13)) without the need to determine the many other parameters of the three-equation continuum model.

Nevertheless, the complete model was presented for the readers interested in more complex situations requiring the consideration of the full set of equations.

III. PORE NETWORK MODEL

As in [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF], a simple cubic network is considered (Fig. 2b). The distance between two adjacent nodes in the network is the lattice spacing, denoted by a. In this model, the pore bodies located at the nodes of the cubic grid are cubes of size dp with dp varying in the range [0.675, 0.725] according to a uniform probability distribution function, noting that the lengths in the PNM are made dimensionless using the lattice spacing a as reference length. The pore throats are channels connecting the pore bodies. The throat size dt is distributed in the range [0.075, 0.125] according to a uniform distribution law. The drying algorithm is the one presented in [START_REF] Prat | Percolation model of drying under isothermal conditions in porous media[END_REF]. As discussed in [START_REF] Prat | Recent advances in pore-scale models for drying of porous media[END_REF][START_REF] Prat | Pore network models of drying, contact angle and films flows[END_REF] or [START_REF] Metzger | Pore-network models: A powerful tool to study drying at the pore level and understand the influence of structure on drying kinetics[END_REF] this algorithm applies to the isothermal drying situation where capillary effects are dominant and corner film flows [START_REF] Prat | On the influence of pore shape, contact angle and film flows on drying of capillary porous media[END_REF][START_REF] Chauvet | Three periods of drying of a single square capillary tube[END_REF] can be neglected. The interested readers can refer to the afore-mentioned articles for details on the algorithm and additional information on the pore network modelling of the drying process.

It must be recalled that that the viscous effects are not explicitly considered in the liquid phase in the version of the algorithm considered in the present paper. As discussed in [START_REF] Prat | Recent advances in pore-scale models for drying of porous media[END_REF], various drying regimes can be actually distinguished depending on the competition between the capillary forces, the gravity and the viscous forces. Here, the capillarity dominant regime is considered. Both the effects of gravity and the viscous forces are assumed to be negligible compared to the capillarity. This regime is referred to as the capillary regime. As a matter of fact, the special case when the viscous effects can be neglected compared to the capillary effects even when the main cluster becomes very ramified as the irreducible saturation is approached is considered. For this reason, this regime is referred to as the "asymptotic" capillary regime.

Additional details on this regime are given later in the paper.

IV. PORE NETWORK SIMULATIONS

PNM simulations of the drying process were performed with a N × N×N cubic network, where Fig. 3a shows the variation of the computed evaporation rate (normalized by the evaporation rate at t =0) as a function of Snet, i.e. the overall network liquid saturation. Note that these data as well as the other PNM data presented in the paper are averages over 15 realizations of the network unless otherwise mentioned. The classical evolution [START_REF] Van Brakel | Mass transfer in convective drying[END_REF][START_REF] Brutsaert | Desorption and the two stages of drying of natural tallgrass prairie[END_REF] is retrieved with a first period, referred to as stage 1, in which the evaporation rate varies weakly over a significant range of saturations (corresponding to Snet varying between about 0.8 and 0.31). This period is referred to as stage 1. As can be seen in Fig. 3a, stage 1 ends when Snet is about equal to 0.31, which corresponds to the vertical dashed line in Fig. 3a Then the evaporation rate drops. This corresponds to stage 2. However, contrary to the classical experimental results [START_REF] Van Brakel | Mass transfer in convective drying[END_REF], an initial period where the evaporation rate drops can be observed before the quasi constant rate period starts in stage 1 (this approximately corresponds to Snet in the range [0.9 -1] in Fig. 3a). This period is discussed below with the results on the saturation profiles. This choice is motivated by the fact that the computations of macroscopic parameters from PNM simulations, such as for instance the local porosity as illustrated in Fig. 6, indicate that an averaging volume of size 10a can be considered as a reasonable Representative Elementary Volume (REV). Thus a running averaging procedure is used considering slices of size N × N × 10. The computed values are affected to the center of the slices. It can be also noted that slice saturations are actually typically considered in the experiments [START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF][START_REF] Thiery | Drying regimes in homogeneous porous media from macro-to nanoscale[END_REF]. The selected slice thickness explains why the profiles start at z = 5 (in lattice spacing unit) and ends at z = 25 in Fig. 5. The slice averaged saturations are denoted by Ssl, S1sl and S2sl, respectively.

N
The liquid phase saturation profiles in Fig. 5 present several differences compared to the profiles typically obtained in the experiments. For the capillary regime considered in the present paper, the latter are typically flat, e.g. [START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF][START_REF] Thiery | Drying regimes in homogeneous porous media from macro-to nanoscale[END_REF]. By contrast, edge effects are noticeable at both ends of each profile in Fig. 5. The profiles are flat only in the region away from the edges. Also, they become flat in this region only when Snet is sufficiently low, i.e. when Snet ~ 0.6 -0.7. These differences with the experimental profiles are discussed in detail in [START_REF] Talbi | Coupling between internal and external mass transfer during stage 1 evaporation in capillary porous media: interfacial resistance approach, submitted to[END_REF]. The initial period where the profiles are not flat in the central region is associated with a finite size effect and corresponds to the period where Snet varies from 1 to SBT, where SBT is the saturation at breakthrough. Since the saturation at breakthrough, i.e. when the gas phase reaches for the first time the network bottom, scales as 1 -𝑆 𝐵𝑇 ∝ 𝑁 -𝛼 where  = 0.48 in 3D according to the percolation theory [START_REF] Wilkinson | Invasion percolation: a new form of percolation theory[END_REF][START_REF] Stauffer | Introduction to Percolation Theory[END_REF], this initial period becomes negligible for a sufficiently large network, and thus cannot be seen typically in the experiments. By contrast, the edge effect size is found to be independent of the network size and on the order of a few lattice spacing [START_REF] Talbi | Coupling between internal and external mass transfer during stage 1 evaporation in capillary porous media: interfacial resistance approach, submitted to[END_REF].

As a result, the corresponding variations of the saturation in the edge regions is indiscernible in the experiments since the size of the edge effect regions is typically very small compared to the sample size in the experiments. For these reasons, the main objective of continuum models as regards the saturation profiles should be to predict the flat profile evolutions since the finite size effect impacted period is negligible in most experiments as well as the relative extension of the edge effect regions. For the sake of comparison with the continuum model, we therefore As can be seen, the percolating phase saturation decreases over stage 1 whereas the nonpercolating liquid phase saturation increases. The variations of both saturation is much faster and more important as the end of stage 1 is approached. As discussed in [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF], the end of stage 1 corresponds to the situation when the percolating liquid phase is about to cease percolating. In the considered asymptotic capillary regime, this corresponds to the situation when the irreducible saturation is about to be reached. From Fig. 3a, it can therefore be considered that Sirr ≈ 0.31.

V. CONTINUUM MODEL SOLUTION

In order to solve the three equation continuum model, Eqs. (4-17), the following parameters must be in principle determined: 𝜀, 𝐷 𝑒𝑓𝑓 (𝑆), 𝑎 𝑙1𝑔 , 𝑎 𝑙2𝑔 , 𝛽,𝜂,𝑆 𝑖𝑟𝑟 ,𝑘, 𝑘 𝑟1 , 𝑃 𝑐1 (𝑆 1 ) noting that the parameters 𝐷 𝑒𝑓𝑓 , 𝑎 𝑙1𝑔 , 𝑎 𝑙2𝑔 , , 𝑘 𝑟1 , 𝑃 𝑐1 are non-linear functions of the saturation. As exemplified in [START_REF] Moghaddam | Evaporation in capillary porous media at the perfect piston-like invasion limit: Evidence of nonlocal equilibrium effects[END_REF][START_REF] Ahmad | Non-local equilibrium continuum modeling of partially saturated drying porous media: Comparison with pore network simulations[END_REF][START_REF] Lu | Transport parameters of macroscopic continuum model determined from discrete pore network simulations of drying porous media[END_REF], these parameters can be determined from PNM simulations. Also, because of the non-linearity, Eqs.(4-17) must be generally solved using a numerical method. However, it can be observed from Fig. 5 that the profiles are spatially uniform (when as discussed before the edge effects are not taken into consideration and the initial period affected by the finite size effect is discarded). Under these circumstances, the numerical solution can be greatly simplified and the determination of most of the three equation continuum model parameters from specific PNM numerical simulations can be avoided. In [START_REF] Talbi | Coupling between internal and external mass transfer during stage 1 evaporation in capillary porous media: interfacial resistance approach, submitted to[END_REF], it is shown that Sbulk can be determined from the following equation,

𝜀𝜌 𝑙 𝐻 𝑝𝑚 𝜕𝑆 𝑏𝑢𝑙𝑘 𝜕𝑡 = -𝑗 ( 20
)
where 𝑗 is the evaporation flux and 𝐻 𝑝𝑚 is the porous medium height (Fig. 2).

Since 𝑆 𝑏𝑢𝑙𝑘 = 𝑆 1𝑏𝑢𝑙𝑘 + 𝑆 2𝑏𝑢𝑙𝑘 , one obtains from Eq.( 20)

𝜀𝜌 𝑙 𝐻 𝑝𝑚 𝜕𝑆 1𝑏𝑢𝑙𝑘 𝜕𝑡 = -𝑗 -𝜀𝜌 𝑙 𝐻 𝑝𝑚 𝜕𝑆 2𝑏𝑢𝑙𝑘 𝜕𝑡 (21)
The governing equation for 𝑆 2 , namely Eq.( 17), can be simplified in the bulk as

𝜀𝜌 𝑙 𝜕𝑆 2𝑏𝑢𝑙𝑘 𝜕𝑡 = -𝑚 ˙21 (22)
because, as illustrated in Fig. 8, the liquid-vapor phase change is negligible in the network outside the top edge effect region.

This leads to express Eq.( 21) as

𝜀𝜌 𝑙 𝜕𝑆 1𝑏𝑢𝑙𝑘 𝜕𝑡 = -𝑗 𝐻 𝑝𝑚 + 𝑚 ˙21 (23)
Where, as mentioned before, it is proposed to express 𝑚 ˙21 as We have tested Eqs. [START_REF] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[END_REF][START_REF] Medici | Understanding water transport in polymer electrolyte fuel cells using coupled continuum and porenetwork models[END_REF][START_REF] Belgacem | Coupled continuum and condensation-evaporation pore network model of the cathode in polymer-electrolyte fuel cell[END_REF] from the PNM data. Eqs.( 22) and ( 23) together with Eq.( 24) were solved using a first order finite difference scheme to express the derivatives with respect to time, i.e. expressions of the form

𝑚 ˙21 = -𝑚 ˙12 = 𝜂𝜀𝜌 𝑙 ( 𝑆 2𝑏𝑢𝑙𝑘 -𝑆 𝑖𝑟𝑟 𝑆 𝑖𝑟𝑟 -𝑆 𝑏𝑢𝑙𝑘
𝜕𝑆 2𝑏𝑢𝑙𝑘 𝜕𝑡 = 𝑆 2𝑏𝑢𝑙𝑘 (𝑡+𝛿𝑡)-𝑆 2𝑏𝑢𝑙𝑘 (𝑡)

𝛿𝑡

.The method was explicit, i.e.

the values of the saturations involved in the expression of the source term 𝑚 ˙12 were taken at the previous time step. The saturation 𝑆 𝑏𝑢𝑙𝑘 was obtained from 𝑆 𝑏𝑢𝑙𝑘 = 𝑆 1𝑏𝑢𝑙𝑘 + 𝑆 2𝑏𝑢𝑙𝑘 . Since the initial period affected by the finite size effect is not considered, the simulations started with the following initial conditions imported from the PNM simulations: 𝑆 𝑏𝑢𝑙𝑘 =0.6833,𝑆 1𝑏𝑢𝑙𝑘 =0.678, 𝑆 2𝑏𝑢𝑙𝑘 =0.00511. Since the focus is on the evolution of the percolating and non-percolating phases, the evaporation rate was considered as an input data. Thus, the evaporation flux j computed from the PNM simulations was used. This led to the results depicted in Fig. 7, which were obtained with 𝜂 = 0.1and 𝑆 𝑖𝑟𝑟 = 0.3138. As can be seen, the three equation continuum model leads to a quite reasonable agreement with the PNM data. This is an interesting confirmation in the context of drying of the approach proposed in [START_REF] Doster | Numerical solutions of a generalized theory for macroscopic capillarity[END_REF].

VI. SOLUTE CONCENTRATION EVOLUTION IN THE HYPERDIFFUSIVE LIMIT

As pointed out in [START_REF] Ahmad | From micro-scale to macroscale modeling of solute transport in drying capillary porous media[END_REF] or in [START_REF] Diouf | Locus of first crystals on the evaporative surface of a vertically textured porous medium[END_REF], the evolution of the concentration of a solute in a drying porous medium generally results from two main effects: i) the solute convective transport in the percolating liquid phase which leads to the accumulation of solute in the porous medium top surface region (for the configuration depicted in Fig. 2), ii) the fact that the volume occupied by the liquid phase shrinks during drying whereas the total amount of solute in this volume does not change (precipitation or wall deposit phenomena being assumed negligible). It was shown in [START_REF] Ahmad | From micro-scale to macroscale modeling of solute transport in drying capillary porous media[END_REF] that the commonly used continuum model of solute transport [START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF][START_REF] Guglielmini | Drying of salt solutions in porous materials: Intermediate-time dynamics and efflorescence[END_REF][START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF] did not lead to a good agreement with PNM simulations. This was attributed to the fact that the classical approach does not distinguish between the percolating liquid phase and the non-percolating liquid phase. The objective in what follows is to explore whether the three equation continuum model can help alleviating the discrepancies between the PNM simulations and the continuum approach. To this end, we focus on the second mechanism as regards the variation of the solute concentration. This mechanism is referred to as the liquid phase shrinking effect since the increase in the concentration due to this mechanism results from the decrease in the volume occupied by the liquid phase in the porous medium. We only consider the limiting situation where convective effect on the solute transport can be neglected, i.e. the very low Péclet number case [START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF]. This corresponds for instance to a very low evaporation rate. This limiting situation is referred to as the hyperdiffusive limit since the solute concentration is consistently assumed to be spatially uniform in each liquid cluster. Therefore, in this limit, there is no need to compute the velocity field in the liquid phase. The solution is dilute so that impact of the solute on the surface tension, the liquid density or the equilibrium vapor pressure at the menisci can be neglected.

A. PNM computations

Initially, the concentration is uniform in the liquid phase and denoted by C0. Then the concentration is updated in each cluster according to the mass conservation equation, 𝐶 𝑖 (𝑡 + 𝛿𝑡)𝑉 𝑖 (𝑡 + 𝛿𝑡) = 𝐶 𝑖 (𝑡)𝑉 𝑖 (𝑡) [START_REF] Nandjou | A pseudo-3D model to investigate heat and water transport in large area PEM fuel cells-Part 1: Model development and validation[END_REF] where 𝐶 𝑖 is the concentration in cluster # i and 𝑉 𝑖 is the volume of liquid cluster # i. The time step 𝛿𝑡 in Eq.( 25) is the time step of the PNM drying algorithm. Eq.( 25) is for a shrinking cluster. A cluster can also split into two smaller clusters as the result of the invasion of a pore by the gas phase. In this case, the concentration 𝐶 𝑖1 𝐶 𝑖2 = 𝐶 𝑖1 𝑖 ⁄ 2 in the newly formed cluster i1 and i2 is computed from the equation,

𝐶
𝐶 2𝑠𝑙 (𝑡, 𝑧) = ∑ 𝐶 2𝑖 (𝑡)𝑉 𝑖𝑠𝑙 (𝑡) 𝑖=𝑛 𝑖=2 ∑ 𝑉 𝑖𝑠𝑙 (𝑡) 𝑖=𝑛 𝑖=2 = ∑ 𝐶 2𝑖 (𝑡)𝑉 𝑖𝑠𝑙 (𝑡) 𝑖=𝑛 𝑖=2 𝜀𝑆 2𝑠𝑙 𝐴ℎ ( 28
)
where𝑉 𝑖𝑠𝑙 is the volume of cluster #i present in the considered slice, ℎ = 10𝑎 is the slice thickness. It is recalled that label #1 is for the percolating cluster. The concentration in the later is spatially uniform over the whole network. Since 𝐶 1 is uniform of the network, it is obvious that the slice averaged concentration𝐶 1𝑠𝑙 is spatially uniform and equal to 𝐶 1𝑛𝑒𝑡 .

The evolution of the three slice average concentration profiles, namely Csl, C1sl and 𝐶 2𝑠𝑙 , during stage 1 is depicted in Fig. 9. As expected the concentration increases during stage 1 because of the liquid cluster dynamics, i.e. the fact that clusters split and shrink. However, it can be noticed that the fragmentation process occurs essentially in the main cluster, i.e. the cluster forming the percolating liquid phase, when only the bulk region is considered since the evaporation rate of the isolated clusters located within the bulk region is quite small (Fig. 8). Interestingly, the concentration profiles are flat. This was of course expected for C1sl since the concentration is uniform in the percolating liquid phase but somewhat less obvious as regards the nonpercolating liquid phase slice averaged concentration, i.e.𝐶 2𝑠𝑙 . This point is discussed further below and explained in the discussion section (section VII). Also, the results shown in Fig. 9 indicate that the concentration is greater in the percolating liquid phase than in the non-percolating liquid phase and thus greater than the average concentration in the liquid phase. This is better illustrated in Fig. 10 showing the variations of

Where i = 1 is for the percolating cluster and n is the number of liquid clusters present in the considered slice. Similarly, the average solute concentration in the non-percolating clusters is computed as

𝐶 2𝑏𝑢𝑙𝑘 = ∑ 𝐶 𝑖 (𝑡)𝑉 𝑖𝑠𝑙 (𝑡) 𝑖=𝑛 𝑖=2 ∑ 𝑉 𝑖𝑠𝑙 (𝑡) 𝑖=𝑛 𝑖=2 (30)
Note again that the concentration is uniform over the percolating cluster and therefore the slice concentration 𝐶 1𝑏𝑢𝑙𝑘 is equal to the concentration over the whole percolation cluster.

Fig. 10 clearly illustrates why it is important to distinguish the percolating and non-percolating liquid phases. Not making this distinction, as in the commonly used approach, i.e. [START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF][START_REF] Guglielmini | Drying of salt solutions in porous materials: Intermediate-time dynamics and efflorescence[END_REF][START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF], leads to underestimate the concentration and thus for instance to overestimate the time corresponding to the onset of crystallization when C corresponds to a dissolved salt concentration [START_REF] Hidri | Porous medium coffee ring effect and other factors affecting the first crystallisation time of sodium chloride at the surface of a drying porous medium[END_REF]. Then the question arises as to whether the evolution depicted in Fig. 10 can be captured by an extended version of the three equation continuum model considering also the presence of the solute. This is studied in the section that follows.

B. Continuum approach

Within the framework of the three equation continuum model, the solute transport equation in the percolating liquid phase is expressed as

𝜀 𝜕𝑆 1 𝐶 1 𝜕𝑡 = 𝛻. (𝜀𝑆 1 𝐷 𝑠1 𝛻𝐶 1 ) -𝑚 ˙12𝑠 (31)
Where 𝑚 ˙12𝑠 (kg/m 3 /s) is the solute mass exchange term between the percolating liquid phase and the non-percolating liquid phase. Eq.( 31) is similar to the solute transport equation in the classical approach [START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF][START_REF] Guglielmini | Drying of salt solutions in porous materials: Intermediate-time dynamics and efflorescence[END_REF][START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF]. The difference lies in the fact that the considered liquid phase is the percolating liquid phase and not the whole liquid phase.

For the non-percolating liquid phase, the solute conservation equation is simply expressed as

𝜀 𝜕𝑆 2 𝐶 2 𝜕𝑡 = 𝑚 ˙12𝑠 (32)
since the convective transport is assumed negligible in the isolated clusters.

Integrating Eqs. ( 31) and ( 32) over the porous medium height leads to

𝜀𝐻 𝑝𝑚 𝜕𝑆 1 𝐶 1 𝜕𝑡 = -∫ 𝑚 ˙12𝑠 𝐻 𝑝𝑚 0 𝑑𝑧 (33)
𝜀𝐻 𝑝𝑚 𝜕𝑆 2 𝐶 2 𝜕𝑡 = ∫ 𝑚 ˙12𝑠 𝐻 𝑝𝑚 0 𝑑𝑧 (34)
Since the profiles are approximately spatially uniform in the bulk (as shown in Figs. 5 and9), Eqs. ( 33) and ( 34) can be expressed as

𝜀 𝜕𝑆 1𝑏𝑢𝑙𝑘 𝐶 1𝑏𝑢𝑙𝑘 𝜕𝑡 = -𝑚 ˙12𝑠 (35)
𝜀 𝜕𝑆 2𝑏𝑢𝑙𝑘 𝐶 2𝑏𝑢𝑙𝑘 𝜕𝑡 = 𝑚 ˙12𝑠 (36)
The exchange term 𝑚 ˙12𝑠 in Eqs. [START_REF] Chauvet | Three periods of drying of a single square capillary tube[END_REF][START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF] is modelled using an expression similar to the one for 𝑚 ˙12 since the solute mass transfer between the percolating liquid phase and the non-percolating liquid phase is due to clusters separating from the main cluster. Furthermore, since the solute concentration in the main cluster is 𝐶 1 it is reasonable to consider that 𝑚 ˙12𝑠 should be proportional to 𝐶 1 . This finally leads to express 𝑚 ˙12𝑠 as,

-𝑚 ˙12𝑠 = 𝜂𝜀𝐶 1𝑏𝑢𝑙𝑘 ( 𝑆 2𝑏𝑢𝑙𝑘 -𝑆 𝑖𝑟𝑟 𝑆 𝑖𝑟𝑟 -𝑆 𝑏𝑢𝑙𝑘 ) 𝜕𝑆 𝑏𝑢𝑙𝑘 𝜕𝑡 (37)
Eqs. ( 35) and ( 36) combined with Eq.( 37) were solved using a method similar to the one used for solving Eqs.( 22) and ( 23) with, as for the saturation problem, 𝜂 = 0.1and 𝑆 𝑖𝑟𝑟 = 0.3138. A first order finite difference scheme to express the derivatives with respect to time was used, i.e. expressions of the form

𝜕𝑆 1𝑏𝑢𝑙𝑘 𝐶 1𝑏𝑢𝑙𝑘 𝜕𝑡 = (𝑆| |1𝑏𝑢𝑙𝑘𝐶 1𝑏𝑢𝑙𝑘 )(𝑡 + 𝛿𝑡) - (𝑆| |1𝑏𝑢𝑙𝑘𝐶 1𝑏𝑢𝑙𝑘 )(𝑡)

𝛿𝑡

.The method was explicit, i.e. the values of the saturations and the value of 𝐶 1𝑏𝑢𝑙𝑘 involved in the expression of the source term 𝑚 ˙12𝑠 were taken at the previous time step. Then 𝑆 𝑏𝑢𝑙𝑘 𝐶 𝑏𝑢𝑙𝑘 was obtained from 𝑆 𝑏𝑢𝑙𝑘 𝐶 𝑏𝑢𝑙𝑘 = 𝑆 1𝑏𝑢𝑙𝑘 𝐶 1𝑏𝑢𝑙𝑘 + 𝑆 2𝑏𝑢𝑙𝑘 𝐶 2𝑏𝑢𝑙𝑘 .

This led to the results depicted in Fig. 11. As can be seen, the three equation continuum model leads here again to a quite reasonably good agreement with the PNM data. As explained before, only the range of saturations not significantly affected by the initial finite size effect is considered for the comparison between the continuum model and the PNM data. This corresponds to the overall network saturations lower than approximately Snet =0. As expected the product 𝑆 𝑏𝑢𝑙𝑘 𝐶 𝑏𝑢𝑙𝑘 is constant over the period of interest. This simply means that the total mass of solute in the bulk region is conserved. The total mass of solute in the percolating phase, i.e. 𝑆 1𝑏𝑢𝑙𝑘 𝐶 1𝑏𝑢𝑙𝑘 , decreases. This is because the percolating liquid phase loses mass as the result of the fragmentation process. The corresponding mass loss corresponds to the mass of solute gained by the non-percolating liquid phase, which therefore increases during stage 1. The latter corresponds to 𝑆 2𝑏𝑢𝑙𝑘 𝐶 2𝑏𝑢𝑙𝑘 in Fig. 11.

From the computation of 𝑆 𝑏𝑢𝑙𝑘 𝐶 𝑏𝑢𝑙𝑘 , 𝑆 1𝑏𝑢𝑙𝑘 𝐶 1𝑏𝑢𝑙𝑘 and 𝑆 2𝑏𝑢𝑙𝑘 𝐶 2𝑏𝑢𝑙𝑘 depicted in Fig. 11 . The corresponding results are compared to the PNM simulation data in Fig. 10. As can be seen, the three equation continuum model predicts quite well the variations of the three considered solute concentrations. In particular, the important fact that the solute concentration is greater in the percolating liquid phase is well captured.

VII. DISCUSSION

The classical one equation continuum model [START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF][START_REF] Guglielmini | Drying of salt solutions in porous materials: Intermediate-time dynamics and efflorescence[END_REF][START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF] predicts that the solute concentration is uniform over the liquid phase in the considered very low Péclet number limit during stage 1.

This solute concentration is simply given by

𝐶(𝑡) = 𝑆 𝑛𝑒𝑡0 𝐶 0 𝑆 𝑛𝑒𝑡 (𝑡) (38)
Where 𝑆 𝑛𝑒𝑡0 and 𝐶 0 are the initial saturation and solute concentration in the network respectively. 𝐶(𝑡) corresponds to Cbulk in the results presented in the previous section. The ) 𝐶 𝑖 𝑖=𝑛 𝑖=1 [START_REF] Lu | Transport parameters of macroscopic continuum model determined from discrete pore network simulations of drying porous media[END_REF] where 𝑉 𝑖 is the volume of liquid cluster # i and Ci is the solute concentration in cluster # i, n is the number of liquid clusters at the considered time.

In other words, although diffusion is the dominant transport mechanism in the very low Péclet number limit, this does not mean that the solute concentration is spatially uniform in the liquid phase. To illustrate further this feature, the standard deviation of the solute concentration over ) 𝐶 𝑖 𝑖=𝑛 𝑖=2 [START_REF] Diouf | Locus of first crystals on the evaporative surface of a vertically textured porous medium[END_REF] Noting that i = 1 corresponds to the main cluster. Thus i in the range [2, n] corresponds to the isolated clusters (n-1 is thus the total number of isolated clusters). The standard deviation of the solute concentration over the non-percolating liquid phase is then computed as,

𝜎 𝑐2 = √∑ 𝑉 𝑖 𝑖=𝑛 𝑖=2 (41)
The variation of 𝜎 𝑐2 over stage 1 is shown in Fig. 12.

FIG. 12. Variation of the standard deviation𝜎 𝑐2 relative to the mean of the solute concentration distribution in the non-percolating liquid phase over the network bulk region during stage 1.

The inset shows the distribution of the concentration in the isolated clusters for Snet = 0.4.

As can be seen, the standard deviation relative to the mean increases during most of stage 1.

The decrease toward the end of stage 1 is due to the fact that the mean increases faster than the standard deviation because of the formation of many clusters of higher concentration toward the end of stage 1 (Fig. 4). As illustrated in the inset in Fig. 12, the variation of the solute concentration over the non-percolating liquid phase is significant with about a factor 2 between the concentration in the clusters of lowest concentration and the clusters of highest concentrations.

The existence of the concentration spatial variability in the considered diffusion dominant regime is directly due to the fact that isolated clusters form from the main cluster all along stage 1 (as illustrated in Fig. 4). The concentration in an isolated cluster in the bulk region is the concentration in the main cluster at the time when the isolated cluster forms. Since the concentration in the main cluster increases during stage 1 (Fig. 11), the later an isolated cluster forms in the bulk region, the greater its concentration is. Based on the main cluster concentration variations depicted in Fig. 11, it can be readily inferred that the concentration over the non-percolating phase varies at a given time in the range [𝐶 0 , 𝐶 1 (𝑡)].

However, this spatial variability due to the historicity of the isolated clusters formation from the main cluster does not imply a spatial variation along the network depth. As illustrated in Fig. 9, the mean concentration profiles are flat. This is due to the fact that the probability of forming a new cluster from the main cluster in the bulk region does not depend of the position [START_REF] Moghaddam | Kinematics in a slowly drying porous medium: Reconciliation of pore network simulations and continuum modeling[END_REF].

The fragmentation of the liquid phase in isolated clusters also leads to a somewhat counterintuitive result after stage 1, i.e. in stage 2. As described in [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF], the percolating liquid phase disappears at the end of stage 1 in the considered "asymptotic" regime purely controlled by capillary effects. This regime is described as asymptotic because it is rarely observed in the standard laboratory experiments due to the viscous effects. In most experiments, the liquid phase is actually percolating up to the receding front forming the lower edge of the dry region developing in the network during stage 2. By contrast, in the considered asymptotic regime, the liquid phase is formed by isolated clusters only and the development of the dry zone results from the gradual evaporation of the isolated clusters [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF]. This asymptotic regime is expected when evaporation is quite low and/or with sufficiently thin systems. In the considered very low Péclet number limit, the concentration profiles could be expected to be flat due to the dominant diffusion transport in the liquid phase.

FIG. 13. Evolution of saturation and concentration profiles in stage 2; z = 0 corresponds to the network top surface (Fig. 2) whereas z = 30 to the network bottom limiting surface.

However, since the isolated cluster evaporation dynamics is not spatially uniform during stage 2, the concentration profiles are actually not flat. This is illustrated in Fig. 13. The evaporation rate at the boundary of the isolated clusters is significantly greater for the clusters in contact with the dry zone. In other words, the cluster evaporation rate rapidly decreases with the increasing depth in the network. Since the concentration increase is due to the clusters shrinking, the concentration increases faster in the region where the cluster evaporation rate is greater. This leads to the remarkably non-linear concentration variations depicted in Fig. 13.

Interestingly, the concentration profiles roughly tend to resemble the exponential like profiles typically resulting from the competition between advection and diffusion effects [START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF]. The mechanism leading to the strongly non-linear shape is of course completely different here and is due to the combination of the fragmentation of the liquid phase in isolated clusters and the screening of the evaporation at the boundary of the clusters located deeper in the network. One can refer to [START_REF] Prat | Isothermal drying of non-hygroscopic capillary-porous materials as an invasion percolation process[END_REF] for more details on the screening phenomenon. The screening phenomenon during stage 1 is illustrated in Fig. 8. A similar rapid decrease in the evaporation flux with the distance from the interface between the dry zone and the shrinking liquid cluster zone exists during stage 2.

Since it has been shown that the liquid films can have a strong impact in drying [START_REF] Prat | On the influence of pore shape, contact angle and film flows on drying of capillary porous media[END_REF][START_REF] Chauvet | Three periods of drying of a single square capillary tube[END_REF][START_REF] Yiotis | Effect of liquid films on the drying of porous media[END_REF][START_REF] Yiotis | Drying in porous media with gravity-stabilized fronts: Experimental results[END_REF],

it must be clear that the results presented in this article as regards the solute distribution are for the situations where the impact of the films on the solute distribution is negligible. For instance, it has been shown [START_REF] Sghaier | On the influence of sodium chloride concentration on equilibrium contact angle[END_REF] that the contact angle for an aqueous solution in the presence of a dissolved salt can be relatively high, on the order of 40° or more. With such values of the contact angle, the development of corner films is significantly hampered. By contrast, for significantly lower contact angles, when the liquid films can develop so as to maintain a hydraulic connection between the "isolated" clusters, some solute transport can occur through the films. This should reduce the concentration variation between clusters. In the case of a sufficiently low evaporation rate for the solute distribution to reach a quasi-steady state during drying, solute diffusion in the films might even lead to a uniform concentration distribution all over the liquid phase.

Hence, as for the drying process in general, i.e. [START_REF] Yiotis | Effect of liquid films on the drying of porous media[END_REF], the consideration of liquid films within the framework of the three equation continuum model would deserve to be studied.

Although the whole set of equations of the three equation continuum model was presented, only a rather simple situation, allowing solving the model quite easily without resorting to the determination of the various transport parameters of the model, was considered. In this respect, it would be interesting to extend the present work by developing a numerical procedure allowing solving the full set of equations. This would notably permit to simulate the full drying process and not only stage 1 and to consider other regimes than the capillary regime and the very low Péclet number regime.

Also, the edge regions, especially the top edge region in the considered drying configuration,

were not studied on the ground that this is the bulk region which is of primary interest for evaluating the continuum models. However, the top region can be of special interest for predicting or analyzing certain phenomena, such as the formation of salt efflorescence at the evaporative surface of porous media for instance, e.g. [START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF][START_REF] Diouf | Locus of first crystals on the evaporative surface of a vertically textured porous medium[END_REF][START_REF] Hidri | Porous medium coffee ring effect and other factors affecting the first crystallisation time of sodium chloride at the surface of a drying porous medium[END_REF][START_REF] Eloukabi | Experimental study of the effect of sodium chloride on drying of porous media: The crusty-patchy efflorescence transition[END_REF][START_REF] Desarnaud | Drying of salt contaminated porous media: Effect of primary and secondary nucleation[END_REF]. Thus, some works should be dedicated to the modelling of the transfers in the top edge region in the future.

It can be also noted that a structured cubic pore network was considered in the present paper whereas realistic networks, extracted from microstructure digital images for instance [START_REF] Dong | Pore-network extraction from micro-computerizedtomography images[END_REF], are unstructured. However, this should not affect the main results of the present paper since the fragmentation process is a general feature inherent to the invasion percolation process at the core of the drying process. The main impact should be on the coefficient 𝜂 (Eq. ( 13)) since the value of this coefficient depends on the particular porous medium considered. Considering various networks, i.e. various microstructures, would be interesting to better characterize this coefficient.

VIII. CONCLUSION

A continuum model of drying in capillary porous media, referred to as the three equation continuum model, was presented. Contrary to the commonly used continuum models of drying, the present model makes an explicit distinction between the percolating and non-percolating liquid phases. From the consideration of the frequently encountered capillary regime, it has been shown that the model was able to predict the variations of the saturation in both the percolating and non-percolating liquid phases during the stage 1. Then, the model was extended so as to predict the evolution of the concentration of a solute in both the percolating and nonpercolating liquid phases. Comparisons with data obtained from pore network simulations were quite satisfactory. In particular, the continuum model consistently predicts that the solute concentration is higher in the percolating phase than in the non-percolating phase. However, only the very low Péclet number regime was considered.

Nevertheless, the results presented in this article also clarify the meaning of the solute concentration predicted by the continuum models, which should be considered as an average concentration over the fragmented liquid phase. Counter-intuitively, it has been shown that the concentration is not necessarily uniform over the liquid phase in the considered very low Péclet number regime. As the result of the liquid fragmentation in numerous clusters, spatial fluctuations of the concentration are expected even in the very low Péclet number limit.

In summary, the three equation continuum model allows significantly more accurate predictions of the solute distribution during drying compared to the classical one equation model. In the more general context of the macroscopic theory of biphasic flow in porous media, the present paper presents both an extension to drying, with and without the presence of a solute, of the model proposed in [START_REF] Hilfer | Capillary pressure, hysteresis and residual saturation in porous media[END_REF][START_REF] Hilfer | Macroscopic capillarity without a constitutive capillary pressure function[END_REF][START_REF] Hilfer | Macroscopic capillarity and hysteresis for flow in porous media[END_REF][START_REF] Doster | Numerical solutions of a generalized theory for macroscopic capillarity[END_REF] and an additional validation of this model, in particular as regards the formulation of the mass exchange term between the percolating and non-percolating liquid phases.

  FIG.2. a) Basic drying situation referred to as macroscopically 1D drying, b) Sketch of pore network representation with external diffusive layer on top.

  the non-local equilibrium (NLE) two-equation model can be expressed as

FIG. 3 .

 3 FIG.3. Evolution of : a) the evaporation rate (drying curve); the vertical dashed line indicates the end of stage 1, b) the overall network liquid saturation Snet in the range [0.31-1]. The reference time is the time when Snet = 0.31, which corresponds to the end of stage 1. The vertical line corresponds to Snet = 0.7, which is the maximum saturation considered for the comparison with the continuum model (see text).

  FIG.4. Number of liquid clusters in the network as a function of time up to the end of stage 1.The vertical line corresponds to Snet = 0.7.

FIG. 7 .

 7 FIG.7. Variation of 𝑆 𝑏𝑢𝑙𝑘 , 𝑆 1𝑏𝑢𝑙𝑘 and 𝑆 2𝑏𝑢𝑙𝑘 as a function of time during stage 1. Comparison between the PNM data and the results from the three equation continuum model. The comparison is performed for overall saturations lower than 0.7 for the comparison to be not hampered by the initial significant finite size effect (see text). The overall saturation 0.7 corresponds to the vertical dashed line in the figure.

16 FIG. 8 .

 168 FIG. 8. Variation of evaporation rate from percolating cluster (J1) and non-percolating clusters (J2) in the slices as a function of z for various network saturations (the colors in the inset corresponds to the same overall saturations as in the main figure).

𝑖1 𝑖 ⁄ 2 (FIG. 9 .

 29 FIG.9. Profiles of Csl, C1sl and 𝐶 2𝑠𝑙 in the network corresponding to Snet =0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3; z = 0 corresponds to the network top surface (Fig.2) whereas z = 30 to the network bottom limiting surface.

FIG. 10 .

 10 FIG. 10. Variation of 𝐶 𝑏𝑢𝑙𝑘 , 𝐶 1𝑏𝑢𝑙𝑘 and 𝐶 2𝑏𝑢𝑙𝑘 in the network as a function of time during the stage 1. The reference time is the time at the end of stage 1 defined as the time when Snet =0.31. The vertical dashed line indicates the beginning of the period of interest for the comparison with the three equation continuum model (see text).

  FIG. 11. Variation of 𝑆 𝑏𝑢𝑙𝑘 𝐶 𝑏𝑢𝑙𝑘 , 𝑆 1𝑏𝑢𝑙𝑘 𝐶 1𝑏𝑢𝑙𝑘 and 𝑆 2𝑏𝑢𝑙𝑘 𝐶 2𝑏𝑢𝑙𝑘 in the network as a function of time during stage 1. The reference time is the time at the end of stage 1 defined as the time when Snet =0.31. The vertical dashed line corresponds to the time when Snet = 0.7 (as explained in the text, only the times greater than the time corresponding to Snet = 0.7 are considered for the comparison between the continuum model and the PNM data).

  results presented in the previous section clearly show that the solute concentration predicted by this quite classical model in the very low Péclet number limit should not be understood as the solute concentration all over the liquid. The liquid phase is actually fragmented and the solute concentration can vary from one liquid cluster to another. As a result, the solute concentration computed with the classical model must be interpreted as an average concentration over the various liquid clusters, more exactly as a weighted average concentration where the weights are the cluster volume fractions, i.e.

  the various clusters can be computed from the PNM results. Since the percolating and the non -percolating liquid phases are distinguished with the three equation continuum model, the spatial variability of the solute concentration is illustrated considering the non-percolating liquid phase. The average solute concentration in the non-percolating liquid phase is computed as

  𝐶 𝑏𝑢𝑙𝑘 , 𝐶 1𝑏𝑢𝑙𝑘 and 𝐶 2𝑏𝑢𝑙𝑘 . As for the saturations𝑆 𝑏𝑢𝑙𝑘 , 𝑆 1𝑏𝑢𝑙𝑘 and 𝑆 2𝑏𝑢𝑙𝑘 , 𝐶 𝑏𝑢𝑙𝑘 , 𝐶 1𝑏𝑢𝑙𝑘 and 𝐶 2𝑏𝑢𝑙𝑘 are the averaged concentrations over the slice located at z = 16. 𝐶 𝑏𝑢𝑙𝑘 is computed as,

	𝐶 𝑏𝑢𝑙𝑘 (𝑡) =	∑ 𝑖=𝑛 𝑖=1 ∑ 𝑖=𝑛 𝐶 𝑖 (𝑡)𝑉 𝑖𝑠𝑙 (𝑡) 𝑉 𝑖𝑠𝑙 (𝑡) 𝑖=1

  and the computation of 𝑆 𝑏𝑢𝑙𝑘 , 𝑆 1𝑏𝑢𝑙𝑘 and𝑆 2𝑏𝑢𝑙𝑘 depicted in Fig.7, the solute concentrations 𝐶 𝑏𝑢𝑙𝑘 , 𝐶 1𝑏𝑢𝑙𝑘 and 𝐶 2𝑏𝑢𝑙𝑘 were obtained as 𝐶 𝑏𝑢𝑙𝑘 =

	𝑆 𝑏𝑢𝑙𝑘 𝐶 𝑏𝑢𝑙𝑘 𝑆 𝑏𝑢𝑙𝑘	, 𝐶 1𝑏𝑢𝑙𝑘 =	𝑆 1𝑏𝑢𝑙𝑘 𝐶 1𝑏𝑢𝑙𝑘 𝑆 1𝑏𝑢𝑙𝑘	, and 𝐶 2𝑏𝑢𝑙𝑘 =
	𝑆 2𝑏𝑢𝑙𝑘 𝐶 2𝑏𝑢𝑙𝑘
	𝑆 2𝑏𝑢𝑙𝑘

Acknowledgments

Financial support from joint project''Drycap'' funded by GIP ANR (project16-CE92-0030-01) and DFG (project TS28/10-1) is gratefully acknowledged.