Faeez Ahmad

Arman Rahimi

Evangelos Tsotsas

Marc Prat

Abdolreza • Kharaghani

From micro-scale to macro-scale modeling of solute transport in drying capillary porous media

Keywords:

Introduction

The drying of capillary porous media saturated with saline so lution is an important process as it has several environmental and industrial applications. The convective transport of dissolved solute towards the evaporating surface driven by liquid capillary pumping induces a solute concentration gradient inside the porous medium which is simultaneously counteracted by a diffusive back flow [START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF][START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF][START_REF] Guglielmini | Drying of salt solutions in porous materials: Intermediate-time dynamics and efflorescence[END_REF]. Once the solute concentration sufficiently exceeds the solubility limit. the sait precipitates and forms solid crystals. The fonnation of crystals can have adverse effects such as mechanical damage to the porous materials, changes in the drying kinetics and transport properties of the porous medium such as pore size distribution, permeability and vapor pressure of the dried fluid [START_REF] Scherer | Stress from crystallization of salt[END_REF][START_REF] Naillon | Ion transport and precipitation kinetics as key aspects of stress generation on pore walls induced by salt crystallization[END_REF][START_REF] Larsen | Groundwater salinity influenced by holocene seawater trapped in incised valleys in the red river delta plain[END_REF][START_REF] Börnhorst | Influence of pore structure and impregnation-drying conditions on the solid distribution in porous support materials[END_REF]. The un derstanding of the intricate transport phenomenon of solute inside drying porous media is also required to address real world prob lems like sait weathering in building materials [START_REF] Goudie | Salt weathering hazards[END_REF] and soit salin ization [START_REF] Chen | Evaporation from a salt-encrusted sediment surface: Field and laboratory studies[END_REF] that affects crop production. Consequently, the mode!-

• Corresponding authors.

E-mail addresses: mprat@imft.fr (M. Prat), abdolreza.kharaghani@ovgu.de (A.

Kharaghani).

ing of the solute transport during drying of saline porous media has recently gained a lot of attention.

