Marouane Talbi

Marc Prat
email: mprat@imft.fr

Coupling between internal and external mass transfer during stage 1 evaporation in capillary porous media: interfacial resistance approach

The coupling boundary condition to be imposed at the evaporative surface of a porous medium is studied from pore network simulations considering the capillary regime. The study highlights the formation of a thin edge effect region of smaller saturation along the evaporative surface. It is shown that this thin region forms in the breakthrough period at the very beginning of the drying process. The size of this region is studied and shown to be not network size dependent. This region is shown to be the locus of a non-local equilibrium effect. The features lead to the consideration of a coupling boundary condition involving an interfacial mass transfer resistance and an external mass transfer resistance. Contrary to previous considerations, it is shown that both resistances depend on the variation of the saturation, i.e. the fluid topology, and the size of the external mass transfer layer, i.e. the mass transfer rate. This is explained by the evolution of the vapor partial pressure distribution at the surface which becomes increasingly heterogeneous during evaporation and depends on both the evolving fluid distribution in the interfacial region and the mass transfer rate. However, the geometric effects due to the configuration of the fluids can be separated from rate effects that arise due to the nonequilibrium mass transport.

I. INTRODUCTION

Predicting evaporation from a porous medium is important in many applications such as the evaporation from soils [START_REF] Or | Advances in soil evaporation physics-A review[END_REF] or the drying steps occurring in many industrial processes, e.g. [START_REF] Mujumdar | Handbook of Industrial Drying[END_REF].

The topic has motivated numerous studies but is still an active research area because predicting the evaporation process is still challenging. The commonly used models actually involve one or more adjustable parameters such as for instance the critical saturation or critical liquid water pressure, a somewhat controversial concept, or the heat and mass transfer coefficients at the 2 porous medium surface [START_REF] Li | Evaluation of model concepts to describe water transport in shallow subsurface soil and across the soil-air interface[END_REF]. For this resaon, it is often considered that the boundary condition at the evaporative surface is still a somewhat unresolved issue, [START_REF] Li | Evaluation of model concepts to describe water transport in shallow subsurface soil and across the soil-air interface[END_REF]. By commonly used models, we mean the classical models based on the continuum approach to porous media, [START_REF] Whitaker | Simultaneous heat, mass, and momentum transfer in porous media: A theory of drying[END_REF]. For the simple situation where the temperature variations can be neglected, such a drying model typically takes the form of a non-linear diffusive equation combining the liquid mass balance equation and the vapor mass balance equation into a single equation, [START_REF] Pel | Analytic solution for the non-linear drying problem[END_REF]. This model is referred to as a local equilibrium (LE) model because it is based on the assumption that the relationship between the liquid saturation and the vapor partial pressure measured under equilibrium conditions can still be used when evaporation takes place. However, this widely used model has been criticized [START_REF] Li | Evaluation of model concepts to describe water transport in shallow subsurface soil and across the soil-air interface[END_REF][START_REF] Lozano | Experimental study and modelling of the water phase change kinetics in soils[END_REF][START_REF] Moghaddam | Evaporation in capillary porous media at the perfect piston-like invasion limit: Evidence of nonlocal equilibrium effects[END_REF] and it has been argued that considering a non-local equilibrium (NLE) model was more relevant. The most common approach to evaluate a drying continuum model is through comparisons with experimental data, [START_REF] Li | Evaluation of model concepts to describe water transport in shallow subsurface soil and across the soil-air interface[END_REF], [START_REF] Chen | A mathematical model of drying processes[END_REF][START_REF] Kaviany | Funicular state in drying of a porous slab[END_REF]. However, an alternative is to proceed via comparisons with macroscopic data obtained from pore scale numerical simulations. This latter approach has notably been developed using pore network model (PNM) simulations [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][START_REF] Ahmad | From micro-scale to macroscale modeling of solute transport in drying capillary porous media[END_REF][START_REF] Lu | Transport parameters of macroscopic continuum model determined from discrete pore network simulations of drying porous media[END_REF]. Naturally, it is then expected that the PNM simulations are "sufficiently" representative of the drying process for the comparison between the continuum models and the PNM simulations to be also insightful as regards the performance of the continuum models with respect to the experiments. In this respect, it can be noted that PNM simulations favorably compare with microfluidic experiments [START_REF] Laurindo | Numerical and experimental network study of evaporation in capillary porous media. Phase distribution[END_REF][START_REF] Vorhauer | Temperature gradient induced double stabilization of the evaporation front within a drying porous medium[END_REF] but have never been quantitatively compared to experiments with classical porous media (glass beads, soil column, etc). Regarding the latter, it has been pointed out that the drying PNMs reproduce quite nicely several important experimental features [16 -18] but this remains a qualitative comparison. There are at least two reasons explaining this situation. Firstly, most of PNM simulations are performed with structured cubic networks whereas the pore network is unstructured in the experiments. A second problem lies in the size of the networks. Consider for instance the experiments with glass beads reported in [START_REF] Kaviany | Funicular state in drying of a porous slab[END_REF]. The size of the particles forming the porous medium was on the order of 0.2 mm. The glass bead rectangular container diameter was 15 cm in width and length and 5 cm deep. Thus there was about 250 particles over the container height and 750 over the width or the length. Thus the size of the network should be on the order of 250 × 750 × 750 to be representative of the experiment. 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], thus significantly smaller. It must actually be noted that most PNM simulations have been performed with even smaller networks because of computational time and memory size issues. The network should be several orders of magnitude larger for being representative of many experiments at the laboratory scale, not to mention the field situations. In brief, the networks are much smaller that the size that would be necessary for performing direct comparisons with experimental results. Although the network size is an issue for performing direct quantitative comparisons with most available experimental results, one of the major objectives of the paper is to show that this is not the case when the objective is to evaluate a continuum model by comparison with PNM simulations provided that the impact of the network size, i.e. the finite size effects, are adequately taken into account.

Another issue lies in the drying regime and the possible temperature variations for the comparison between experiments, PNM simulations and continuum modelling to be insightful.

The regime and the temperate variations, if any, should be comparable. For simplicity, only the situation where the temperature variations are negligible is considered. This situation is referred to as (quasi-) isothermal drying. Various regimes of isothermal drying are identified and discussed in [START_REF] Prat | Recent advances in pore-scale models for drying of porous media[END_REF] depending on the competition between the capillary forces, the gravity forces and the viscous forces, see also [START_REF] Lehmann | Characteristic lengths affecting evaporative drying of porous media[END_REF]. We focus in the present study on presumably the simpler regime, namely the capillary regime, where the capillary forces are dominant. The main features of this regime are presented in section IV. As in many previous works, the archetypical drying situation sketched in Fig. 1 is considered.

The porous medium is homogeneous and fully saturated by a volatile liquid initially. Only the top surface of the porous sample is in contact with the ambient air whereas the other limiting surfaces are sealed. This situation is deemed to correspond to a 1D transfer situation as regards the continuum modelling.

As mentioned previously, relying on PNM simulations to analyze the drying process in conjunction or not with continuum models is not a novelty. In [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF][START_REF] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[END_REF], the focus was on the analysis of the PNM simulations. No attempt was made to compare the PNM simulations with a continuum model. The focus in [START_REF] Yiotis | Coupling between external and internal mass transfer during drying of a porous medium[END_REF] was on the impact of liquid films. Here again no comparison with a continuum model was reported. By contrast, we focus in the current paper on the situation where the effect of the liquid films can be neglected with an emphasis on the continuum modelling of the "PNM drying". Comparisons between a continuum model and PNM simulations are reported in [START_REF] Moghaddam | Evaporation in capillary porous media at the perfect piston-like invasion limit: Evidence of nonlocal equilibrium effects[END_REF] but the emphasis was not on the capillary regime but on the viscous-capillary regime. Furthermore, the parametrization of the mass transfer at the evaporative surface was not addressed. The focus in [START_REF] Moghaddam | Kinematics in a slowly drying porous medium: Reconciliation of pore network simulations and continuum modeling[END_REF] was on the velocity field induced in the liquid phase by the evaporation and thus different than in the present paper. Comparisons between PNM simulations and a continuum model were presented in [START_REF] Lu | Transport parameters of macroscopic continuum model determined from discrete pore network simulations of drying porous media[END_REF] but only as regards the saturation profiles. As in [START_REF] Moghaddam | Evaporation in capillary porous media at the perfect piston-like invasion limit: Evidence of nonlocal equilibrium effects[END_REF], the regime was the capillary-viscous regime whereas the focus in on the capillary regime in what follows. Also, the size of the external mass transfer layer was not varied whereas it appears to be an important factor.

In summary, the objective of the paper is threefold: i) to clarify the procedure for comparing PNM simulations on small networks and continuum models, ii) to explore the coupling between the external mass transfer and the transfer in the porous medium, iii) to perform the comparison between PNM simulations and a continuum model solution. In relation with point ii), the modelling of transfer in the interfacial region between a porous medium and a free fluid has actually been the subject of many studies, mostly as regards the single phase flow modelling, a problem often referred to as the Beavers and Joseph problem, e.g. [START_REF] Valdés-Parada | Diffusive mass transfer between a microporous medium and an homogeneous fluid: Jump boundary conditions[END_REF][START_REF] Ochoa-Tapia | Heat transfer at the boundary between a porous medium and a homogeneous fluid[END_REF] and references therein. Other transfers have also been studied, e.g. [START_REF] Beavers | Boundary conditions at a naturally permeable wall[END_REF][START_REF] Chandesris | Jump conditions and surface-excess quantities at a fluid/porous interface: a multi-scale approach[END_REF]. However, the drying problem considered in the present paper is presumably significantly more complex due to the evolution of the liquid and gas distribution in the interfacial region during drying.

The paper is organized as follows: The continuum model is presented in section II. The PNM drying model is summarized in Section III. The main features of the capillary regime are presented in Section IV. The saturation profiles are analyzed in Section V. The drying kinetics is discussed in Section VI. The non-local equilibrium (NLE) effect is discussed in Section VII.

The coupling between the internal and external transfers is studied in Section VIII. The solution of the continuum model for the considered drying regime is presented in Section IX. This is followed by Section X which proposes a discussion. Section XI consists of the main conclusions of the study.

II. CONTINUUM MODEL

As discussed in previous works, e.g. [START_REF] Ahmad | Non-local equilibrium continuum modeling of partially saturated drying porous media: Comparison with pore network simulations[END_REF]