In this context. it is vital to understand the transport phe nomenon that leads to the crystallization of solute. Researches have already been conducted to study the complex processes that play a rote in the regime preceding the occurrence of crystal lization, e.g., see [START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF][START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF][START_REF] Guglielmini | Drying of salt solutions in porous materials: Intermediate-time dynamics and efflorescence[END_REF][START_REF] Hidri | Porous medium coffee ring effect and other factors affecting the first crystallisation time of sodium chloride at the surface of a drying porous medium[END_REF]. These studies are based essentially on continuum modeling (CM) (11 J and consist of solution of the advective-diffusive transport of solute. The commonly used macro scopie advective-diffusive transport equation (ADE) led to a quite good agreement with the measurements for the so-called evapora tion wicking situation [START_REF] Pel | A simplified model for the combined wicking and evaporation of a NaCl solution in limestone[END_REF]. As discussed in several previous works, e.g. [START_REF] Pel | A simplified model for the combined wicking and evaporation of a NaCl solution in limestone[END_REF][START_REF] Diouf | Locus of first crystals on the evaporative surface of a vertically textured porous medium[END_REF][START_REF] Veran-Tissoires | Discrete salt crystallization at the surface of a porous medium[END_REF][START_REF] Veran-Tissoires | Evaporation of a sodium chloride solution from a saturated porous medium with efflorescence formation[END_REF][START_REF] Hidri | Stagnation points as loci of solute concentration extrema at the evaporative surface of a random porous medium[END_REF][START_REF] Lazhar | Combined wicking and evaporation of NaCl solution with efflorescence formation: The efflorescence exclusion zone[END_REF], the later corresponds to a situation where the porous sample is fully saturated by the solution as the result of capillary rise. Afthough experimental results in qualitative agreement with the predictions of the ADE have been reported [START_REF] Pel | Ion transport and crystallization in inorganic building materials as studied by nuclear magnetic resonance[END_REF], it has been ar gued from other comparisons with experimental data, e.g. [19.20), that the classical macroscopic ADE might not lead to good results in the drying situation. The later refers to the commonly encoun tered situation where the solute transport occurs together with a decrease in liquid saturation of the porous medium resulting from evaporation.

For this reason, the focus of this study is on the drying situation. The objective is to assess and evaluate the commonly used model based on the classical macroscopic ADE [START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF][START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF][START_REF] Guglielmini | Drying of salt solutions in porous materials: Intermediate-time dynamics and efflorescence[END_REF]. Instead of laboratory experiments, we use pore network model (PNM) simulations as a benchmark. PNM simulations provide pore-scale insights analogous to that obtained in laboratory experiments using advanced imaging techniques. The effectiveness of PNM simulations to capture the behavior of drying porous media with various complex physical phenomena is well established in the literature e.g. [START_REF] Laurindo | Numerical and experimental network study of evaporation in capillary porous media. phase distributions[END_REF][START_REF] Prat | Recent advances in pore-scale models for drying of porous media[END_REF][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]. Recently, comprehensive pore network models for the drying induced solute transport inside porous media have been presented in [START_REF] Dashtian | Pore-network model of evaporation-induced salt precipitation in porous media: The effect of correlations and heterogeneity[END_REF][START_REF] Rahimi | Discrete modeling of drying induced ion transport and crystallization in porous media[END_REF]. In this manuscript, we employ a 3D PNM which is based on [START_REF] Rahimi | Discrete modeling of drying induced ion transport and crystallization in porous media[END_REF].

The PNM simulations make clear that the liquid distribution during drying is characterized by the formation of many liquid clusters [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF]. The commonly used macroscopic model does not explicitly take into account the liquid phase fragmentation process. The question thus arises as to whether the liquid phase fragmentation phenomenon has an impact on the solute dynamics during drying that cannot be captured by the traditional macroscopic approach. A related issue is whether the liquid clusters are or not interconnected by liquid films [START_REF] Prat | On the influence of pore shape, contact angle and film flows on drying of capillary porous media[END_REF][START_REF] Chauvet | Three periods of drying of a single square capillary tube[END_REF]. Since the liquid films can provide a pathway for the solute transport between clusters -this is a potentially important issue. However, since the contact angle of a saline aqueous solution is relatively high, the development of liquid films is not favored [START_REF] Sghaier | On the influence of sodium chloride concentration on equilibrium contact angle[END_REF]. For this reason, liquid films are assumed to be negligible in this study.

To explore the above questions, we focus on the capillarydominated regime in a homogeneous porous medium. This is a regime frequently observed in laboratory experiments, e.g. [START_REF] Pel | Ion transport and crystallization in inorganic building materials as studied by nuclear magnetic resonance[END_REF][START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF][START_REF] Coussot | Scaling approach of the convective drying of a porous medium[END_REF][START_REF] Eloukabi | Experimental study of the effect of sodium chloride on drying of porous media: The crusty-patchy efflorescence transition[END_REF][START_REF] Faiyas | Transport of a water-soluble polymer during drying of a model porous media[END_REF]. This is also the regime considered in previous works on the solute distribution during drying [START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF][START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF][START_REF] Guglielmini | Drying of salt solutions in porous materials: Intermediate-time dynamics and efflorescence[END_REF]. As reported for instance in [START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF], this regime is characterized by a first drying period, also referred to as the constant rate period (CRP), in which the evaporation rate is controlled by transport of vapor in the external boundary layer and is almost constant. In this regime, the saturation profiles are flat, that is spatially uniform. In this study, we focus on the CRP, which is also a period in which the liquid phase in the porous medium is connected to the porous medium evaporative surface.

As reported in several previous works, e.g. [START_REF] Lazhar | Combined wicking and evaporation of NaCl solution with efflorescence formation: The efflorescence exclusion zone[END_REF][START_REF] Shokri | Pore-scale dynamics of salt transport and distribution in drying porous media[END_REF][START_REF] Shokri-Kuehni | Iodine k-edge dual energy imaging reveals the influence of particle size distribution on solute transport in drying porous media[END_REF][START_REF] Eloukabi | Experimental study of the effect of sodium chloride on drying of porous media: The crusty-patchy efflorescence transition[END_REF], salt crystallization is frequently observed at the porous medium surface during the CRP. This is fully consistent with the semi-analytical or numerical solutions of the commonly used ADE; [START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF][START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF][START_REF] Guglielmini | Drying of salt solutions in porous materials: Intermediate-time dynamics and efflorescence[END_REF], which predict that the locus of the maximum solute concentration is at the drying surface. Since crystallization occurs when the solution concentration reaches a critical value, the most likely place of crystallization is the surface. As reported in [START_REF] Eloukabi | Experimental study of the effect of sodium chloride on drying of porous media: The crusty-patchy efflorescence transition[END_REF], the crystallization process at the surface and the subsequent development of a salt crust on top of the porous medium (a classical situation where the evaporative surface is on top of the porous sample is assumed) can have a severe impact on the evaporation rate. In this study, the focus is on the situation prior to the crystallization, i.e. on the solute distribution before the onset of crystallization.

The paper is organized as follows: In Sections 2 and 3, we briefly describe the classical CM formulation and PNM algorithm, respectively. Results of PNM simulations are presented in Section 4. A comparison between PNM results and CM solution is presented in Section 5. Sections 6 and 7 are aimed at deepening the understanding of the transport at the evaporative surface using PNM simulations. This is followed by Section 8 which consists of a discussion linking the various aspects of the study. Section 9 consists of the main conclusions of the study and recommendations for future work.

Continuum model

As exemplified in [START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF][START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF][START_REF] Guglielmini | Drying of salt solutions in porous materials: Intermediate-time dynamics and efflorescence[END_REF] , the continuum drying model in the presence of solute aims at predicting the evaporation rate, the liquid distribution, i.e. the saturation profiles, and the solute concentration distribution within the sample. While the corresponding model generally requires a numerical solution [ 10 , 13 ] for both the evolution of the saturation and the solute concentration, a much simpler approach is possible when only the CRP in the capillary regime is considered in the classical configuration where a homogeneous porous medium is limited by impervious walls at the bottom surface and along the lateral sides with evaporation only at the top surface. First, the evaporation flux j is constant in the CRP and therefore just an input. Second, the fact that the saturation, S , is spatially uniform can be exploited [START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF] to obtain its evolution as,

S = 1 - j ρ l Lε , ( 1
)
where ρ l is the liquid mass density, ε the porosity and L is the height of the porous sample. The advective-diffusive equation classically used for modeling the solute transport during the drying process in the considered 1D macroscopic configuration reads [START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF] ∂ ρ l εSC ∂t

+ ∂ ∂z ( ρ l εUSC ) = ∂ ∂z ρ l εSD * s ∂C ∂z , ( 2
)
where C is the solute concentration, z is the vertical spatial coordinate, U is the average liquid velocity, also referred to as the interstitial velocity, and D * s is the effective solute diffusion coefficient. The boundary condition needed to solve Eq. ( 2) is expressed by zero flux at both boundaries (top evaporative surface and the boundary on the opposite end of the surface) of the 1D domain expressed as

ρ l ε USC -ρ l ε SD * s ∂C ∂z = 0 . ( 3
)
The effective solute diffusivity D * s is a function of liquid saturation S . The respective function is expressed using the classical formulation [ 34 , 35 ], which depends on the pore-scale solute diffusivity, porosity and saturation percolation threshold S c ,

D * s = S -S c 1 -S c D * s ( 1 ) , (4)
where D * s [START_REF] Sghaier | On ions transport during drying in a porous medium[END_REF] is the value of effective solute diffusivity for fully saturated region, i.e., when S = 1. It is generally expressed as a function of porosity and tortuosity of the porous medium and porescale solute diffusion coefficient D s . Since D * s (1) depends on the microstructure of the considered porous medium, it needs to be determined for our particular porous medium, i.e. our pore network model. This is performed in Section 5.1 . The value of D s is considered as 1 × 10 -9 m 2 /s which is close to that of the common salts such as sodium chloride.

Like the saturation, the velocity field U(z) in Eq. ( 2) can also be obtained analytically [START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF] as

U = j ε ρ l S 1 - z L , ( 5
)
where z = 0 is at the porous domain bottom and z = L at the evaporative surface. Eq. ( 2) is solved numerically with MATLAB using the analytical solutions according to Eqs. (1) and ( 5) for the evolution of saturation and velocity.

Pore network modeling

As illustrated in Fig. 1 , the pore network model consists of a three-dimensional regular network of cylindrical throats which account for the void volume. These throats are connected to each other through nodes referred to as pores. The pores serve as computational points only and do not have a volume. The length of the cylindrical throats in the network is uniform, whereas their radii vary according to a normal distribution based on a mean radius and standard deviation. The network is isolated from the sides and bottom, whereas it is exposed to an external boundary layer at the top. Inside the network, the mass transport in the liquid phase is dominated by capillary forces, meaning that the influence of effects caused by the liquid viscosity are negligible and liquid can be pumped between interconnected interfacial liquid throats in the network regardless of the distance between the considered throats. The capillary pressure in the interfacial throats is computed from the Young-Laplace equation as

P c = 2 σ cos θ r ( 6
)
where P c , σ , θ and r denote the capillary pressure, liquid surface tension, equilibrium contact angle and throat radius, respectively. As explained in previous works, e.g. [START_REF] Laurindo | Numerical and experimental network study of evaporation in capillary porous media. phase distributions[END_REF][START_REF] Prat | Recent advances in pore-scale models for drying of porous media[END_REF][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] , the Young -Laplace equation dictates the throat to be invaded at the periphery of each liquid cluster (see below). According to Eq. ( 6) , this is the throat of greater radius at the cluster boundary. The vapor transport is based on Stefan's flow expressed by

J v ,k = π r 2 t ˜ M v ˜ R T L t P atm D v a ln P v ,i -P atm P v , j -P atm , ( 7
)
where J v ,k is the vapor mass flow rate between pores i and j which are connected through throat k, D va is the pore-scale vapor diffusion coefficient and P v indicates the partial vapor pressure. For the computation of vapor pressure field, a system of equations which is obtained by applying mass balance at each pore is solved numerically. The network is initially fully saturated with aqueous solution in which the concentration of dissolved solute is uniform. As drying goes on, because of capillary pumping the liquid phase gets distributed into progressively diminishing clusters. These clusters are accounted for by labelling them through a variant of the Hoshen-Kopelman algorithm [START_REF] Metzger | Erratum: Extension of hoshen-kopelman algorithm to non-lattice environments[END_REF] . The drying process is stopped when the first drying period ends, i.e. when the evaporative surface essentially becomes unsaturated. Since this is not the main focus of this study, we neglect the influence of dissolved solute on the transport properties of water. Detailed description on the PNM drying algorithm can be found in [ 22 , 23 ].

The advective-diffusive transport of the dissolved solute in the liquid phase is modeled at the pore scale within the throats. In the PNM algorithm that does not account for solute transport the drying process is discretized according to time for emptying of one meniscus throat, here referred to as global time step. During each discretized time step, the boundary conditions for vapor pressure field are assumed to be stationary. However, in the PNM that accounts for solute transport, the time step is discretized according to the stability criterion of the advection-diffusion equation with an explicit solution scheme [START_REF] Licsandru | Dissolution-precipitation-driven upward migration of a salt crust[END_REF] . Solute transport calculations are carried out within these time increments while assuming constant capillary flow rates until the global time step has elapsed. An adaptive sub-discretization algorithm, based on the local Péclet number, is utilized for the throats to improve numerical accuracy of the model while keeping the computational time in check. For more details on the PNM algorithm for modeling of solute transport, see [START_REF] Rahimi | Discrete modeling of drying induced ion transport and crystallization in porous media[END_REF] .

When the solute concentration in a throat reaches a threshold referred to as the critical concentration ( C cr ), crystallization occurs and the solution in the vicinity of the crystals rapidly becomes close to the solubility, also referred to as the saturation concentration ( C sat ), e.g. [ 5 , 38 ]. As a result, the solute concentration value does not supersede the threshold value of C sat after the time step has been elapsed. As mentioned earlier, since the focus in this study is on the solute distribution prior to crystallization, no particular threshold value for crystallization is considered. The concentration can thus increase as dictated only by the competition between advection and back-diffusion and the changes in liquid saturation.

Also, the influence of the solute concentration on the properties of the fluid is neglected in this study for simplicity. In this respect, the important feature is to impose the same fluid properties for the PNM simulations and the CM. Based on this, the value of initial solute concentration can be arbitrarily chosen and its influence on other parameters is not investigated in this work.

The macroscopic variables in the CM are interpreted within the framework of the volume averaging method [START_REF] Whitaker | Simultaneous heat, mass, and momentum transfer in porous media: A theory of drying[END_REF] , where the macroscopic variables are expressed as volume averages of the corresponding microscopic variables over a representative averaging volume. For the comparison of the PNM results with the CM solution, the data obtained from PNM simulations are therefore transformed into macroscopic data by volume averaging. Hence, the three-dimensional PNM domain is vertically split into local averaging volumes or slices where each slice consists of vertical throats and their bottom side neighbor horizontal throats. An exemplary 3D pore network is illustrated in Fig. 1 , where a slice of thickness z is indicated. For more details on the method of slicing and the macroscopic parameter or variable extraction from PNM results, see [ 39 , 40 ] and references therein.

PNM simulations

The three-dimensional PNM simulations are carried out on a network that consists of 25 computational nodes in each spatial direction. The pore size distribution is characterized by a mean throat radius of 250 μm with a standard deviation of 25 μm and a uniform throat length of 1 mm. The network porosity is 0.59 and drying occurs at room temperature and pressure. The average total drying flux ˙ j (which stays nearly constant with time as our study is limited to CRP) in the PNM and CM simulations is ing the total network height as the characteristic length L and the initial effective solute diffusivity D s * (1). In the simulations, Pe = 0.2 at the beginning of drying. This is relatively small but sufficient to have a noticeable advection effect. The advective transport becomes more dominant with time because the effective solute diffusivity decreases with the decrease in saturation.

We perform multiple PNM simulations with different realizations of throat size distribution. Though the network saturation at the end of CRP may vary among simulations, it can be reliably considered that the CRP lasts in fact down to a network saturation of 0.7. Hence, in the CM simulations, the analysis is limited to a network saturation of 0.7. We have considered an arbitrarily chosen value of 10 kg/m 3 as the initial solute concentration for all simulations in this work. In the capillary-dominated regime, the depletion of liquid from the network occurs such that the continuity of the liquid phase decreases as the initially connected bulk liquid phase is split into clusters and isolated single throats [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF]. At any instant during the drying process, the largest liquid cluster with respect to volume is characterized as the main cluster. All other liquid clusters are characterized as isolated clusters, because they are isolated from each other and more importantly from the main cluster which spans over the entire network height. Each isolated liquid cluster eventually breaks down to several single isolated liquid throats which are the smallest structural element of the liquid phase. This increasing discontinuity of the liquid phase directly influences the solute transport because solute transport can only occur between the connected liquid phase elements (as liquid films are not considered). Hence, it is important to represent the solute concentration with respect to the liquid phase elements, i.e. the main cluster, the isolated clusters and isolated single throats. In Fig. 2 we present the instantaneous mapping of the solute concentration for a 2D pore network of 100 × 100. The solute concentra-tion gradient in the main cluster can be seen in Fig. 2 a. The progressive fragmentation of the liquid phase from the main cluster into isolated clusters is also shown in the cut-outs in Figs. 2 b-d, where the solute concentration increase that develops as a result of shrinking and further fragmentation of an isolated cluster can be clearly observed.

In Fig. 3 we present the instantaneous local solute concentration profiles and the corresponding saturation profiles for regular intervals of network saturation ranging between 0.9 and 0.7. The profiles are obtained using the slice averaging method discussed briefly in Section 3 . For the sake of simplicity, the isolated single throats are considered together with isolated clusters as part of the isolated liquid phase. The total liquid phase saturation in a slice is computed as

S = i = n i =1 S ti V ti i = n i =1 V ti , ( 8
)
where n is the number of throats in the slice, S ti is the saturation in throat i ( S ti = 1 if the throat is fully filled with liquid, S ti = 0 when the throat is empty), V ti is the volume of throat i. Similarly, the main cluster saturation is computed as

S mc = i = n mc i =1 S ti V ti i = n i =1 V ti , ( 9
)
where n mc is the number of throats of the main cluster in the considered slice. The isolated cluster saturation can be simply computed as S ic = S -S mc . The slice averaged solute concentration C in the total liquid is computed as

C = i = n i =1 C i S ti V ti i = n i =1 V ti . ( 10
)
The slice averaged concentration in the main cluster C mc is computed as While the slice averaged solute concentration in isolated clusters C ic is computed as

C mc = i = n mc i =1 C i S ti V ti i = n mc i =1 V ti . ( 11
)
C ic = i = n ic i =1 C i S ti V ti i = n ic i =1 V ti , ( 12
)
where n ic refers to the number of throats of isolated clusters in the considered slice. The saturation profiles in Fig. 3 reveal that the proportion of liquid in the isolated liquid elements increases for lower network saturations. At the network saturation of 0.7, the proportion of liquid in isolated liquid elements is approximately equal to that in the main cluster except for a few slices furthest away from the evaporative surface. Another important phenomenon is the presence of strong edge effect [ 41 , 42 ] (a sharp drop in saturation profiles) adjacent to the evaporative surface. A less marked edge effect is also visible in the region adjacent to the bottom surface. The focus being on what happens at the evaporative surface, the edge effect only refers to the vicinity of the evaporative surface in this study. In the early stages of drying, the number of meniscus throats per slice is the highest in the surface slice compared to all other slices in the network, since the network is open to evaporation only from the top. Consequently, it is more likely that the meniscus with the largest radius is located in the surface slice, leading to higher degree of fragmentation of the liquid phase at the surface. Additionally, the higher evaporation rate at the surface and the lack of liquid compensation by capillary pumping leads to fast shrinking of the isolated clusters and single throats. It must be noted that such an edge effect is not expected to be visible in the saturation profile measurements reported in the literature. In the pore network simulations, the edge effect is typically in a region with a thickness around 4-5 lattice spacings (the lattice spacing is the distance between two nodes in the network). This is much too small to be detected by the NMR set-ups used to measure the saturation profiles, e.g. [START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF] . In addition, the interpretation of the NMR signal is much less obvious in the vicinity of the porous medium surface, where measurement artefacts might affect the saturation measurements. Another important difference in Fig. 3 with the experimental profiles, e.g. [ 30 , 43 ], is that the PNM profile is not flat for high network overall saturations, i.e. in Fig. 3 a ( S net = 0.90) and Fig. 3 b ( S net = 0.85), whereas the profiles are typically flat for similar global saturations in the experiments [ 30 , 43 ]. As discussed in [START_REF] Bray | Three-dimensional pore network simulation of drying in capillary porous media[END_REF] , this is a finite size effect due to the network small size. The profiles tend to get flat (except in the region of the edge effect) only after the breakthrough (BT), which is the moment when the gas phase reaches for the first time the porous medium opposite side. As discussed in [START_REF] Wilkinson | Percolation effects in immiscible displacement[END_REF] , the liquid saturation at BT scales as 1 -S BT ∝ N -0 . 48 , where N is the network size ( = the number of pores in one spatial direction). In our simulation N = 25 and S BT ≈ 0.9, whereas in the experiments N >> 25 and thus S BT is expected to be very close to 1. Thus, the period before breakthrough is typically not documented in the experimental works and is typically not captured by the CM which predicts a flat saturation profile right from the beginning of the capillary regime. The edge effect is not predicted by the commonly used continuum model either, i.e. Eq. ( 1) .

As can be seen from Fig. 3 , the edge effect in the saturation profiles can be observed in the main cluster as well as in the isolated liquid elements. The instantaneous solute concentration profiles in the main cluster and for the total liquid phase have the typical exponential-like shape resulting from a dominant convective effect (the back-diffusion is not sufficient to spatially equalize the concentration) [START_REF] Huinink | How ions distribute in a drying porous medium: A simple model[END_REF] . As discussed in some details in [START_REF] Pel | A simplified model for the combined wicking and evaporation of a NaCl solution in limestone[END_REF] , two main effects contribute to the concentration build-up in the porous medium top region: the advective transport of the ions toward the surface where they accumulate and the saturation decrease (the concentration in a shrinking liquid volume necessarily increases). In the main cluster, both effects contribute to the concentration build-up since the main cluster spans the network, but the situation is actually subtler because the change in the main cluster saturation is due to two effects: the loss of water due to evaporation and the fragmentation. Isolated clusters initially form by detachment from the main cluster. The loss of water due to evaporation increases the concentration in the main cluster since again the same amount of salt is contained in a smaller volume of solution. The fragmentation affects only marginally the concentration in the main cluster (i.e. only in the vicinity of the throat whose invasion leads to an isolated cluster detachment). By contrast, the convective accumulation effect is expected to be limited in the isolated liquid clusters since there is no transport over a large distance. It is therefore expected that the change in concentration in the isolated clusters is mostly due to the cluster shrinking effect due to evaporation (as illustrated in Figs. 2 c and2 d). This is consistent with the shape of the isolated cluster concentration profiles in Fig. 3. Owing to the so-called screening effect, i.e. the fact that the evaporation rate is very low inside the medium and only significant in the network top region, there is limited change in isolated cluster concentrations away from the top surface. This can also be observed in Fig. 2 a on the left where the concentration in the isolated clusters that are screened from evaporation is the lowest (as these got fragmented from the main cluster in the early stages of drying and then got screened from evaporation by the main cluster on the right side and the isolated clusters and single liquid throats on the top). The concentration increase in the screened liquid clusters is mainly due to the birth of new isolated clusters with slightly higher concentration from the main cluster.

The abrupt increase in isolated cluster concentration in the network top region in Fig. 3 results from a combination of at two effects. First, this is the evaporation active region. As a result, the isolated clusters shrink (as illustrated in Figs. 2 c and2 d), which again leads to a concentration increase. Second, new isolated clusters are formed from the main cluster with a relatively high concentration (as illustrated in Figs. 2 b and2 c), since the concentration in the main cluster is high in this region. In this respect, the edge effect contributes to the concentration build-up since the saturation decrease is faster in the top region than elsewhere in the medium. Another interesting aspect lies in the fact that the concentration in the main cluster is slightly higher than the average concentration in the total liquid phase, except in the edge effect region where the main cluster concentration is comparable to or lower than in the total liquid phase. In particular, the concentration maximum at the surface is observed in the isolated clusters. This is a clear indication that the crystallization must start in isolated throats or clusters connected to the surface. However, since the corresponding volume of solution can be small, this does not lead necessarily to a significant crystal development. This point is addressed in more details further in the paper.

The concentration profiles depicted in Fig. 3 are computed as indicated by Eqs. ( 9) to [START_REF] Whitaker | Simultaneous heat, mass, and momentum transfer in porous media: A theory of drying[END_REF] and are therefore instantaneous solute concentration profiles. In [START_REF] Moghaddam | Kinematics in a slowly drying porous medium: Reconciliation of pore network simulations and continuum modeling[END_REF], it was shown that the spatial average, i.e. the slice averaging, should be combined with a time average to obtain meaningful macroscopic velocity profiles from the PNM velocity field. Thus, the question arises as to whether a similar time-averaging procedure should be combined with the slice averaging procedure as regards the concentrations. We do not think that this is necessary because the instantaneous concentration profiles actually integrate the time-averaging effect on the velocity field, i.e. the evolution of the velocity field in the network as a function of time up to the time the three concentration profiles are computed. Nevertheless, one can argue that the instantaneous profiles can be misleading in situations where the solute concentration in a single isolated throat increases sharply as the result of the liquid mass loss due to evaporation. However, such a single isolated throat with very high solute concentration usually corresponds to little liquid volume that consequently empties rapidly. Hence this instantaneous phenomenon does not result in a sharp fluctuation in the slice concentration. Since the local solute concentration is computed by considering liquid contained in all the throats in a slice, the contribution of the liquid and solute mass in a single throat is not significant. Moreover, the time for evaporation of the single isolated throats is relatively small, hence its influence on the time-averaged solute concentration profile would be small as well.

Solution of classical CM in comparison with PNM results

The only missing parameter for solving the commonly used CM summarized in Section 2 is the effective solute diffusion coefficient D s * . The method used for computation of this parameter for pore network simulations is explained in this section.

Computation of effective solute diffusivity from PNM simulations

In this section, we compute the macroscopic solute diffusivity using the volume averaging method described in Section 3 . In PNM simulations, the pore-scale solute diffusivity is fixed; however, similar to the classical CM formulation of effective diffusivity, the effective solute diffusivity computed from PNM simulations is a function of local saturation. Macroscopic solute diffusivity is computed from PNM simulation results based on

D * s ( S ) = j s,di f f ρ l εS∂ C/∂ z , ( 13
)
where j s,diff , S , ∂ C/∂ z denote the macroscopic solute diffusive flux, the slice average saturation and the gradient of the slice average concentration along the spatial coordinate z , respectively. The porescale solute diffusive flow rate across a vertical throat k is computed by

J s,di f f,k = A k D s ∂C ∂z , ( 14
)
where A k , D s , and ∂C ∂z denote the cross-section area of the vertical throat k , the pore-scale solute diffusivity and the pore-scale concentration gradient across throat k , respectively. Using Eq. ( 14) , the solute diffusive flow rate across all the vertical throats in a slice is computed and summed up to obtain the macroscopic diffusive flow through the slice. This macroscopic diffusive flow rate is divided by the macroscopic cross section area of the slice to obtain the macroscopic solute diffusive flux j s,diff which is then used to compute the macroscopic solute diffusivity based on the macroscopic concentration gradient through the slice in consideration ( Eq. ( 13) ).

As Fig. 4 illustrates, the effective solute diffusivity computed from PNM simulation results is nearly a linear function of local saturation, consistently with the usual functional form given by Eq. ( 4) . The data points in Fig. 4 have been computed by averaging the values of effective solute diffusivity for intervals of slice average saturation equivalent to 0.005. In fact, there are values of effective solute diffusivity for local saturation smaller than 0.7, because as shown in the saturation profiles of Section 4 , the local saturation is smaller in the region near the evaporative surface as compared to the rest of the network. However, in CM the local saturation is assumed to be uniform in space, hence, at any instant, the network saturation and local saturation are equal. Moreover, as 4) ). mentioned early in Section 4 , we consider a network saturation of 0.7 as the lower threshold value of CRP. Therefore, in Fig. 4 , we do not consider values smaller than the threshold value of 0.7.

Fig. 4 also shows the profile obtained from fitting the parameters using the classical formulation expressed by Eq. ( 3) . This fitted profile is based on D * s (1) = ε 0 . 86 D s and S c = 0.69. As can be seen, Eq. ( 4) is quite well adapted to represent the PNM computed effective diffusion coefficient.

Solution of CM

In this section, we employ the effective solute diffusivity obtained from PNM simulation results to solve the CM. The solution of CM is compared with the PNM slice averages in Fig. 5 . The CM solution complies well with the concentration profiles obtained from PNM simulation results except for the region near the evaporative surface. For the comparison with CM solution, PNM results are considered for the whole bulk liquid regardless of the distinction between the main cluster and the isolated clusters.

As can be seen, the CM solution underestimates the concentration at the surface. This means, for instance, that the CM model should predict the crystallization onset later than the PNM solution. This might explain why the time of first crystallization determined from CM simulation is longer than in the experiments [START_REF] Hidri | Porous medium coffee ring effect and other factors affecting the first crystallisation time of sodium chloride at the surface of a drying porous medium[END_REF] . Based on the non-compliance in CM solution adjacent to and at the evaporative surface, it is evident that in order to further improve the CM solution, an investigation dedicated to studying the dynamics of solute transport at the surface must be performed. In this regard, the pore scale insight offered by PNM simulations can be exploited. To this end, we present a detailed analysis focused on transport at the porous medium surface in the following sections.