Where 𝐷 𝑙 (𝑆) is the liquid diffusivity, 𝜀 is the porosity, 𝐷 𝑒𝑓𝑓 is the vapor effective diffusivity;

𝑚 ˙ is the NLE phase change term, which is expressed as [START_REF] Ahmad | Non-local equilibrium continuum modeling of partially saturated drying porous media: Comparison with pore network simulations[END_REF],

𝑚 ˙≈ 𝑎 𝑔𝑙 𝑀 𝑣 𝑅𝑇 𝛽(𝑃 𝑣𝑠 -𝑃 𝑣 ) (3)
where agl is the specific interfacial area between the liquid and gaseous phases and 𝛽 is a coefficient, Pvs is the saturated vapor pressure. The missing conditions to solve this set of equations are the conditions at the top where the coupling between the mass transfer in the porous medium and the mass transfer in the surrounding atmosphere takes place. This coupling is analyzed later in the paper in Section VIII from comparisons with PNM simulations.

Furthermore, this type of model is generally solved numerically, e.g. [START_REF] Chen | A mathematical model of drying processes[END_REF], using a numerical discretization technique. To obtain the numerical solution, the various effective parameters, such as the liquid diffusivity or the vapor effective diffusivity must be specified. These parameters are non-linear functions of the saturation. Our approach is different. In order to focus on the coupling at the top, we wish to avoid the complexity and the approximations associated with a numerical approach and the unavoidable uncertainties introduced by the computations of the effective parameters 𝐷 𝑙 (𝑆) and 𝐷 𝑒𝑓𝑓 on small networks, e.g. [START_REF] Moghaddam | Evaporation in capillary porous media at the perfect piston-like invasion limit: Evidence of nonlocal equilibrium effects[END_REF]. To this end, as mentioned in the introduction, the capillary regime is considered. As we shall see, a simple analytical solution can be obtained in this case without resorting to the specifications of the effective parameters.

The various macroscopic variables in the continuum model, i.e. 𝑃 𝑣 , S, are classically interpreted as spatial averages over an averaging volume (AV), e.g. [START_REF] Whitaker | Simultaneous heat, mass, and momentum transfer in porous media: A theory of drying[END_REF]. This will allow to perform comparisons between the continuum model and the PNM by volume averaging the PNM simulation results. This view of the macroscopic variables allows to clarify the concept of nonlocal equilibrium (NLE). Consider an averaging volume containing some liquid and neglect Kelvin and adsorption phenomena. Then, under equilibrium conditions, the vapor partial pressure is uniform over the averaging volume and equal to the vapor partial pressure at the menisci, i.e. the saturated vapor pressure 𝑃 𝑣𝑠 . As a result, the volume averaged vapor partial pressure is also equal to 𝑃 𝑣𝑠 . In drying, the vapor partial pressure is not spatially uniform in the pore space in the regions where mass transfer occurs in the gaseous phase (vapor partial pressure gradients are generated). As the result, the volume averaged vapor partial pressure can deviate from 𝑃 𝑣𝑠 . When this happens in regions where liquid is present, the deviation is referred to as a non-local equilibrium (NLE) effect. As can be seen from Eq.( 3), the above model is based on a near-equilibrium approximation considering a linear variation of the phase change term with the vapor partial pressure deviation. Actually, it will be shown that the NLE effect is confined in the interfacial region in the capillary regime. This localization of the NLE effect will be harnessed to derive a simple solution to the continuum model without the need to specify parameters agl and 𝛽 in Eq.(3).

III. PORE NETWORK DRYING MODEL

As in many previous works, e.g. [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][START_REF] Lu | Transport parameters of macroscopic continuum model determined from discrete pore network simulations of drying porous media[END_REF][START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF][START_REF] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[END_REF][START_REF] Yiotis | Coupling between external and internal mass transfer during drying of a porous medium[END_REF], a simple cubic network is considered (Fig. 2). 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 of square cross-section connecting the pore bodies. The throat diameter dt is distributed in the range [0.075, 0.125] according to a uniform distribution law. The PNM 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],

this algorithm applies to the isothermal drying situation where capillary effects are dominant and corner film flows [START_REF] Yiotis | Coupling between external and internal mass transfer during drying of a porous medium[END_REF][START_REF] Prat | On the influence of pore shape, contact angle and film flows on drying of capillary porous media[END_REF] can be neglected. This algorithm essentially combines invasion percolation (IP) rules [START_REF] Wilkinson | Invasion percolation: a new form of percolation theory[END_REF] to deal with the capillary effects and a finite-volume like numerical solution of the vapor diffusion problem in the region occupied by the gas phase, considered as a binary mixture of vapor and an inert gas (air).

FIG. 2. Sketch of pore network model with external mass transfer layer on top.

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. As sketched in Fig. 2, the coupling between the internal mass transfer and the external mass transfer is taken into account by setting computational nodes in an external layer of size H, where the vapor transfer is governed by (Fick's) diffusion. More details on this approach can be found 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] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[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]. Additional information on the evaporation flux computation in direct relation with the present paper is given in Appendix A.

As discussed in several previous works, e.g. [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF][START_REF] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[END_REF], the PNM simulations show that the drying process results in the formation of many liquid clusters, i.e. separate groups of interconnected pores and throats filled with liquid, a process referred to as the liquid phase fragmentation process. It is convenient to distinguish the main cluster (MC), from the smaller liquid clusters, also referred to as the isolated clusters. The main cluster is the liquid cluster with the greatest number of pores and throats. Initially the main cluster is percolating whereas the other clusters that form are not percolating. As in classical percolation [START_REF] Stauffer | Introduction to Percolation Theory[END_REF], percolating means that the cluster is connected to both the top and bottom sides of the porous medium. The main cluster gradually shrinks during the drying process, becomes very ramified and eventually ceases to be percolating. In relation with this evolution of the main cluster, one can define the first stage of drying, namely stage 1 evaporation, as the period over which the main cluster is percolating [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF][START_REF] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[END_REF]. As discussed in the next section, stage 1 evaporation is more classically defined considering the drying kinetics rather than the evolution of the main cluster. However, the PNM simulations indicate that the stage 1 evaporation definition based on the MC evolution and the one based of the drying kinetics coincides. In other words, the end of stage 1 evaporation corresponds to the disconnection of the main cluster from the top surface.

IV. COMPARISON BETWEEN EXPERIMENTAL AND PNM RESULTS

In order to perform the comparison between the continuum model and the PNM simulations, a first issue must be addressed. It has to do with the fact that the gas phase forms a percolating cluster as the result of the gradual invasion of pores and throats by the gas. The situation when the gas phase forms for the first time a percolating cluster, i.e. reaches for the first time the bottom of the porous domain, is called the breakthrough (BT). The gas phase at BT forms a percolating cluster, which has the same fractal properties as a percolating cluster at the percolation threshold in classical percolation [START_REF] Lenormand | Invasion percolation in an etched network: measurement of a fractal dimension[END_REF]. On the other hand, the continuum model is based on the classical length scale separation concept, i.e. a representation of the medium as a collection of unit cells whose size is small compared to the porous domain size [START_REF] Battiato | Applicability regimes for macroscopic models of reactive transport in porous media[END_REF]. Such a representation is inadequate in the vicinity of a percolation threshold owing to the fractal nature of the fluid distribution. As a result, the comparison between the continuum model and the PNM simulations makes sense only when the system is sufficiently away from BT for the length scale separation to hold. This issue is addressed and illustrated in what follows starting from some comparisons between experimental results and PNM simulations.

Typical drying experiments aim at measuring i) the cumulative evaporation, i.e. the evaporated mass as a function of time, or equivalently the drying kinetics, i.e. the variation of the evaporation rate as a function of the overall saturation model, ii) the saturation profiles. Fig. 3 show a few representative results from the literature together with PNM results obtained using the algorithm described in [START_REF] Prat | Percolation model of drying under isothermal conditions in porous media[END_REF]. These results are deemed to be representative of the capillary regime.

EXPERIMENT [START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF] PNM (this work) FIG. 3. Top: comparison between typical experimental and drying PNM saturation profiles in the capillary regime. Bottom: comparison between a typical experimental drying kinetics and a PNM drying kinetics. MCD indicates the main cluster disconnection [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF][START_REF] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[END_REF]. BT indicates the breakthrough.

As mentioned earlier, several important features are well captured by the PNM approach: drying can be described in two main periods, referred to as stage 1 and stage 2, the saturation profiles are flat during stage 1 whereas a drying zone develops during stage 2 leading to a drop in the evaporation rate. However, interesting differences are also clearly visible:

1) the experimental saturation profiles are flat in the experiments during stage 1 right from the first measured profiles. Flat profiles are also obtained in the PNM simulations but after a while, i.e. when the overall saturation is significantly lower than in the experiments (S ~ 0.6 -0.7 in Fig. 3).

2) The PNM profiles are not flat in both edge regions, namely the top and bottom ends of the network. Similar edge effects are not visible in the experimental profiles.

3) The evaporation rate decreases sharply at the beginning of the drying process, i.e. at the beginning of stage 1, in the PNM simulations. This is not seen in the experiments.

These various points are discussed in the next sections in order to define the conditions for performing a meaningful comparison between continuum models and PNM simulations on small networks.

V. SATURATION PROFILES

A. Saturation profile flatness and finite size effects

The computation of saturation profiles from PNM simulations requires defining an averaging volume (AV). As in several previous works [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][START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF], horizontal slices of thickness a (where a is the lattice spacing) are considered. Each slice contains N × N pore bodies located in a horizontal plane and the half of the vertical throats connecting the slice pore bodies to the pore bodies in the two adjacent horizontal planes. The thin-slice averaged saturation is thus computed as,