Liquid phase connectivity and solute transport at the surface

In Fig. 6 , we present the instantaneous solute concentration and the liquid connectivity mapping for the liquid in vertical surface throats for various network saturation values. Note that the results presented in this section are based on one PNM simulation, as opposed to the results presented in Section 4 which are averages of multiple simulations, each with a different realization of throats size distribution.

The liquid connectivity plots for various network saturation values presented in Fig. 6 reveal that the majority of surface throats belong to the main cluster. Moreover, there is no distinct trend in the variation of the proportion of isolated liquid phase with decrease in network saturation. This is because, as mentioned in Section 4 , the surface throats that belong to isolated liquid phase get evaporated quickly because of their limited volume (and hence limited supply of liquid through capillary pumping) and higher evaporation rate on the surface (no evaporation screening effect because of direct exposure to the boundary layer). The proportion of surface throats that belong to isolated liquid phase is an important factor as the solute concentration values in the isolated clusters and isolated single throats are relatively high as compared to that in the main cluster, regardless of the network saturation value.

The solute concentration maps presented in Fig. 6 reveal that the solute concentration fluctuates spatially in the main cluster at the surface. This is more visible for the network saturations S net = 0.9 and 0.8. For the lower saturations, the concentration fluctuations in the main cluster are hidden by the quite large concentration variation range due to the very high concentrations reached in a few isolated throats. For this reason, the concentration distribution in the main cluster at the surface for S net = 0.75 and 0.7 are more specifically shown in Fig. 7 . As discussed in [START_REF] Veran-Tissoires | Evaporation of a sodium chloride solution from a saturated porous medium with efflorescence formation[END_REF] for the simpler evaporation-wicking situation, two main factors explain the concentration fluctuations in the main cluster. The first one is referred to as the internal disorder effect. Due to the internal disorder of the pore space, i.e. the fact that the throat sizes vary randomly, the pore-scale velocity field also fluctuates spatially. This means that the mean velocity in a throat varies from one throat to the other. Also, the velocity field in the throats fluctuates in time due to the main cluster mass loss and fragmentation. As a result of these fluctuations, the concentration, which greatly depends on the advective transport, also fluctuates. The second factor is referred to as the surface disorder effect. As a result of the throat size variation at the surface and the formation of wet patches of different sizes (a wet patch refers to a 2D cluster of saturated throats [START_REF] Moghaddam | A pore network study of evaporation from the surface of a drying non-hygroscopic porous medium[END_REF] at the surface which is surrounded by gas throats), the evaporation rate varies from one wet throat to the other at the surface. These evaporation rate variations induce velocity variations in the vertical wet throats connecting the network to the surface, which, in turn, induce concentration variations. In particular, the spacing between the wet patches at the surface plays a key role on the diffusive vapor flux per pore (see [START_REF] Chauvet | Three periods of drying of a single square capillary tube[END_REF] ). Similarly, a dry patch refers to a 2D cluster of gas throats surrounded by liquid throats. Larger dry patches between the wet patches may result in higher evaporation fluxes through the wet patches, therefore higher solute concentration pores that have higher evaporation rate. The liquid patch on the mid-left of Fig. 7 , for network saturation of 0.75, is an example of the high solute concentration as a result of high evaporation rate because it is surrounded by large dry patches.