𝑆 𝑡ℎ𝑖𝑛-𝑠𝑙𝑖𝑐𝑒 = ∑ 𝑉 𝑝𝑙𝑗+∑ 𝑉 𝑡𝑙𝑗+ 𝑗=𝑛 𝑡 𝑗=1 ∑ 𝑉 𝑝𝑗 +∑ 𝑉 𝑡𝑗 𝑗=𝑛 𝑡 𝑗=1 𝑗=𝑛 𝑝 𝑗=1 𝑗=𝑛 𝑝 𝑗=1 (6)
where 𝑛 𝑝 is the number of pores in the slice (np = N × N), 𝑉 𝑝𝑙𝑗 is the volume of liquid in pore #j, 𝑉 𝑝𝑗 is the volume of pore #j, 𝑛 𝑡 is the number of throats in the slice, 𝑉 𝑡𝑙𝑗 is the volume of liquid in throat #j within the slice, 𝑉 𝑡𝑗 is the volume of throat #j within the slice. The thin slice averaged saturation profiles were computed in drying for various sizes of the external diffusive layer, namely H = 5, 10, 20, 30, 40, 50 (in lattice spacing unit) for a 30 × 30 ×30 network. The profiles were found to be independent of H. On the one hand, this is not surprising since the main cluster invasion is controlled by the distribution of throat sizes (the sequence of throat selection at the boundary of the main cluster is independent of the evaporation rate in the considered capillary regime and follows IP rules). On the other hand, it will be seen later in the paper than the smaller the external mass transfer layer thickness H, the stronger is the NLE effect in the top region of the network. Thus, an impact on the evolution of the isolated clusters forming in this top region could be expected. The simulations thus show that this impact, if any, is indiscernible on the saturation profiles. The thin-slice saturation profiles are depicted in Fig. 4.

As mentioned with point #1, flat profiles are obtained in the PNM simulations but only after a while, i.e. when S ~ 0.6 -0.7 in Fig. 4. Physically, non-flat profiles are actually expected in the beginning of drying since the invasion of the porous medium by the gas phase as the result of evaporation starts in the top region of the network. At some point, the gas phase reaches the network bottom (breakthrough BT). Since the gas phase cluster at BT is fractal, some network size dependence is expected in this initial period and a while after until the saturation becomes uniform over the height of the porous medium. To analyze further this initial period, it can be first noticed that the saturation profiles obtained in drying and in drainage using the classical invasion percolation (IP) algorithm without trapping [START_REF] Wilkinson | Invasion percolation: a new form of percolation theory[END_REF] are very similar. This is illustrated in As the result, the impact of the network size can be more easily studied using the IP algorithm, which is significantly less computationally demanding compared to the drying algorithm. The IP saturation profiles for two network sizes, namely 20 × 20 × 20 and 100 × 100 × 100 are compared in Fig. 5.

The comparison between the two network sizes in Fig. 5 illustrates the network size dependence of the initial period up to BT and a while after. The profile is getting flat (in the bulk) for a greater network saturation in the largest network. The network size dependence at BT is further illustrated in Fig. 6 with the IP saturation profiles at BT for various network sizes showing that the saturation at BT increases with the network size. Actually, it is known from percolation theory [START_REF] Wilkinson | Invasion percolation: a new form of percolation theory[END_REF][START_REF] Stauffer | Introduction to Percolation Theory[END_REF] that the saturation at BT obeys the following scaling law in 3D,

1 -𝑆 𝐵𝑇 ∝ 𝑁 -0.48 (7)
Accordingly, 𝑆 𝐵𝑇 1 for a sufficiently large network. In other words, the initial period around BT cannot be detected in the typical laboratory experiments which correspond to much larger networks than those considered in the simulations because the saturation variation over the BT period is quite small for large networks. This explains why the saturation profiles are typically flat right from the first measured profiles in the experiments (as illustrated in Fig. 3). It follows that the comparison between continuum models and PNM simulations must exclude the initial period around BT. In other words, the comparison should be made for overall saturations corresponding to "reasonably" flat profiles in the PNM simulations.

B. Edge effects

The second striking difference between the experiments and the PNM simulations lies in the extend of the edge effects. Edge effect refers to the fact that the saturation profiles near the top and bottom ends of the porous medium are not flat. The saturation is lower than in the bulk in the top region whereas this is the opposite in the bottom region where the saturation is greater than in the bulk. In what follows, these edge effects are referred to as the top and bottom edge effects respectively. The edge effects are not really visible in the experimental profiles reported in Fig. 3 [32] because of significant fluctuations around the mean saturation probably due to porous medium heterogeneities in the sample. The experimental profiles reported in [START_REF] Thiery | Drying regimes in homogeneous porous media from macro-to nanoscale[END_REF], see Fig. 1 in [START_REF] Thiery | Drying regimes in homogeneous porous media from macro-to nanoscale[END_REF], are much smoother and almost perfectly flat over almost the full sample height.

The edge effect on the side of the surface open to ambient air (top edge effect) is quite small, if any. On the side of the closed bottom end, a bottom edge effect is visible in Fig. 1 in [START_REF] Thiery | Drying regimes in homogeneous porous media from macro-to nanoscale[END_REF] but the saturation is lower than in the bulk contrary to the PNM simulations. It might be due to some experimental artefact, i.e. the fact that the MRI measurements are impacted by the presence of the wall. Nevertheless, the main observation is that the top edge effect is not visible in the experimental profiles. In other words, the relative extend of the top edge effect region, if any, is less than in the PNM simulations in Fig. 4. Interestingly both edge effects are also visible in the IP profiles displayed in Figs. 4 and5 and quite similar to the ones observed in the drying simulations (Fig. 4). This is an indication that the top edge effect is not primarily due to a stronger evaporation in the network top region in the drying case. Based on the profile evolution displayed in Fig. 4, the top edge effect can be interpreted as a reminiscence of the very first drying period before and around BT, i.e. the BT period. As mentioned before and illustrated in Fig. 4, there is a preferential invasion of the top region before BT is reached. As indicated in the inset in Fig. 4b, almost 80 % of the surface throats are invaded by the gas phase during the BT period. Then the question arises as to whether the edge effect is also network size dependent. A first indication is obtained by plotting in the inset in Fig. 6 the IP slice averaged saturation at the network top at breakthrough as a function of network size. This is the saturation at z = 0 in Fig. 6. Contrary to the saturation deeper in the network, e.g. at z = 0.2 or z = 0.4 in Fig. 6, this saturation varies very weakly, if any, with the network size. In other words, the top edge effect is present regardless of the network size. A similar conclusion holds as regards the bottom edge effect.

Then, the IP thin slice averaged saturation profiles for Snet =0.5 have been plotted in Fig. 7 for various network sizes. This figure also indicates that the extend of the top and bottom edge effect regions is not network size dependent. An estimate of the extend 𝛿of the edge effect regions can be obtained from the plots in Fig. 7. As can be seen, 𝛿 𝑎 ⁄ ≈ 10 both at the bottom and at the top. Because of the evaporation at the network top, one can wonder whether this size is also relevant for the drying case. The comparison of the drying and IP saturation profiles in Fig. 4 does indicate that this is indeed the case. As we shall see, the evaporation is active during stage 1 only in the top edge effect region of size 𝛿. As a result, the saturation is expected to be lower in the drying case in the top edge effect region due to the total or partial evaporation of some isolated clusters in this region but the size of the top edge effect region is actually about the same in drainage (IP profiles) and in drying. In summary, the top edge effect, as well as actually the bottom edge effect, is expected to be present regardless of the network size. The top edge effect is a reminiscence of the invasion period before and around BT whereas the bottom edge effect develops after BT until the profiles becomes essentially flat in the bulk. The BT period is indiscernible in most laboratory experiments because of its network size dependence. In a large network, this period actually corresponds to a very small variation of the overall saturation. Although this initial period cannot be detected in most experiments, it does have an impact regardless of the network size.

The impact is the formation of transition layer at the very top of the network where the liquid distribution and the transfers (as we shall see) are different from the bulk. This transition zone is referred to as the top edge effect region in the paper. This top edge effect region is thin, less than 10 lattice spacing. The existence of this top edge effect region has been overlooked in the modelling of drying according to the continuum approach to porous media and is, as we shall see, a key element in the analysis of the drying process.

VI DRYING KINETICS

The most obvious difference between the drying curve obtained from the PNM simulation and the experimental one in Fig. 3 lies in the evaporation rate sharp decrease at the beginning of the drying process in the PNM simulations. This is not seen in the experiments. In the particular case of the PNM simulation shown in Fig. 3, it can be seen that the mean evaporation rate in most of stage 1 after the BT period is only about 30% of the potential evaporation, defined as the evaporation rate when the network is fully saturated at the very beginning of drying. The initial period of sharp decrease in the evaporation rate corresponds to the "breakthrough period" discussed in section V, i.e. the period up to BT and a while after. Since the top edge effect region forms on top of the sample during the breakthrough period regardless of the network size, the question arises as to whether the 30% decrease in the evaporation rate compared to the potential evaporation is representative of the experiments. To discuss this point, one starts from some typical values of the evaporation flux in experiments at the beginning of drying. For instance, this evaporation flux is ~ 6.6 × 10 -4 kg/m 2 /s in [START_REF] Chandesris | Jump conditions and surface-excess quantities at a fluid/porous interface: a multi-scale approach[END_REF] and ~1.4 × 10 -4 kg/m 2 /s in [START_REF] Wilkinson | Invasion percolation: a new form of percolation theory[END_REF].

Recalling that in slow drying the vapor diffusion transport can be considered as quasi-steady, this flux can be characterized by an equivalent diffusion length, denoted by 𝐻, from the Fickian relationship,

𝑗 = 𝐷 𝑣 𝑀 𝑣 𝑅𝑇 (𝑃 𝑣𝑠 -𝑃 𝑣∞ ) 𝐻 (8) 16
Where 𝐷 𝑣 is the vapor molecular diffusion coefficient in air, 𝑃 𝑣∞ is the vapor partial pressure in the ambience, 𝑃 𝑣𝑠 is the saturated vapor pressure,𝑀 𝑣 , R and T represent the molar mass of the volatile species (water typically in many experiments), universal gas constant, and temperature.

Applied to the experimental data reported in [START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF] and [START_REF] Thiery | Drying regimes in homogeneous porous media from macro-to nanoscale[END_REF], one obtains, 𝐻 ≈ 0.7 mm and 𝐻 ≈ 3.3 mm respectively. These values are representative of many other experiments. Hence, 𝐻 O(1mm). Thus, a first conclusion is that this size is small compared to the size of the sample in the experiments, which is typically on the order of a few cm. On the other hand, the grain / pore size in the experiments is most of the time quite small compared to 𝐻. For instance, the mean pore size in [START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF] is on the order of a few microns. The greater particle size in the experiments reported in [START_REF] Thiery | Drying regimes in homogeneous porous media from macro-to nanoscale[END_REF] is 45 µm. Thus, the experimental situations typically correspond to d << 𝐻