The solute concentration maps for network saturations of 0.75 and 0.70 presented in Fig. 6 show that solute concentration can rise to very high values in the isolated liquid phase compared to the main cluster. This is due to two main factors in case of single isolated throats. First the lifespan of such a throat is quite short since the volume of solution in it is small and evaporation is intense at the surface. Thus, the concentration is already high in such a throat when it forms, since it forms as the result of fragmentation of the main cluster or an isolated cluster where the concentration is high. Then the volume of solution in the throat rapidly decreases due to evaporation, which leads to a sharp increase in solute concentration. A similar process occurs in isolated clusters connected to the surface, with in addition subsequent fragmentations and increases in the concentration due to the solution volume shrinking effect.

Statistical analysis of transport at the surface based on PNM Monte Carlo simulations

As presented in Section 6 , the evolution of solute concentration at the surface is strongly linked to the evolution of liquid structure at the surface. However, the heterogeneity in liquid phase structure is strongly dependent on the throat size distribution which varies in each PNM simulation. Therefore, in this section, we characterize stochastically the evolution of liquid phase structure and solute concentration at the surface. Moreover, as the main purpose of the PNM simulations is to provide a benchmark for the CM solution, a key parameter in this comparison can be the onset of crystallization based on a threshold value of solute concentration. Due to the averaging of data on the scale of local slices and the subsequent averaging of these results among multiple simulations with different realizations of throat size distribution, it is possible that the onset of crystallization on the pore scale in a PNM simulation may not be visible in the resulting instantaneous saturation and solute concentration profiles presented in Section 4 . Therefore, we need to perform a dedicated study at the pore scale. For this purpose, we conduct a thorough analysis based on Monte Carlo simulations by performing 70 PNM simulations with varying realizations of the throat size distribution. Similar to the surface analysis presented in Section 6 , the statistical analysis is focused only on the vertical surface throats.

Based on the data obtained from the Monte Carlo simulations, we characterize the liquid phase heterogeneity with respect to varying network saturation values and present the results in the form of a histogram. As indicated in Fig. 8 , nearly all of the saturated surface throats are likely to be part of the main cluster until S net = 0.8. Also, for network saturation of 0.75, the probability that a saturated surface throat belongs to the main cluster is still 0.97. However, when network saturation reaches 0.7, the probability that a saturated surface throats belong to the main cluster decreases to 0.917, whereas that of belonging to isolated clusters and single isolated throats is 0.052 and 0.031, respectively.

Apart from the stochastic characterization of liquid connectivity of surface throats, Fig. 8 also illustrates the variation of surface slice saturation with network saturation. Interestingly, for the initial 10 % decrease in network saturation, there is a dramatic decrease of approximately 38 % in saturation of the surface slice. Afterwards, the decrease in surface slice saturation is for further 10 % decrease in network saturation, i.e. from S net = 0.90 to 0.80, is relatively small, i.e. S sur decreases from 0.62 to 0.50. As network saturation decreases from 0.80 to 0.70, the surface slice saturation decreases from 0.50 to approximately 0.30. Overall, the decrease in surface slice saturation is the highest as network saturation de- creases from 1 to 0.90. The edge effect is therefore quite significant with significantly less saturation in the top slice compared to the bulk.