where d is the pore or particle size. Since the top edge effect region is only a few lattice spacing, i.e. particle sizes, thick, its size 𝛿 is typically small compared to the external mass transfer layer thickness, i.e. 𝛿 << 𝐻. In the case of the PNM simulation reported in Fig. 3, the situation is different since 𝐻 = 10𝑎. Thus, 𝐻 is comparable to the top edge effect region extend ≈ 6𝑎 -7𝑎. In other words, the question arises as to whether the initial sharp drop in the evaporation flux in the PNM simulation in Fig. 3 is related to the difference in the ratio 𝐻 / 𝛿 between the experiments, where this ratio is presumably quite large, and the PNM simulation where this ratio is O(1). This question is explored by varying the external mass transfer layer size 𝐻 in the PNM simulations.

The results are reported in Fig. 8. As can be seen, the impact of the initial breakthrough period becomes much less as the external mass transfer layer thickness is increased. Since 𝐻 / 𝛿 >> 1 in the experiments, the results displayed in Fig. 8 are consistent with the fact that the evaporation flux in stage 1 must be close to the potential evaporation in most experiments.

Another interesting feature shown in Fig. 8 lies in the fact that the evaporation rate is actually not really constant during stage 1 (after the initial BT period). In fact, the evaporation rate decreases smoothly during stage 1. As noted in [START_REF] Yiotis | Pore-network study of the characteristic periods in the drying of porous materials[END_REF], this is also observed in experiments, e.g.

[9], [START_REF] Coussot | Scaling approach of the convective drying of a porous medium[END_REF]. The variation of the evaporation rate in stage 1 also depends on the ratio 𝐻/d (or equivalently 𝐻/𝛿). The greater this ratio, the smaller is the variation of the evaporation rate during stage 1. As a result, the evaporation rate is practically constant over stage 1 only when 𝐻 / 𝛿 >> 1, i.e. when the thickness of the external mass transfer layer is sufficiently large compared to the grain size or the top edge effect region size.

VII. NLE EFFECT. MASS TRANSFER IN TOP EDGE EFFECT REGION.

In the LE models, the equilibrium is described via a relationship relating the equilibrium vapor partial pressure and the saturation. The equilibrium vapor partial pressure is less than the saturated vapor pressure because of a combination of capillary effects (Kelvin effect) and physisorption phenomena corresponding to the existence of thin liquid films over the pore wall. The equilibrium relationship is thus generally referred to as the desorption isotherm. However, in the PNM drying, physi-sorption phenomena and the Kelvin effect are not considered on the ground that their impact should be negligible in capillary porous media (an argument is that the drying of a circular or a square capillary tube can be quite well predicted [START_REF] Chauvet | Use of infrared thermography for the study of evaporation in a square capillary tube[END_REF] without considering adsorbed films and/or the Kelvin effect, the latter being non negligible in submicronic pores only). In other words, LE equilibrium in the case of the "PNM drying" means that ⟨𝑝 𝑣 ⟩ ≈ 𝑃 𝑣𝑠 if S > 0 where ⟨𝑝 𝑣 ⟩ is the averaged vapor partial pressure over the AV. Using as averaging volume (AV) the thin slices previously considered, the PNM simulations lead to the slice averaged vapor partial pressure profiles displayed in Fig. 9. The figure makes clear that: i) the NLE is noticeable, ii) but only in the top edge effect region, iii) the greater the external mass transfer layer thickness, the less is the NLE effect in the top edge effect region (as shown in the inset in Fig. 9). The results shown in the inset in Fig. 9 indicates that the NLE effect in the top edge effect region cannot be represented by a desorption isotherm like relationship since the deviation of the vapor partial pressure from the saturated vapor pressure depends on both the local saturation [START_REF] Moghaddam | Evaporation in capillary porous media at the perfect piston-like invasion limit: Evidence of nonlocal equilibrium effects[END_REF] and the external mass transfer.

An obvious consequence of the NLE effect in the network top region is that the average vapor partial pressure at the network top surface is less than the saturated vapor pressure during stage 1. However as indicated in Fig. 9, 𝑃 𝑣-𝑡ℎ𝑖𝑛-𝑠𝑙𝑖𝑐𝑒 (0) → Pvs for H >> δ. When the interfacial resistance becomes negligible compared to the external mass transfer resistance, then the approximation 𝑃 𝑣-𝑡ℎ𝑖𝑛-𝑠𝑙𝑖𝑐𝑒 (0) ≈ Pvs becomes acceptable. These concepts of interfacial and external resistances are discussed in more details in the next section.

VIII. MACROSCOPIC TRANSFER COUPLING AT THE POROUS MEDIUM

SURFACE

As mentioned in section II, the full formulation of the continuum model requires to specify boundary conditions at the evaporative surface of the medium, the top surface in our case. Three main different approaches have been considered in the literature. In this section, these approaches are recalled and discussed in order to select the one, if any, which seems the most appropriate to perform the comparison between PNM simulations and continuum modeling.

A. The three main approaches

The simpler approach is based on the concept of critical saturation Scr (or equivalently of critical water pressure) [START_REF] Li | Evaluation of model concepts to describe water transport in shallow subsurface soil and across the soil-air interface[END_REF]. The boundary condition is then expressed as

𝑗 = 𝑗 𝑝𝑜𝑡 when S(0) > Scr S(0) = Scr else (9)
Thus, the evaporation flux is equal to the potential evaporation jpot as long as the saturation computed at the surface is greater than Scr. As soon as S(0) = Scr, the boundary condition switches to the constant saturation condition S(0) = Scr. In the case of the PNM simulation, the saturation at the surface decreases to zero during stage 2 (Fig. 3). Thus, one should take Scr = 0.

However, the end of stage 1 is observed at a much greater saturation ~ 0.3 and 𝑗 ≈ 𝑗 𝑝𝑜𝑡 in stage 1 only when H >> δ. Although the critical saturation concept is commonly used, it is also oftentimes debated [START_REF] Li | Evaluation of model concepts to describe water transport in shallow subsurface soil and across the soil-air interface[END_REF]. In any case, Eq.( 9) is clearly inappropriate to model the PNM drying.

In many studies based on LE continuum models, e.g. [START_REF] Chen | A mathematical model of drying processes[END_REF][START_REF] Kaviany | Funicular state in drying of a porous slab[END_REF] among others, a mass transfer coefficient hm at the surface is introduced and the boundary condition is typically expressed as

𝑗 = ℎ 𝑚 𝑀 𝑣 𝑅𝑇 (𝑃 𝑣𝑖 -𝑃 𝑣∞ ), (10)
where 𝑃 𝑣𝑖 and 𝑃 𝑣∞ are the vapor partial pressure at the porous medium surface and in the external gas away from the porous surface, respectively. Expressing the mass transfer coefficient in terms of external flow parameters, such as the Reynolds or Schmidt numbers, is however generally not sufficient to reproduce the experimental data, e.g. [START_REF] Chen | A mathematical model of drying processes[END_REF]. The mass transfer coefficient is actually fitted assuming that it depends on the saturation at the surface. Actually, since the saturation is the main variable in the LE continuum model, the surface vapor partial 20 pressure 𝑃 𝑣𝑖 should be also expressed as a function of saturation. The commonly used approach is to assume that the LE relationhip corresponding to the desorption isotherm

This boundary is referred to as TBC2. As discussed in previous works [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 | A pore network study of evaporation from the surface of a drying non-hygroscopic porous medium[END_REF], this approach is clearly not applicable in the case of the comparison with the PNM simulations since, as discussed before, the fact that the vapor partial pressure can be lower than the saturated vapor pressure for S > 0 is the signature of a NLE effect and has nothing to do with adsorption phenomena or Kelvin effect. Furthermore, the network top surface is actually the place where the NLE effect is the strongest as illustrated in Fig. 9. A somewhat naïve idea could be to introduce a pseudo-desorption isotherm or more appropriately a NLE relationship 𝜑 𝑁𝐿𝐸 (𝑆) = 𝑃 𝑣 (𝑆) 𝑃 𝑣𝑠 but actually such a simple one to one relationship does not exist since

𝑃 𝑣 (𝑆) 𝑃 𝑣𝑠
can take different values for the same value of the local saturation depending on the overall saturation [START_REF] Moghaddam | Evaporation in capillary porous media at the perfect piston-like invasion limit: Evidence of nonlocal equilibrium effects[END_REF] or the external mass transfer layer thickness H as illustrated in Fig. 9.

An alternative boundary condition (TBC3) used mainly by soil scientists, e.g. [START_REF] Li | Evaluation of model concepts to describe water transport in shallow subsurface soil and across the soil-air interface[END_REF][START_REF] Smits | An evaluation of models of bare soil evaporation formulated with different land surface boundary conditions and assumptions[END_REF] and references therein, consists in expressing the boundary condition at the top surface as

𝑗 = 𝑀 𝑣 𝑅𝑇 (𝑃 𝑣𝑠 -𝑃 𝑣∞ )
𝑟 𝑎 +𝑟 𝑝𝑚 [START_REF] Ahmad | From micro-scale to macroscale modeling of solute transport in drying capillary porous media[END_REF] where 𝑟 𝑎 represents the external mass transfer resistance and 𝑟 𝑝𝑚 is the resistance due to the mass transfer in the porous medium. In this approach, it is assumed that the external mass transfer resistance 𝑟 𝑎 can be determined only from the consideration of the external transfer at the beginning of the drying process. For convenience, this resistance is denoted by 𝑟 𝑎0 since it is supposed to correspond to the external resistance at t = 0 when the porous medium is fully liquid saturated. In other words, this resistance is considered as being independent of the changes occurring in the fluid distribution at the surface during stage 1. By contrast, 𝑟 𝑝𝑚 does depend on the saturation [START_REF] Smits | An evaluation of models of bare soil evaporation formulated with different land surface boundary conditions and assumptions[END_REF].

B. Exploring the resistance approach

At first glance, the approach based on the resistance concept is attractive because of the existence of the top edge effect region in which the NLE effect is present. The mass transfer resistance 𝑟 𝑝𝑚 would be associated with the transfers within the top edge effect region, noting furthermore that the vapor partial pressure is indeed very close to 𝑃 𝑣𝑠 at the bottom of the top edge effect region (Fig. 9). Since the external transfer is purely diffusive in our PNM simulations, 𝑟 𝑎0 can be expressed as

𝑟 𝑎0 ≈ 𝐻 𝐷 𝑣
, where 𝐷 𝑣 is the vapor molecular diffusion in the binary mixture (vapor + air). From the PNM simulations, i.e. the computation of the evaporation flux j, 𝑟 𝑝𝑚 can therefore be computed as As can be seen, the interfacial resistance increases along the drying process during stage 1, i.e.

𝑟 𝑝𝑚 = 𝑀 𝑣 𝑅𝑇 (𝑃 𝑣𝑠 -𝑃 𝑣∞ ) 𝑗 -𝑟 𝑎0 (13
increases as the saturation in the medium decreases. This variation is qualitatively in agreement with the parametrization used by soil scientists [START_REF] Smits | An evaluation of models of bare soil evaporation formulated with different land surface boundary conditions and assumptions[END_REF][START_REF] Aa Van De Griend | Bare soil surface resistance to evaporation by vapor diffusion under semiarid conditions[END_REF], which reads,

𝑟 𝑝𝑚 = 10𝑒 0.3563𝜀(𝑆 𝑟 -𝑆 𝑡𝑜𝑝 ) (14)
where 𝑆 𝑟 is the residual saturation and 𝑆 𝑡𝑜𝑝 is the saturation in the top 1 cm layer of the soil.

The use of Eq.( 14) first requires to specify 𝑆 𝑟 and 𝑆 𝑡𝑜𝑝 . In the case of our PNM simulations, 𝑆 𝑟 can either be specified as 0.3 which corresponds to the saturation when the liquid phase ceases to be percolating (MCD) and becomes formed by isolated clusters (this would correspond to the standard definition of the residual saturation) or as 0 since adsorption phenomena are not considered (as exemplified in Fig. 3, the surface eventually fully dries in the PNM drying).

Both options have been tested. They lead to similar results. Thus, the results obtained considering 𝑆 𝑟 =0.3 are presented in what follows. Then, there is the question of 𝑆 𝑡𝑜𝑝 . One option could be to consider the average saturation in the top edge effect region, which could play a role similar to the soil top layer. However, the boundary condition expressed by Eq.( 12)

is to be used in conjunction with continuum models. Since the top edge region is quite thin, we argue that the target of continuum models for the considered capillary regime is to predict the saturation in the bulk and not the rapid saturation variation in the top edge effect region. For this reason, we took 𝑆 𝑡𝑜𝑝 = Sbulk. Then the direct application of Eq.( 14) led to quite poor results with values about two orders of magnitude too low compared to the values determined from the PNM simulations. However, if one adjusts the two coefficients in the equation, one obtains the results shown in Fig. 11 (red curve), which indicates that the parametrization of the interfacial resistance with an expression similar to Eq.( 14) is fair for the considered case, namely H =10, at least in the range of saturation [0.4-0.7].

𝑟 𝑝𝑚 𝑟 𝑎0 = 2.38𝑒 -0.43(𝑆 𝑟 -𝑆 𝑏𝑢𝑙𝑘 ) - (15)
with 𝑆 𝑟 = 0.3. However, it can be seen from Fig. 11 that a linear fit is actually as good over the considered range of saturation. 16)).