Fig. 9 illustrates the probability of reaching of specified threshold concentration value, denoted by C th . In this example, C th is equal to 50 kg/m 3 whereas the initial concentration is as mentioned before C 0 = 10 kg/m 3 . Thus, the threshold value corresponds here to an increase of the concentration by a factor 5 compared to the initial concentration. The line plot (corresponding to the y-axis on the right) represents the absolute probability of reaching the threshold value in a saturated surface throat (absolute probability does not distinguish between the main cluster and isolate liquid phase). As can be observed, the probability of reaching the threshold value for S net values of 0.90 and 0.85 is extremely low, while the probability of reaching the threshold value is higher in single isolated throats as compared to the main cluster (as illustrated by the respective histograms). In other words, if one assumes for instance that the ratio C th / C 0 = 5 corresponds to the onset of crystallization, Fig. 9 shows that the probability of crystallization in this example for S net values of 0.90 and 0.85 is quite low and that the crystallization, when it happens, occurs with a much greater probability in an isolated single throat than at the top of the main. As network saturation decreases from 0.75 to 0.70, the absolute probability of reaching the concentration threshold value for all saturated vertical surface throats increases significantly from At intermediate network saturation values, e.g. 0.80, the probability of reaching the concentration threshold value in isolated liquid phase approximately matches that of the main cluster. It is obvious that the commonly used CM cannot capture these effects which reflect the impact of heterogeneities at a scale lower than the Representative Elementary Volume (REV) scale typically associated with the macroscopic modeling.

Discussion

At first glance, the results reported in Fig. 9 could be seen as in contradiction with the solute concentration profiles depicted in Fig. 3 . Fig. 3 shows that the greater concentration at the surface is observed in the isolated cluster for all the considered network saturations. On the other hand, Fig. 9 shows that the concentration threshold value at the surface is reached with a greater probability in the main cluster when the network saturation is sufficiently low, for example S net = 0.7, in the case considered. The latter result simply reflects the fact that, as illustrated in Fig. 6 , there are much more surface throats belonging to the main cluster than to isolated clusters. For instance, if the concentration threshold was selected such that the solute concentration is greater than the considered concentration threshold value in all liquid throats at the surface, the probabilities plotted in Fig. 9 would be simply given by μ abs -mc = number of sur f ace throats in main cluster number of liquid sur f ace throats [START_REF] Veran-Tissoires | Evaporation of a sodium chloride solution from a saturated porous medium with efflorescence formation[END_REF] μ abs -ic = number of sur f ace throats in isol ated cl usters number of liquid sur f ace throats [START_REF] Hidri | Stagnation points as loci of solute concentration extrema at the evaporative surface of a random porous medium[END_REF] with μ rel-mc = μ abs -mc , μ rel-ic = μ abs -ic for this particular case.

In terms of probability of crystallization during the CRP, the results can be summarized as follows. As shown in Fig. 3 , the crystallization should systematically start in the isolated clusters. This is because the fast evaporation of the isolated clusters connected to the surface leads to a rapid increase in the solute concentration. However, these clusters can be small or even can simply be liquid bridges (the single isolated throats in our simulations). As a result, this crystallization is likely to correspond to a relatively small amount of crystals. As a result, the formation of these crystals could only marginally affect the evaporation process. After a while, the crystallization concentration should be reached in the main cluster at the surface (as indicated by Fig. 3 and Fig. 9 ). Then, the development of crystals (efflorescence) can be expected to be much more important since the main cluster corresponds to a much greater amount of salt compared to an isolated cluster. A detailed analysis of the above dynamics would however require to include the crystallization process in the modelling and to couple the efflorescence development with the modeling of evaporation and solute transport processes. This is a quite challenging objective, which has not yet been addressed in the framework of pore network modeling.

As mentioned in Section 3 , we consider the pores as volumeless computational nodes. This assumption leads to an overestimation of the porous medium porosity (through the overlap region of all neighboring throats at each pore). However, as we consider the capillary dominated regime, the consideration of pore body volume does not impact the liquid fragmentation dynamics because based on the geometric constraints, the radius of the pore body is always larger or equal to its neighbor throat radii. This means that the invasion percolation dynamics is the same in PNM with and without pore body volumes. As the liquid fragmentation dynamics is not affected, the use of PNM without the pore body volume can be considered to be adequate in the capillary dominated regime. On the other hand, considering the pore body volumes can lead to lower irreducible saturation and a longer CRP if the mean throat radius is sufficiently high (see [START_REF] Lu | Transport parameters of macroscopic continuum model determined from discrete pore network simulations of drying porous media[END_REF] ). Also, the mixing process taking place in the pores might depend on the pore body volume. Perfect mixing, i.e. a uniform solute concentration, is assumed in the pore bodies in our model. Although, we do not believe that the consideration of a pore network with pore body volumes would change the main conclusions of the present work, this aspect would deserve further investigations.

Conclusions

The drying pore network simulations reported in this study show two striking differences compared with the predictions of the commonly used continuum model as regards the saturation profiles during the CRP in the capillary dominant regime, which was the focus of this study. First the saturation profiles become flat in the bulk only after an initial period which is typically not seen with the CM, nor in the experiments, e.g. [START_REF] Gupta | Paradoxical drying due to salt crystallization[END_REF] . Second a quite significant edge effect is observed leading to a significantly reduced saturation in the top region of the network compared to the bulk. The ratio between the saturation in the bulk and at the surface can be greater than 2. The first effect is a finite size effect due to the small size of the considered network compared to the typical size of the porous medium in experiments. Therefore, this effect needs not specific consideration since it is expected to be negligible with usual porous media. However, it might need consideration in applications involving thin porous media, such as in fuel cells for example [START_REF] Prat | Thin porous media, Chapter 4[END_REF] . The second effect, ref erred to as the edge effect, is more problematic. Unlike the first effect, preliminary simulations (not shown in this study) suggest that it is not network size dependent.