The dashed line marks the saturation below which the comparison between the PNM simulations and the continuum models is considered as meaningful (see text).

Nevertheless, both fits, i.e. Eq.( 15) and the linear fit, fail to correctly represent the greater increase of 𝑟 𝑝𝑚 as the saturation ~0.3 marking the end of stage 1 is approached. Interestingly, the overall shape of the 𝑟 𝑝𝑚 curve in Fig. 11 looks very much like a standard capillary pressure curve, e.g. [START_REF] Dullien | Porous media: fluid transport and pore structure[END_REF]. This type of curve is commonly represented using a mathematical relationship proposed by Van-Genuchten [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF]. Adapted here, this leads to express

𝑟 𝑝𝑚 𝑟 𝑎0
as

𝑟 𝑝𝑚 𝑟 𝑎0 = 𝑟 𝑝𝑚-𝑟𝑒𝑓 𝑟 𝑎0 [( 1-𝑆 𝑟 𝑆-𝑆 𝑟 ) 1 𝑚 ⁄ -1] 1 𝑛 ⁄ (16)
where 𝑚 = 1 -1 𝑛 ⁄ . As can be seen, Eq.( 16) depends on three parameters: n, 𝑆 𝑟 and 𝑟 𝑝𝑚-𝑟𝑒𝑓 .

As before, 𝑆 𝑟 is taken equal to 0.3 (saturation approximately corresponding to the end of stage 1) whereas 𝑟 𝑝𝑚-𝑟𝑒𝑓 and n have been adjusted. As can be seen in Fig. 11 (VG curve), Eq.( 16)

with n=18 and 𝑟 𝑝𝑚-𝑟𝑒𝑓 𝑟 𝑎0 ⁄ = 2.03 leads to a satisfactory representation of the 𝑟 𝑝𝑚 curve over the whole range of considered saturations, in particular as the end of stage 1 is approached.

In summary, TBC3 is by far the most satisfactory option, compared to TBC1 or TBC2, to express the boundary condition for the continuum modelling of the PNM drying. However, it is not easy to use since the use of TBC3 requires a careful characterization of its variation with both the saturation and the external mass transfer resistance, i.e. H in our model.

C Interfacial resistance and effective surface approach

Although TBC3 appears as the most satisfactory, or perhaps the least bad, approach to model evaporation during stage 1, it can be noted that the interpretation of Eq.( 12) is not as obvious as it may at first appear. Eq.( 12

in the top edge effect region zone (where 𝑆 𝑡𝑒 and 𝐷 𝑒𝑓𝑓-𝑡𝑒 are the average saturation and the effective diffusion coefficient of the vapor in the top edge effect region, respectively) and combining Eq.( 17) and ( 18) leads to Eq.( 12 . This approach to Eq.( 12) seems in line with the mention that𝑟 𝑝𝑚 is linked to the transfers within the soil 1 cm top layer in [START_REF] Smits | An evaluation of models of bare soil evaporation formulated with different land surface boundary conditions and assumptions[END_REF][START_REF] Aa Van De Griend | Bare soil surface resistance to evaporation by vapor diffusion under semiarid conditions[END_REF]. Then, it would make sense to apply this condition in our case at z = 𝛿 (bottom of the top edge effect region) rather than at z = 0. Since very often 𝛿 << Hpm this should actually have no serious significance for the modelling of the numerous situations where the constraint 𝛿 << Hpm is satisfied. However, in the case of the PNM simulations discussed in the current article, the situation is different since 𝛿 ≈ 6a is not negligible compared to the height of the network Hpm ≈ 30a. More importantly, there is an inconsistency between the derivation leading to Eq.( 12) from Eqs.( 17) and ( 18) and the PNM results. The derivation is actually based on the assumption that the evaporation flux j is constant over the top edge effect region. In other words, the evaporation flux at the bottom of the top edge effect region is assumed to be identical to the evaporation flux at the top of the top edge effect region. This is not in agreement with the PNM simulations which shows that the vapor diffusive flux actually vanishes over the top edge effect region, i.e. jv = 0 at z = 𝛿. This is illustrated in Fig. 12. Therefore, this observation suggests that TBC3 must be applied at z = 0 and not at z = 𝛿. Thus, the interpretation of rpm must be different than simply a vapor diffusion resistance within the top edge effect region. Actually, evaporation takes place within the top edge effect region, i.e. at menisci located within the top edge effect region, as well as at the upper limiting surface of the top edge effect region where a fraction of the pores is filled with the liquid. Thus, the modelling of the evaporation flux with Eq.( 12) should rather be interpreted as follows in relation with the concept of effective surface, e.g. [START_REF] Chandesris | Jump conditions and surface-excess quantities at a fluid/porous interface: a multi-scale approach[END_REF], [START_REF] Guo | Effective surface and boundary conditions for heterogeneous surfaces with mixed boundary conditions[END_REF] among others. Within this view of the continuum modelling approach, the layer formed by the top edge effect region is replaced by a fictitious "bulk" layer of same thickness, thus a zone where the internal evaporation is negligible and the saturation is equal to the bulk saturation and thus does not vary spatially over this fictitious layer, at least in the case of the considered capillary regime. The boundary condition TBC3, i.e. Eq.( 12), is applied on top of the domain formed by the fictitious layer and the remaining porous medium. Thus, in this approach the total height of the porous medium Hpm is not modified and is the same as in the experiments (or the pore network in our case).

Then, the resistance rpm is interpreted as an interfacial resistance at the porous medium top surface accounting for the complex transfer occurring actually within and at the top of the top edge effect region. This view of the interfacial resistance rpm is further discussed in the next section.

IX. CAPILLARY REGIME CONTINUUM MODEL SOLUTION

A. Assuming a constant external resistance over time

Focusing on stage 1, the boundary condition at the top is expressed using the concept of fictitious effective surface discussed in the previous section by Eq.( 12), namely,

𝜌 𝑙 𝐷 𝑙 (𝑆) 𝜕𝑆 𝜕𝑧 = 𝑗 = 𝑀 𝑣 𝑅𝑇 (𝑃 𝑣𝑠 -𝑃 𝑣∞ ) 𝑟 𝑎0 +𝑟 𝑝𝑚 = at z = 0 ( 18
)
As explained, in the previous section, the NLE effect is confined in the top edge effect region.

Elsewhere, 𝑃 𝑣 = 𝑃 𝑣𝑠 . The mass transfer in this region is taken into account via the interfacial resistance 𝑟 𝑝𝑚 . Consistently with the interfacial resistance approach, Eq.( 2) can therefore be discarded. The problem to be solved over stage 1 is therefore

𝜀𝜌 𝑙 𝜕𝑆 𝜕𝑡 = 𝜕 𝜕𝑧 (𝜌 𝑙 𝐷 𝑙 (𝑆) 𝜕𝑆 𝜕𝑧 ) (19)
together with Eqs.(4) (zero flux condition at the bottom) and [START_REF] Yiotis | Coupling between external and internal mass transfer during drying of a porous medium[END_REF]. The solution to this problem is discussed in Appendix B, where it is shown that this solution indeed corresponds to quasiflat saturation profiles when the capillary number, which is a small parameter in the problem, is sufficiently small. The saturation is then simply given by the equation,

𝜀𝜌 𝑙 𝐻 𝑝𝑚 𝜕𝑆 0 (𝑡) 𝜕𝑡 ≈ -𝑗 (20)
where 𝑆 0 (𝑡) represents the leading term in the development of the saturation as a function of the capillary number (see Appendix B). Furthermore, in Appendix C, we argue that 𝑆 0 (𝑡) can be replaced by Sbulk in order to compare the continuum model solution with the PNM simulations on a small network. Hence,

𝜀𝜌 𝑙 𝐻 𝑝𝑚 𝜕𝑆 𝑏𝑢𝑙𝑘 𝜕𝑡 ≈ -𝑗 (21)
Eq.( 21) is to be solved in conjunction with Eq.( 18 As can be seen in Fig. 13, the agreement between the PNM data (Sbulk-PNM) and the CM model (red and blue dashed lines) is reasonably good. The discrepancy is consistent with the fact that top and bottom edge region compensation effect (see Appendix C)is not perfect, i.e. Sbulk > Snet (Fig. A2 in Appendix C). The good agreement is not surprising when the interfacial resistance is directly estimated from the PNM data (option #1) because the evaporation flux is of course very well predicted in this case since the interfacial resistance has been determined from the PNM data using Eq.( 13), i.e. in fact adjusted so that Eq.( 13) exactly reproduces the evaporation flux. The evaporation flux variation as a function of time is depicted in the inset in Fig. 13. The result obtained with option #1 is actually indiscernible from the PNM computed evaporation flux. Somewhat surprisingly at first glance, using the proposed parametrization of rpm, i.e.

option #2, also leads to results as good as with option#1 regarding the prediction of Sbulk as shown in Fig. 13. However, the impact of the approximation is clearly visible on the evaporation flux (inset in Fig. 13), consistently with the fact that the parametrization deviates from the "exact" interfacial resistance when the end of stage 1 is approached in Fig. 11. Nevertheless, it can be considered that these results are quite encouraging

B. The external mass transfer resistance is time dependent

Based on these results, it is tempting to conclude that the interfacial resistance concept, as introduced with Eq.( 12), is fully adequate to perform the coupling at the surface.

However, there is a flaw in the interfacial resistance approach used so far. To be fully consistent, this approach should lead to a fair estimate of the mean vapor partial pressure at the surface.

From Eq.( 12), the mean vapor partial pressure at the porous medium could therefore be expected to be given by

𝑃 𝑣𝑠𝑢𝑟𝑓 𝑃 𝑣𝑠 = 1 - 𝑗𝑟 𝑝𝑚 𝑝 𝑣𝑠 𝑅𝑇 𝑀 𝑣 (22)
The results obtained using Eq.( 22) are compared in Fig. 14 with the spatially averaged vapor partial pressure at the surface computed from the PNM simulations. The latter is computed as

𝑃 𝑣𝑠𝑢𝑟𝑓 = ⟨𝑝 𝑣𝑖 ⟩ 𝑠𝑢𝑟𝑓 = ∑ 𝑑 𝑡𝑖 2 𝑝 𝑣𝑖 𝑖=1,𝑛 ∑ 𝑑 𝑡𝑖 2 𝑖=1,𝑛 (23)
of the interfacial throats connecting the network to the external mass transfer layer at the network top surface, 𝑝 𝑣𝑖 is the vapor partial pressure at the entrance of surface throat # i. As can be seen from Fig. 14, Eq.( 23) (which corresponds to the black curves in Fig. 14) significantly underpredicts the mean vapor partial pressure at the surface.

The discrepancy is attributed to the fact that the vapor partial pressure is not spatially uniform at the surface except at the very beginning of the drying process. In other words, 𝑝 𝑣𝑖 in Eq. [START_REF] Ochoa-Tapia | Heat transfer at the boundary between a porous medium and a homogeneous fluid[END_REF] varies from one interfacial throat to the other at the surface with 𝑝 𝑣𝑖 = 𝑃 𝑣𝑠 at the entrance of the interfacial throats filled by liquid and 𝑝 𝑣𝑖 < 𝑃 𝑣𝑠 at the entrance of the gaseous interfacial throats. This spatial variability has been documented in previous works, [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF][START_REF] Moghaddam | A pore network study of evaporation from the surface of a drying non-hygroscopic porous medium[END_REF] but only for a single value of the external mass transfer layer thickness and is illustrated again in Fig. 15, which shows the standard deviation of the distribution of 𝑝 𝑣𝑖 at the surface divided by the mean vapour partial pressure at the surface (Eq.( 37)), a quantity referred to as the coefficient of variation As can be seen from Fig. 15, 𝜎 𝑝 𝑣𝑖 first increases rapidly during the BT period (during which Sbulk ~1). Once the top edge effect region is formed, the variation of 𝜎 𝑝 𝑣𝑖 is less but it can be seen that 𝜎 𝑝 𝑣𝑖 𝑝 𝑣𝑖 > 𝑠𝑢𝑟𝑓 increases during stage 1 and can be as high as about 0.15 for the lower mass transfer external layer thickness considered in Fig. 15. This indicates quite significant variations of 𝑝 𝑣𝑖 over the surface.

The variations of 𝑝 𝑣𝑖 over the network top surface are further illustrated in Fig. 16. For clarity, the vapor partial pressure values corresponding to each pore at the surface are presented in ascending order. This figure clearly shows the changes in the vapor partial pressure distribution at the surface. The figure suggests that the gradual invasion of the surface pores is the main factor modifying the distribution since the shape of the curves for the values of the normalized vapor partial pressure less than 1 are similar. 676 pores out of 900 (75%) are invaded in this example during the BT period as indicated by the curve for Snet = 0.9 in Fig. 16. Then, only 115 additional pores (13% of the surface pores) are invaded during the remaining period of stage 1.

Although this represents a relatively small number of pores, the impact is noticeable since the vapour partial pressure is less in the invaded pores compared to the pores filled by the liquid. This impact is illustrated in Fig. 16b The variations in the vapor partial pressure distribution at the surface necessarily affects the external mass transfer because the structure of the vapor partial pressure field near a surface where the vapor partial pressure is not spatially uniform is different from the structure for a spatially uniform distribution of the vapor partial pressure at the surface. This means that is erroneous to consider the external mass transfer resistance 𝑟 𝑎 as constant in Eq.( 12) for a given H. Thus, we consider Eq.( 24) where the crucial difference compared to Eq.( 12) lies in the fact that both 𝑟 𝑎 and 𝑟 𝑝𝑚 are deemed to depend on Sbulk for a given H. This leads to the results reported in Fig. 17 for both 𝑟 𝑎 and 𝑟 𝑝𝑚 . As can be seen from Fig. 17a, both 𝑟 𝑎 and 𝑟 𝑝𝑚 increases during stage 1. It can be noted that both 𝑟 𝑎 and 𝑟 𝑝𝑚 vary with H. In Fig. 17b, the values of 𝑟 𝑎 and 𝑟 𝑝𝑚 for a given H were divided by the corresponding values for a reference H, here H = 50. The results in Fig. 17 suggest a possible variable separation with 𝑟 𝑎 (𝑆 𝑏𝑢𝑙𝑘 , 𝐻) = 𝑓 𝑟 𝑎 (𝐻)𝑟 𝑎-𝐻 𝑟𝑒𝑓 (𝑆 𝑏𝑢𝑙𝑘 ) and 𝑟 𝑝𝑚 (𝑆 𝑏𝑢𝑙𝑘 , 𝐻) = 𝑓 𝑟 𝑝𝑚 (𝐻)𝑟 𝑝𝑚-𝐻 𝑟𝑒𝑓 (𝑆 𝑏𝑢𝑙𝑘 ), at least over the range of saturation over which the CM can be considered as relevant, i.e. sufficiently away from the BT period.

The results depicted in Fig. 17 were for a single realization of the network. Averages over 50 realizations of the 30 × 30 × 30 network for H = 10 are shown in Fig. 18. The shape of the curves for both 𝑟 𝑝𝑚 and 𝑟 𝑎 suggests that a parametrization using a Van-Genuchten like relationship, i.e. Eq.( 16), can be used again. Both resistances increase linearly with a decreasing saturation over a relatively large range of saturations whereas the increase is greater and nonlinear toward the end of stage 1 and over the BT period at the very beginning of drying.

FIG. 18. Variation of the external mass transfer resistance and interfacial resistance (inset) with

Sbulk for H = 10 averaged over 50 realizations of the 30 × 30 × 30 network; 𝑟 𝑎0 = 𝐻 𝐷 𝑣 ⁄ is the external mass transfer resistance for the considered value of H at the very beginning of drying when the vapor partial pressure is spatially uniform at the surface.

In summary, the present study suggests that none of the coupling boundary conditions considered in previous works is fully satisfactory. However, the study indicates that expressing the coupling conditions via the concept of interfacial and external resistances can be a satisfactory approach provided that the impact of the evolution of the vapor partial pressure distribution at the surface during stage 1 is adequately taken into account. The most consistent approach implies to consider that both the interfacial and external resistances vary during stage 1 and not only the interfacial resistance as commonly considered. The fact that the simulations 34 indicate a possible variable separation in the dependence of the interfacial resistances with the saturation and H, i.e. the mass transfer rate, is an indication that the geometric effects (due to the configuration of the fluids) can be separated from rate effects that arise due to the nonequilibrium mass transport behavior. This opens up the route toward the in depth study of the dependence of the interfacial resistances with the fluid topology as briefly evoked in the next section.

X. DISCUSSION

Several issues still need clarification as regards the modelling of drying within the framework of the continuum approach to porous media. The present paper proposes a contribution from the consideration of the special case of the PNM drying in the capillary regime. The underlying idea was to test the continuum approach in the case of the drying process obtained via PNM simulations. In the case of the considered PNM drying, the NLE effect, the fact that the vapor partial pressure in less than the saturated vapor pressure in the presence of liquid, is a major feature, which should not be confused with adsorption or Kelvin effects since these effects are the evaporation flux at the surface can be separated into a contribution from the dry pores at the surface and a contribution from the wet pores, then the partition of the total evaporation flux into a dry pore contribution and a wet pore contribution should be parametrized [START_REF] Thiery | Drying regimes in homogeneous porous media from macro-to nanoscale[END_REF][START_REF] Lu | Insights into evaporation from the surface of capillary porous media gained by discrete pore network simulations[END_REF], noting that information on this partition is typically not available from the usual experimental data.

For the particular case considered in the present study, i.e. stage 1 in the capillary regime, this difficulty has been circumvented through the concept of interfacial resistance and the fact that the NLE effect is actually confined in a very thin region on top of the porous medium. In this approach, the two equation CM model reduces to a one equation model. 24) seems a promising approach as regards the long standing issue of the modelling of the coupling between the transfers in the porous medium and the external transfers in drying, many questions remain open and should be studied further. Often, stage 1 evaporation ends in practice not only because of the fragmentation of the liquid phases in many non-percolating clusters as in our PNM simulations but because of the impact of viscous effects in the liquid phase, thus because of a transition from the capillary regime to the viscous-capillary regime. In this context, it is unclear whether the parametrization of the interfacial resistance rpm and the external resistance ra determined from the consideration of the capillary regime still fully applies in the viscouscapillary regime. Also, the modelling of the full drying process is to be clarified, especially at the stage 1 -stage 2 transition. It is tempting to extend the concept of interfacial resistance to the transition zone between the dry zone and the wet zone during the receding front period (RFP) occurring in stage 2 but it can be argued that the size of the NLE zone in the vicinity of the receding front is perhaps not very thin anymore, which could question the applicability of the interfacial resistance concept to the other drying periods. Naturally, it is expected that the consideration of the continuous modelling of the PNM drying is also insightful for the drying process observed in the experiments. In this respect, it has been pointing out in previous works that the PNM simulations are consistent with many features observed experimentally.

Nevertheless, a recurrent question is the impact of liquid films, [START_REF] Yiotis | Coupling between external and internal mass transfer during drying of a porous medium[END_REF][START_REF] Prat | On the influence of pore shape, contact angle and film flows on drying of capillary porous media[END_REF][START_REF] Yiotis | Drying in porous media with gravitystabilized fronts: Experimental results[END_REF] that can form during the drying process. The liquid films were not considered in the PNM version used for the present study, which therefore is representative of the situations where the pore geometry and the contact angle are such that liquid films can be ignored. The liquid films can potentially greatly reduce the interfacial resistance if they are present within the top edge effect region since they help maintaining the vapor partial pressure closer to the saturation vapor pressure in this region compared to the situation in the absence of films [START_REF] Yiotis | Coupling between external and internal mass transfer during drying of a porous medium[END_REF][START_REF] Prat | On the influence of pore shape, contact angle and film flows on drying of capillary porous media[END_REF]. The detailed study of the liquid film impact on the interfacial and external resistances is thus an interesting open question.

One must also keep in mind that PNM introduce geometric simplifications compared to the complexity of real microstructures. It can be argued that the pore-structure represented by a PNM is less heterogeneous compared to most real microstructures. This should impact the fluidfluid surface area, which in turn will have impact on the mass transfer process in the interfacial region. A possibility to explore this issue could be to rely on some direct numerical simulations [START_REF] Panda | Influence of thermal gradients on the invasion patterns during drying of porous media: A lattice Boltzmann method[END_REF][START_REF] Qin | Lattice Boltzmann modeling of drying of porous media considering contact angle hysteresis[END_REF] on microstructure digital images using the results of the present article as guidelines.

However, it should be clear that this would represent a quite heavy task due to the computational resources needed to perform this type of simulations.

In our simulations, only cubic networks were considered and no particular effort was made to characterize the topology of the fluid fluid interface. However, recent works, see [START_REF] Mcclure | Modeling Geometric State for Fluids in Porous Media: Evolution of the Euler Characteristic[END_REF] and references therein, indicate that the fluid topology can be measured based on the Euler characteristic. Furthermore, it is also suggested, see again [START_REF] Mcclure | Modeling Geometric State for Fluids in Porous Media: Evolution of the Euler Characteristic[END_REF] and references therein, that a unique relationship exists between saturation, interfacial area, curvature and Euler characteristic. Then, a possible interesting approach could be to explore in depth the dependence of the interfacial resistances with the fluid topology in the interfacial region via the Euler characteristic for instance, along similar lines of the approach presented in [START_REF] Mcclure | Modeling Geometric State for Fluids in Porous Media: Evolution of the Euler Characteristic[END_REF]. In this respect, the fact that the resistances can be represented using Van-Genuchten relationships (Fig. 18) should not be a pure coincidence since the resistance variation is directly related to the changes in the fluid distribution in the interfacial region, and thus the capillary pressure. This is an additional indication in favor of future in depth studies on the dependence of the interfacial resistances with the fluid topology.

Nevertheless, the interfacial resistance approach must be considered as heuristic since Eq.( 24)

has not been derived mathematically but merely postulated. In this respect, some of the techniques used to model transfers between a porous medium and a free fluid in previous works, e.g. [START_REF] Beavers | Boundary conditions at a naturally permeable wall[END_REF][START_REF] Chandesris | Jump conditions and surface-excess quantities at a fluid/porous interface: a multi-scale approach[END_REF][START_REF] Valdés-Parada | Diffusive mass transfer between a microporous medium and an homogeneous fluid: Jump boundary conditions[END_REF][START_REF] Ochoa-Tapia | Heat transfer at the boundary between a porous medium and a homogeneous fluid[END_REF], could help setting Eq.( 24) on a more firm theoretical basis. However, compared to the interfacial problems in these references, the drying problem presents the additional difficulty of an evolving heterogeneity since the distribution of the menisci in the interfacial region, i.e. the fluid topology, varies with time. This makes perhaps a possible connection with previous works addressing the impact of other types of evolving heterogeneities, e.g. [START_REF] Cushman | Fickian dispersion is anomalous[END_REF][START_REF] Cushman | Nonlocal dispersion in media with continuously evolving scales of heterogeneity[END_REF].

As a result, one can wonder whether the interfacial coupling process studied in the present paper would not require fundamentally more sophisticated approaches than the simple parametrization represented by Eq. [START_REF] Ochoa-Tapia | Heat transfer at the boundary between a porous medium and a homogeneous fluid[END_REF]. Such an endeavor might be an interesting objective for future works. Nevertheless, it can be noted that the diffusion problem governing the mass transfer can be considered as quasi-steady in slow drying, which is an interesting simplification compared to the time dependent processes studied in [START_REF] Cushman | Fickian dispersion is anomalous[END_REF][START_REF] Cushman | Nonlocal dispersion in media with continuously evolving scales of heterogeneity[END_REF].

XI. CONCLUSION

The often debated problem of the modelling of the drying process during stage 1 evaporation drying period within the framework of the continuum approach to porous media was studied from comparisons with PNM simulations focusing on the capillary regime in the classical situation of a porous medium column open at the top.

The study first shows that the drying process results in the formation of two edge effect regions, at the top and at the bottom. The top edge effect region is of special interest since this region is where the liquid-vapor phase change, i.e. the evaporation, takes place. It has been shown that this region is a reminiscence of the breakthrough (BT) period occurring at the very beginning of the drying process. It has been argued that this initial period is very difficult to detect in the experiments because it typically corresponds to a very small mass loss compared to the liquid mass initially present in the porous medium. The size of the edge effect regions has been characterized and has been found to be on the order of a few lattice spacing only, independently of the network size. Also, it has been shown that the non-local equilibrium (NLE) effect is confined within the top edge effect region.

Since this region is quite thin, it has been argued that the concept of interfacial resistance in which the top edge effect region is represented by a surface of zero thickness (a representation of the interfacial zone referred to as a Gibbs representation) was adapted for describing the coupling between the external mass transfer and the transfer in the porous medium at the top of the porous medium. The interfacial resistance was computed from PNM simulations and shown to increase along stage 1. For a given external mass transfer resistance, the interfacial resistance can be parametrized over a significant range of saturations using an expression similar to the one used in previous works by soil scientists. However, the numerical factors are different and the latter expression actually reduces to a linear variation. Also, this expression is not adapted to parametrize the interfacial resistance over the initial breakthrough period and when the end of stage 1 is approached. A very classical relationship used to represent the capillary pressure curve can be adapted to represent the interfacial resistance over the whole stage 1. More importantly, the interfacial resistance was found to depend not only on the saturation, i.e. the fluid topology, but also on the external mass transfer resistance, i.e. the mass transfer rate. Also, the interfacial resistance was parametrized using the "bulk" saturation and not the saturation in the top edge effect region on the ground that the continuum model cannot predict the saturation variation in this very thin region. Most importantly, it has been shown that both the interfacial resistance and the external mass transfer resistance must be parametrized as functions of both the bulk saturation and the external mass transfer layer thickness. However, our results indicate a possible variable separation. In other words, the geometric effects due to the configuration of the fluids can be separated from rate effects that arise due to the non-equilibrium mass transport.

In summary, the main outcome of the present article is to clarify the boundary condition to be specified at the (evaporative) top surface via the concept of effective interfacial resistance.

In addition, the present paper clarifies how to handle the comparison between PNM simulations and continuum models so as to limit the impact on the comparison of the finite size effects inherent in the PNM simulations. In brief, the small size of the networks typically considered in the PNM simulations lead to a BT period much longer than in the experiments. This period cannot be captured by the commonly used continuum models. As a result, the comparison should start a while after the end of the BT period so that the length scale separation concept underlying the continuum approach is verified.

As mentioned before, all this has been developed considering only stage 1 and the capillary regime in the absence of noticeable temperature variations. Much work is still needed to explore whether the concept of interfacial resistance is still pertinent beyond the situation considered in the present paper and can be extended to other drying regimes as well as to the modelling of the full drying process.

Furthermore, it must be noticed that the formulation of the interfacial coupling via the resistance concept was merely postulated based on similar formulations in previous works. Therefore, there is need for future works aiming at setting the coupling in a rigorous framework. Also, Where 𝑢 𝑒𝑣. is the evaporation velocity (𝑢 𝑒𝑣. = 𝑗 𝜌 𝑙 ⁄ . Eq.(A6) makes clear that 𝜂 is actually a capillary number also taking into account the ratio between the macroscopic length 𝐻 𝑚𝑝 , which is the relevant length for estimating the viscous pressure drop and the pore size r (it is recalled that 𝑘 ∝ 𝑟 2 , 𝑟 𝑘 ∝ 𝑟 -1 ) , which is the relevant length for characterizing the capillary pressure.

capillary numbers Ca, that is in the so-called capillary regime. In summary the saturation profile is given by 𝑆 = where 𝑟 𝑝𝑚 is a function of 𝑆 0 (𝑡). This solution can thus be considered as the asymptotic solution to the continuum model in the limit of sufficiently small capillary numbers (note that Ca varies during stage 1 with the definition of Ca given by Eq.(A6)).

C Bottom and top edge regions compensation effect

Based on the elements presented in the main text, it is first considered that the target of the continuum model must be the period after the BT period for two main reasons: i) the BT period corresponds to the formation of a fractal object, a percolation cluster, which is incompatible with the length scale separation concept underlying the continuum approach to porous media, ii) the BT period is actually indiscernible in most practical situations of interest. Second, since the top edge effect region is expected to be extremely thin compared to the typical total height of porous samples in practice and typically thinner than the size of a REV, the objective of a continuum model should not to capture the rapid saturation variation in the top edge effect region, nor actually in the bottom edge effect region which is also very thin. In other words, the main objective of a continuum model is to predict the (quasi-)uniform distribution of the saturation over the height of the porous sample observed after the very short initial BT period and over most of the sample height in conjunction with the evolution of the evaporation rate.

The latter would be obtained from Eq.( 24) via the adequate parametrization of the interfacial and external resistances. In what follows, we argue that the above continuum model can be tested from comparisons with the PNM simulation results despite the small size of the considered network.

To this end, the NLE continuum model, i.e. Eqs.

FIG. 1 .

 1 FIG.1. Sketch of the archetypical drying situation considered in the study.

  and references therein, the NLE continuum model reads (for the 1D situation considered in the current paper),

  The saturation is supposed to be spatially uniform initially S = S0. The zero flux boundary conditions at the bottom read -𝜌 𝑙 𝐷 𝑙(𝑆)

Fig. 4 where

 4 Fig.4where Snet is the overall saturation in the network.

FIG. 5 .

 5 FIG.5. IP thin slice averaged saturation profiles for Snet = 0.3, 0.4, 05, 0.6, 0.7, 0.8, 0.9 and N = 20 (averaged over 100 realizations) and N = 100 (averaged over 20 realizations). The overall saturation at breakthrough is Snet ≈ 0.83 for N = 20 and Snet ≈ 0.92 for N =100.

FIG. 6 .

 6 FIG. 6. IP saturation profiles at breakthrough for various network sizes. The inset shows the saturation at network top in the first slice (z = 0) as a function of network size. Profiles are averaged profiles over 150 (N = 20), 100 (N=30), 80 (N=40), 60 (N=50), 30 (N=100), 10 (N=200) realizations, respectively

FIG. 7 .

 7 FIG. 7. Top and bottom edge effect region size. IP thin slice averaged saturation profiles for Snet = 0.5 for various network sizes. The inset in the figure on the left shows a zoom over the top region. The figure on the right shows the overlap of the profiles for the various network sizes in the bottom region. The inset in the figure on the right shows a zoom over the bottom region. Profiles are averaged profiles over 100 (N=20), 80 (N=40), 60 (N=50), 30 (N=80), 20 (N=100) realizations, respectively

FIG. 8 .

 8 FIG.8. PNM simulation drying curve over stage 1 for one realization of the 30 × 30 × 30 network for various external mass transfer layer thickness H (measured in lattice spacing unit a). Jref is the evaporation rate at the very beginning of drying for H = 5. Snet is the overall saturation in the network.

FIG. 9 .

 9 FIG. 9. Thin slice averaged vapor partial pressure profiles for one realization of the 30 × 30 × 30 network for Snet = 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 and H = 5. The inset shows the variation of the thin slice averaged vapor partial pressure on top of network as a function of H for Snet = 0.4.

  𝜑(𝑆) -𝑃 𝑣∞ ).

)FIG. 10 .

 10 FIG. 10. Interfacial resistance from Eq.(13) as a function of Sbulk for various external mass transfer layer thicknesses H. Results are for one realization of the 30 × 30 × 30 network.