The comparison between the predictions of the commonly used continuum model (CM) and the pore network model as regards the solute transport shows that the CM underpredicts the concentration at the evaporative surface compared to the PNM simulations. This can be mainly attributed to the edge effect since the saturation variations are one of the important factors controlling the evolution of solute concentration. The greater saturation variation in the edge region leads to greater solute concentration compared to the predictions of the CM. This is an indication that the traditional continuum model should overpredict the crystallization onset time at the porous medium surface.

Another important feature, which is not captured by the traditional CM, is related to the liquid phase structuration in liquid clusters. The PNM simulations show that the dominant mechanisms controlling the solute concentration variations in the porous medium top region are different in the percolating main cluster and the isolated clusters. The solute concentration variations are mainly due to the saturation variations in the isolated clusters and the combination of advective transport and saturation variations in the main cluster. The solute concentration at the surface is greater in the isolated clusters than in the main cluster. This suggests that an improvement should be to develop a drying continuum model making an explicit distinction between the main cluster, i.e. the percolating liquid phase, and the isolated cluster, i.e. the non-percolating liquid phase. This type of model has been proposed for the modelling of two-phase flow in porous media, e.g. [START_REF] Doster | Numerical solutions of a generalized theory for macroscopic capillarity[END_REF] , but has not yet been extended to the drying situation.

In this article, the macroscopic data were obtained from the PNM simulations by volume averaging considering relatively thin slices as averaging volumes. Then the question arises as to whether the slices can be considered as a truly representative averaging volume. In this respect, it would be interesting to consider thicker averaging volumes and perform simulations over larger networks. This is not possible with the current version of the drying PNM code used in this article. Although the traditional length scale separation criterion is therefore not fully met in our simulations, we however consider our work to give valuable insights into the shortcomings of the conventional continuum model.

Fig. 1 .

 1 Fig. 1. A sample three-dimensional pore network illustration with evaporative surface on top. The boundary layer of thickness L bl is shown on the top. A slice of thickness z (local averaging volume) that consists of one row vertical throats and the interconnected bottom side horizontal throats is also illustrated. The throats filled with liquid are represented in blue, the gas-filled throats in gray and all pores in red color. Three single isolated liquid throats at the surface are also highlighted.

s 3 .

 3 025 × 10 -6 kg/m 2 /s. The competition between the advective and diffusive solute transport is characterized by the Péclet number defined as Pe = j ˙L D * . The initial Péclet number is computed by usε ρ l

Fig. 2 .

 2 Fig. 2. A 2D PNM representation of the instantaneous distribution of solute in main cluster and isolated clusters for a network of 100 × 100 pores. The evaporative surface lies at the top of the network. The gas phase is shown in white in all plots, whereas the solid phase is shown in grey in (a) and in white in (b), (c) and (d). The subplot (b) is a cut-out from (a) and illustrates a section of the main cluster. The subplot (c) represents the splitting of an isolated cluster from its parent main cluster illustrated in (b), whereas (d) indicates the subsequent fragmentation and shrinking of the isolated cluster shown in (c). The color bar on the left indicates the range of solute concentration. Note that results of 3D PNM simulations are presented in the rest of the paper.

Fig. 3 .

 3 Fig. 3. Left: Saturation profiles for main liquid cluster, isolated liquid clusters including single isolated throats, and total liquid. Right: Solute concentration profiles for the respective liquid elements. The legends illustrated in (a) apply to the rest of the rows. Plots in (a) to (e) represent results for network saturation values of 0.90, 0.85, 0.80, 0.75 and 0.70, respectively. The profiles are averaged for data obtained from 15 PNM simulations, each with different realization of throat size distribution.

Fig. 3 .

 3 Fig. 3. Continued

Fig. 4 .

 4 Fig. 4. Black square symbols represent the effective diffusivity obtained from PNM simulation results. The dashed blue line indicates fitted profile based on parameter adjustment of the classical formulation of effective diffusivity ( Eq. (4) ).

Fig. 5 .

 5 Fig. 5. CM solution based on effective diffusivity extracted from PNM simulation results compared with corresponding PNM simulation results. The profiles are presented for network saturations of 0.9, 0.85, 0.80, 0.75 and 0.70. The arrow indicates the trend of the profiles with decrease in overall network saturation S net .

Fig. 6 .

 6 Fig.6. The column on the left shows the liquid connectivity for vertical surface throats, where the red color represents throats that belong to the main cluster, green color represents the isolated single throats and yellow color represents the throats that belong to isolated liquid clusters. The column on the right maps the solute concentration for the corresponding throats, where the legend (vertical bar) indicates the range of solute concentration values in the throats. Plots in (a) to (d) represent results for network saturation values of 0.90, 0.80, 0.75 and 0.70, respectively. The size of the throats is scaled up for ease of visualization.

Fig. 7 .

 7 Fig. 7. Solute concentration maps for main cluster surface throats for network saturations of 0.75 and 0.70. The saturated throats that belong to isolated liquid phase are not shown. Legend (vertical bar) indicates the range of solute concentration values in the throats.

Fig. 8 .

 8 Fig. 8. Statistical analysis of the liquid phase heterogeneity based on Monte Carlo PNM simulations. The bars quantify the probability of a saturated surface throat being part of the main cluster, part of an isolated cluster or a single isolated throat. The line plot (corresponding to y-axis on the right) quantifies the surface slice saturation with respect to different network saturation values.

Fig. 9 .

 9 Fig. 9. Statistical analysis quantifying the probability of reaching a concentration threshold value equal to 5 times the initial concentration. The line plot (corresponding to y-axis on the right) quantifies the absolute probability of reaching the concentration threshold value in surface throats with respect to different network saturation values. The bars (corresponding to the y-axis on the left) quantify the probability of reaching the concentration threshold value within the respective liquid phase elements with respect to the total number of throats in which the concentration threshold value is reached.

Acknowledgement

This work was financed by joint project "Drycap" funded by DFG (project TS28/10-1) and GIP ANR (project 16-CE92-0030-01).

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

CRediT authorship contribution statement