FIG. 11 .

 11 FIG. 11. Interfacial resistance as a function of Sbulk. The thick black curve is obtained from Eq.(13) and corresponds to the average over 50 realizations of the 30 × 30 × 30 network. The thin black line corresponds a linear fit over the range of saturation [0.4 -0.7]. The VG blue curve is obtained adjusting the parameters of a Van-Genuchten like expression (see text and Eq.(16)).

  ) looks like the combination of the steady-state diffusive transfer in the top edge effect region and the external mass transfer layer. Expressing the evaporation flux as mass transfer layer and as 𝑗 = 𝜀𝑆 𝑡𝑒 𝐷 𝑒𝑓𝑓-𝑡𝑒 (𝑆 𝑡𝑒 )

FIG. 12 .

 12 FIG. 12. Distribution of vapor diffusive flux over the top edge effect region from PNM simulations. Results are for one given realization of the 30 × 30 × 30 network.

= 2 .

 2 FIG.13. Variation of 𝑆 𝑏𝑢𝑙𝑘 as a function of time during stage 1. Comparison between the PNM data (averages over 50 realizations of the 30 ×30× 30 network) and the results from the continuum model. The comparison is performed for overall saturations lower than 0.7 (part of the figure on the right side of the vertical dashed line) for the comparison not to be hampered by the initial significant finite size effect associated with the BT period (see text). The reference time is the time at the end of stage 1. The inset shows the variation of the evaporation flux for the two considered options (see text).

FIG. 14 .

 14 FIG. 14. Variation of the mean vapor partial pressure at the porous surface 𝑃 𝑣𝑠𝑢𝑟𝑓 as a function of Sbulk for various external mass transfer layer thicknesses H with comparison with PNM results. The black curves (CM) are obtained using Eq.(23). Results are for one given realization of a 30 × 30 × 30 network.

FIG. 15 .

 15 FIG.15. Standard deviation of the vapor partial pressure distribution at the surface as a function of Sbulk for various mass transfer external layer thicknesses. Results are for one given realization of the 30 × 30 × 30 network.

  FIG. 16. a) Vapor partial pressure variations over the network top surface 900 pores for various overall network saturations. The vapor partial pressures are sorted in ascending order in the main figure. The inset shows the vapor partial pressure variations for two overall network saturations when the surface pores are scanned row by row; b) Evaporation rates from the 900 network top surface pores for various overall network saturations. The evaporation rates are sorted in ascending order in the main figure. The inset shows the evaporation rates for two overall network saturations when the surface pores are scanned row by row; 𝐽 𝑝-𝑟𝑒𝑓 = 𝐽 0 𝑁×𝑁 where𝐽 0 is the evaporation rate when the network is fully liquid saturated. Results are for one given realization of a 30 × 30 × 30 network.

  𝑎 is estimated using the PNM data from the equation, . a)Variation of the external mass transfer resistance and interfacial resistance (inset) with Sbulk for various external mass transfer layer thicknesses H; 𝑟 𝑎0 = 𝐻 𝐷 𝑣 ⁄ is the external mass transfer resistance for the considered value of H at the very beginning of drying when the vapor partial pressure is spatially uniform at the surface; b) External mass transfer resistance rescaled by its value for H = 50 as a function of Sbulk. The inset shows the interfacial resistance rescaled by its value for H = 50 as a function of Sbulk The vertical dashed lines correspond to Sbulk = 0.7 considered as the saturation range upper bound over which the comparison with the CM should be performed. Results are for one given realization of a 30 × 30 × 30 network.

  not taken into account in the considered drying PNM. In principle, the consideration of the NLE effects implies to consider a two equation continuum model (CM) as presented in Section II. Then the question arises which boundary conditions must be imposed at the porous medium open surface. The boundary condition at the top is indeed presumably much trickier than for the commonly used one equation LE model since in principle a boundary condition for each equation, i.e. Eqs. (1) and (2), must be imposed at the medium top surface. If one considers that

  The important difference compared to the commonly used one equation LE CM model is that the remaining equation takes into account the transport in liquid phase only whereas in the usual LE CM model the liquid mass balance equation and the vapor mass balance equation are combined to obtain a single equation. Although the concept of interfacial resistance via Eq.(

1 𝐷 2 𝜂( 1 - 2 𝜂

 1212 FIG. A1. Schemaric of network and external mass transfer layer.

( 1 )

 1 FIG. A2. Mean saturation𝑆 > 𝑡𝑧 in the top edge effect region, mean saturation in the bottom edge effect region 𝑆 > 𝑏𝑧 and overall network saturation Snet as a function of the saturation in the bulk. These results are averages considering 200 realizations of the 30 × 30 × 30 network (𝑆 > 𝑏𝑧 from IP simulations) or 50 realizations (𝑆 > 𝑡𝑧 using the drying PNM algorithm) respectively.

Acknowledgments

Financial support from joint project''Drycap'' funded by GIP ANR (project16-CE92-0030-01) and DFG (project TS28/10-1) is gratefully acknowledged. Part of the numerical simulations has been performed with the computing resources of CALcul en Midi Pyrénées (CALMIP, Toulouse, France). We are also grateful to an anonymous referee for suggesting several interesting ideas to analyze and extend the present work.

APPENDIX A Computation of evaporation flux in the PNM simulations

A schematic of the computational domain in the PNM approach is shown in Fig. A1. A 2D mesh is considered for simplicity but the actual mesh is 3D. Also, only 2 vertical rows of pores are shown for simplicity. The mesh cells in the external mass transfer layer (external diffusive layer in Fig. A1) are shown as well as the computational nodes. For a given fluid configuration, the steady state diffusion equation ∆𝑝 𝑣 = 0, where∆ is the Laplacian operator, is solved using a standard finite volume approach, e.g. [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF]. This gives the vapor partial pressure in each node of the computational domain (i.e. all nodes in the external diffusive layer and the nodes occupied by the gas phase in the network). Then, the evaporation rate from the network can be computed at the interface between the network and the diffusive layer or, more simply, at the top of the diffusive layer (under the considered quasi-steady approximation, the mass transfer rate is the same through each horizontal plane in the external diffusive layer),

Where 𝑗 𝑣 corresponds to the vapor mass flux from the liquid free pores at the porous medium top surface and 𝛿 𝑡 is the size of the top edge effect region. Integrating Eq. ( 1) leads to

where 𝑗 𝑙 (0) is the evaporation flux from the menisci located at the entrance of the liquid filled pores at the porous medium top surface and 𝑗 𝑙 (𝛿 𝑡 ) is the liquid flow rate per unit surface area in the liquid phase at the bottom of the top edge effect region. Combining Eq.( A8) and (A9)

and noting that 𝑗 𝑙 (0) + 𝑗 𝑣 = j yield

Since evaporation is negligible below the top edge effect region (Fig. 12), integrating Eq.(

over the "bulk" region and the bottom edge effect region leads to