Effects of model boundary conditions on simulated drying kinetics and inversely determined liquid water permeability for cement-based materials

Abstract The accurate simulation of moisture transport in cement-based materials is an important step to predict many durability processes of concrete structures. This paper studied the effects of boundary conditions on modeling moisture transport, in particular for simulating drying kinetics and the inverse determination of liquid water permeability. Three boundary conditions were investigated, namely, the Dirichlet (imposed) condition, the flux boundary condition, and the boundary supplied by moisture movement in the ambient environment (coupling with a convection-diffusion moisture transport model). By comparing with the measured drying mass loss curves (without airflow) for ordinary Portland cement pastes, we found that the Dirichlet boundary condition provides slightly better simulated drying kinetics than the flux boundary condition, while results of the coupling model are less promising. As for inversely determined liquid permeability, the studied boundary conditions lead to the max. 36% difference. If considering the airflow in the environment, both the flux boundary condition and the coupling model with computational fluid dynamics (CFD) show that the airflow does not have an effect on moisture evaporation for the studied cement-based materials, because of their low permeability, resulting in a very short initial stage of drying. Moisture evaporation for these materials is mainly controlled by the internal moisture transport. However, for materials with high permeability, the drying kinetics increase with the airflow velocity because the initial stage of drying becomes much longer than the low permeable materials. The implication is that some laboratory methods for drying cement-based materials (specimens stored in a sealed chamber without airflow) may need to be reevaluated because the relative humidity at the material surface is different from the controlled relative humidity.


Introduction
In the natural conditions, the concrete cover of steel reinforced concrete structures is rarely saturated with water, [1] meaning that there are pore spaces available for both liquid transport and gas/vapor diffusion. Liquid water movement may carry aggressive agents, such as chloride ions, into concrete and thus leads to corrosion of steel reinforcing bars (rebars). [2][3][4] Carbon dioxide can diffuse into concrete through the empty pores and cause concrete carbonation if concrete is unsaturated with water, [5] and later causes steel corrosion once water content in concrete reaches a high level. For these reasons, the accurate simulation of moisture transport is crucial for studying the durability of reinforced concrete structures.
The unsaturated moisture transport in porous materials can be described by models with different levels of complexities. The multi-phase model for cement-based materials, including transport of liquid water, water vapor and dry air, has shown a good ability to provide the simulated drying mass loss curves matching well the measured ones. [6,7] The multi-phase model can be simplified to a two-phase model, which only considers liquid water advection and vapor diffusion. [8,9] The two-phase model was further reduced to a single-phase model with transport of liquid water. [7] However, it was found that this simplification only works for moisture transport at high relative humidity (RH). [9] At low RH, the single-phase model underestimates the amount of moisture transport. [10,11] In addition to the complexities of models, other factors that can affect the simulation accuracy were investigated in the literature, such as the microstructure of cement-based materials and the liquid water permeability k l . The swelling of hydration products during liquid water uptake has been observed by environmental scanning electron microscopy more than two decades ago. [12] The microstructural rearrangement during drying was also reported in some experimental results (e.g., [13] ). Other causes, such as dissolution (during water uptake) and precipitation (during drying) of solids, were also investigated in the literature by using the modeling approaches (see an example in Ref. [14] However, the dynamic change of the microstructure during moisture transport was not fully considered in the conventional models, which results in the low simulation accuracy for some measurements, in particular for the long-term measurements. These results that cannot be simulated by conventional models have been termed as anomalous moisture transport. [15,16] The current approaches in the literature considering the microstructural effect on moisture transport is directly to write k l as a function of time, in which, the theory is that transport coefficients are directly related to materials' microstructure and the change microstructure is also a function of time; this group is known as the time-dependent models. [17][18][19][20] Another group considers the different transport speeds of moisture in different sizes of pores, so-called dual-permeability models. [21] Both of them can provide high accuracy simulation results. The liquid water permeability k l , generally consisting of two parts, intrinsic permeability K l and relative permeability k rl (see Ref. [7] and Section 2.2), is a key factor for a moisture transport model to predict the mass change and water content in a cementitious material exposed to the natural environmental conditions, in which the material generally endures relatively high water content. [22] These conditions could lead to more serious durability issues than the dry conditions. [23] Previous studies tried to feed the moisture transport model with the directly measured K l by the traditional steady state water flow method [24][25][26][27] and indirectly estimated K l by Katz-Thompson [28] and Kozeny-Carman equations [29,30] and sorptivity method. [31] However, the simulation results were found to overestimate moisture transport compared with experimental observations. [7,32,33] The indirect measurement methods, such as the beam bending method [34][35][36] and thermopermeametry, [37,38] which calculate K l based on the time-dependent deformation of a specimen by a theoretical model, were found to be able to provide K l with better results than the other methods, because these two methods employ saturated specimens whose original microstructure can be better maintained than Katz-Thompson and Kozeny-Carman equations and the sorptivity method. When K l is the only unknown in a moisture transport model, it can be determined by the "inverse analysis" (the inverse modeling from the numerical modeling point of view), which was firstly used for cement-based materials by Coussy [6] and now can be commonly found in the literature (e.g., [1,9,[39][40][41] ). This approach utilizes a numerical moisture transport model to back calculate transport coefficients based on the measured data, such as drying mass loss, moisture distribution, etc.
However, the effects of boundary conditions on moisture transport simulation are rarely discussed in the literature. A boundary condition describes the interactions between the material and the ambient environment. The use of different boundary conditions is expected to change the simulation accuracy and further affect K l if the inverse analysis method is used. There are two common boundary conditions in the literature: the Dirichlet boundary condition (imposed RH at the boundary) and the flux boundary condition. Previous studies have shown that it is important to consider the moisture transport in the ambient environment for drying porous materials such as soils, [42] woods, [43,44] etc. However, for the mature cement-based materials, no study has been performed to examine whether the consideration of the external moisture transport has an effect on moisture transport inside cement-based materials.
The present study is devoted to evaluate the effects of boundary conditions on modeling moisture transport, in particular for the simulated drying mass loss curves and the determination of K l by the inverse analysis method. This study first considers the case without airflow on the material surface. Experimental data, including sorption isotherms and drying kinetics, were obtained for Ordinary Portland cement pastes with three water-to-cement ratios (w/c) to calibrate the model. Secondly, moisture transport in cement pastes is coupled with moisture movement in ambient environment simulated by a convection-diffusion model. A 2 D model is created to consider the effect of airflow, which is provided by a computational fluid dynamics (CFD) model. Finally, the effect of airflow on moisture transport in the cement-based material is discussed.
2. Review of moisture transport modeling approaches

Boundary conditions
A boundary condition is used to take into account the moisture interaction between the porous material and its surroundings. The Dirichlet boundary highly simplifies this interaction for numerical implementation. However, this boundary condition can cause a sharp initial vapor pressure gradient and a sudden change of mass loss. These phenomena are often observed for initially high water content materials. [45] The flux boundary condition views the environment as a thin boundary layer. This layer has the same properties as the porous material and thus the transport parameters from the material part can be used to determine the moisture state in the layer. Even though the flux boundary condition is able to provide good simulation results, it is different to the real condition where moisture flows and diffuses in the air close to the material surface. In addition, experimental measurements revealed that a vapor pressure gradient clearly exists in the environmental region adjacent to the material. [46] This implies that the use of the externally controlled RH as boundary RH may be inappropriate. However, most moisture transport models for cement-based materials in the literature only focus on the mass transport taking place within the material. In the real cases, it is important to couple moisture transport between the porous material and its surrounding. The coupling involves several mechanisms that affect the macroscopic transport, including phase change between liquid and vapor, mass transport by diffusion and convection, etc. Unlike the flux boundary condition only considering effects of external RH on moisture transport in the material (one-way effect), the coupling method can take into account mutual interactions between the material and its surroundings.
Several coupling methods can be found in the literature to take into account the ambient moisture transport. The 1 D-HAM (Heat-Air-Moisture) PCprogram [47,48] created a fictitious surface film that is used to account for the material surface resistance. This film is simulated as a material layer without additional airflow resistance. This method is similar to the flux boundary condition. Building energy models, such as EnergyPlus, [49] TRNSYS, [50] HAM and simplified models (e.g., EMPD), assume that state variables are uniform in the air, so that only one node is considered in the environment part to couple moisture. [49,50] Laurindo and Prat [51,52] considered the presence of a mass boundary region (MBR) over the external surface of the porous material to estimate the drying rates. As the name implies, the MBR is a region close to the porous material surface but has different moisture content to the bulk air. Thus, the coupled system of the porous material and its surrounding can be viewed as three parts: a) the porous material, b) the MBR and c) the part beyond the MBR in the environment (see Figure 1). Actually, most models in the literature did not distinguish part (b) and part (c), [53][54][55][56] and call them together as an atmospheric boundary layer (e.g., [57] ). However, modeling results by Laurindo and Prat [51,52] concluded that the thickness of the MBR x e has a great effect on drying rate. Yiotis et al. [58] further stated that the drying rate is inversely proportional to x e .
These conclusions encourage us to review the description of the boundary condition associated with experiments. In laboratory drying experiments, a specimen is generally placed in a desiccator with constant RH controlled by a saturated salt solution. Generally, the distance between the specimen surface and the saturated salt solution is 5 cm À 15 cm. In this region, air does not flow and only water vapor diffusion occurs. If vapor diffusion in this region has an influence on moisture transport in the material, K l determined by the inverse analysis method can be also altered. For the natural condition with airflow, the drying rate may be enhanced for most of porous materials, but the initial high mass loss can discontinue liquid water in the pore network and then slow down the mass loss in the later drying stages. [59] 2.1.1. Dirichlet boundary condition The Dirichlet boundary condition is most widely used in the early studies. Relative humidity (RH) of the point on the material surface is imposed to be equal to the ambient RH.
where the superscripts "m, 0" and "e" represent the material side and the environment side at the interface, respectively.

Flux boundary condition
The flux boundary condition (also known as convective condition [60] ) is used to account for an imperfect moisture transport between the environment and the surface of the material. The expression is given as, [61] q where q (kgÁ m -2 s -1 ) is the surface moisture flux and S(-) is the degree of liquid water saturation of the surface of the material. The degree of liquid water saturation S is defined as the ratio of the pore volume occupied by water to the total pore volume determined by the difference between saturated and dried states a specimen (see Section 3.1). This boundary condition includes the material property (porosity /), the ambient environment (external vapor pressure P e v ), moisture state in the material near the surface (P m, 0 v and saturation S) and the interactions between the ambient environment and the material (the coefficient E). The term /S takes into account the reduction of vapor diffusion due to the fact that pores are partially occupied by liquid water in the unsaturated material. The coefficient E (s Á m À1 ) is generally determined from experiments. [61,62] For a given environmental condition (constant RH, temperature and atmospheric pressure), the value of E is affected by the airflow velocity v e (mÁ s À1 ). Note that the direction of v e is along with the material surface (perpendicular to the direction of moisture transport). A Menzel's formula is generally adopted [63,64] and we here take the measured results in Ref. [61] However, some recent studies showed that if the airflow is taken into account, the linear relation of the flux boundary condition is not applicable to describe the interactions between the porous material and its surroundings. [9,65] Eq. (3) was estimated based on the measured evaporation rate of a liquid surface (liquid water or fresh concrete with bleed water on the surface). [61,63,64,66] It can be imagined that, with different initial moisture contents, external RHs and material properties, the validity of Eq. (3) may be questionable.
A model which is able to simulate moisture transport in the environment becomes more useful than the flux boundary condition.

Coupling with external moisture transport
Uno [64] analyzed several traditional equations for the drying evaporation from the surface of hygric concretes, and finally concluded that the drying rate is related to the velocity of ambient air and moisture diffusion coefficient at a given temperature. To fully consider the moisture interactions between the material and its surroundings, both airflow and moisture transport (only in the vapor phase) in the MBR must be considered. Vapor transport in the MBR is represented by a convection-diffusion process with a constant diffusion coefficient.
where D v0 (m 2 Á s À1 ) is the free vapor diffusion coefficient in the air and q v (kgÁ m -3 ) is the vapor density.
For a 1 D model, v represents the airflow velocity with the same direction of moisture transport, which is perpendicular to the direction of v e in the flux boundary condition (see Eq. (3)).
In a laboratory condition, specimens are dried in a chamber without airflow; therefore, the air velocity v can be set as zero. Coupling moisture transport in the porous material with that in the MBR can be done in different ways, such single-domain concept and twodomain concept. [67] Since the moisture transport mechanisms in two domains are different, the two domain coupling concept is chosen in this study. For more porous materials (e.g., soils), a transition zone that mainly results from the rough surface can be considered. [68] Nevertheless, the surface of cement-based materials is much smoother than soils, so that a transition zone is very thin and can be neglected. Hence, the case of a sharp interface is employed here. [67] The continuity of moisture flux at the interface must be fulfilled: One should bear in mind that, in the natural conditions, the airflow at the surface of concrete structures is not avoidable. In the MBR, it needs to simulate not only the moisture transport but also the airflow. The airflow (normally treated as free flow) can be simulated by a computed fluid dynamic (CFD) model by solving the Navier-Stoke equations or only Stokes equations for the cases of low Reynolds numbers. [42] The coupling with a CFD model will be discussed later when analyzing the effect of airflow on drying (see Section 4).

Moisture transport model in the porous material
In order to inversely determine k l , an efficient moisture transport model has to be chosen. At the general REV (representative elementary volume) level, the solid of the porous material is assumed to be rigid. Moisture transport in a partially saturated porous medium like cement-based materials is mainly governed by the transport of liquid-water, water vapor and dry air. In a previous study, it has been shown that dry air has very low contribution to mass transport and only causes fluctuated air pressure in the material. [7] This conclusion was also confirmed by the asymptotic analysis performed by Coussy and Thi ery. [69,70] Additionally, considering that the liquid phase remains incompressible, moisture mass balance equations can be represented by a simplified equation, including liquid water and vapor. Hence, here the model selected for the inverse analysis should include the transport of both liquid water and water vapor. The mass balance equation for unsaturated moisture transport is written as, where q l (kgÁ m -3 ) is the density of liquid water, / (m 3 Á m -3 ) is the porosity of the porous material, g (PaÁ s) represents the dynamic viscosity of liquid water, f ðS, /Þ represents the resistance factor for gaseous diffusion and is related to the tortuosity of the pore network and moisture content, and k l ¼ k rl K l with the relative permeability k rl (-) which is generally treated as a function of S. The quasi-equilibrium between liquid and vapor is assumed in this model. The pressure difference between the gas and liquid (P c ¼ P g À P l ) is defined as capillary pressure which is calculated by the Kelvin's law from the known RH (equation can be found elsewhere [7] ). For a given material, there is a characteristic capillary pressure curve representing the relationship of P c as a function of S. In the case of cement-based materials, this curve can be derived from the measured water vapor sorption isotherm. [71] In this study, the capillary pressure curve is expressed by the two-parameter van Genuchten equation (VG): [72] where a (Pa) and m are two fitting parameters. The relative permeability k rl is calculated by the van Genuchten-Mualem (VGM) model. [72,73] where parameter m is taken from Eq. (7) and ' ¼ 0:5 for general porous materials.
The resistance factor f ðS, /Þ represents the reduction of accessibility for water vapor diffusion due to the presence of the liquid phase and the tortuous path for diffusion, etc. It can be formulated as: [74] f ðS, with x D ¼ 2.74. [33] In this moisture transport model, values of empirical parameters (x D and ') are taken from the literature, which have be proved to work well for cementitious materials. Of course, adjusting these values will change the simulation results, while without sufficient experimental data to reevaluate these parameters, we consider them as constants in this study. Therefore, the only unknown is the intrinsic permeability K l which is determined by the inverse analysis.

Experimental data
Cement pastes with three water-to-cement (w/c) ratios were used to calibrate the described model. The material properties are gathered in Table 1. All materials are made from the same OPC cement (CEM I, according to EN 197-1 European standard). All pastes were cast in cylindrical molds with 7 cm in diameter. After curing in the sealed condition for 6 months, specimens were removed from the molds. Cylinders were cut into around 10 cm-length specimens, and the circumferential surface and one end of the specimen were sealed with self-adhesive aluminum foil sheets. The other end was open for moisture exchanges with the ambient environment. All specimens were stored in desiccators for drying experiments at room temperature (T ¼ 25 C) and RH e ¼ 33 or 71% fixed by using saturated salt solutions (drying of P04 at 71% RH was not performed). The drying duration was 141 days for RH e ¼ 71% and 97 days for RH e ¼ 33%. The porosity was calculated based on the mass difference of a small specimen between the vacuum-saturated state and dried state in an oven. The mass loss curves were measured by weighing specimens at different times. Saturation profiles at the end of drying were calculated based on the density changes measured by the gamma-ray attenuation technique. [61,75] The cylindrical specimen was put on a rotating plate with a constant speed while the gamma rays passed through the specimen. By comparing with gamma ray intensities for a specimen at the vacuum-saturated and oven dried states, the degree of liquid water saturation can be determined. The rotating plate can move vertically so that the saturation profile along the specimen is measured.
For the desorption isotherm measurements, a small piece of paste was crushed and stored in limewater prior to measurements. The measurements were conducted in a Dynamic Vapor Sorption instrument (DVS Adventure from Surface Measurement Systems Ltd., UK), which is able to precisely control temperature and RH and to accurately measure the mass of samples. [76][77][78][79][80] The equipment is placed in an incubator, permitting the constant temperature. The required RH is created by mixing the completely dry nitrogen gas with the nitrogen gas with 100% RH. The relative humidity of the mixed gas flow could be varied between 0 and 98%. [81,82] After samples were taken out from limewater, they were immediately transferred into DVS to avoid the further drying and to reduce the risk of carbonation. The desorption isotherm measurements started from RH ¼ 95% and stepwisely decreased to zero: 95 !90 !80 !70 !50 !30 !20 !10 (constant temperature at 25 C). To get a totally dried sample, after drying at 10% in DVS, samples were moved to a N 2 filled chamber. The equilibrium was defined by the relative mass difference criteria (Dm=Dt ¼ 0:0001 %=min) with a fixed maximum duration (120000 s). [80] Even so, we found that some measured mass change curves have not reached the plateau when criteria were fulfilled. To determine the final mass of the sample at each step, measured curves were fitted by a dual-Weibull function and the final mass can be extrapolated. [83]

Models setup
For the implementation of the Dirichlet and the flux boundary conditions, a 1-D model is sufficient. The left boundary of the simulation domain (porous material) is set with the corresponding boundary condition. The right boundary is set as zero flux exchange with the environment.
To implement the coupling model, the state variable for the porous materials needs to be changed to the density of vapor q v which is the same as the moisture transport in the environment. Therefore, the mass balance equation (see Eq. (6)) is rewritten as where R ¼ 8.314 (JÁ K À1 Á mol -1 ) is the gas constant, T (K) is the absolute temperature and M v (kgÁ mol -1 ) is the molar mass of water molecule. The term dS dq v is calculated from the sorption isotherm. The continuity of moisture flux at the interface (see Eq. (5)) is then rewritten as: is the apparent diffusivity of moisture in the porous material.
Since in this section we only compare the simulated results with experimental data which does not include airflow, the airflow velocity v e is set as zero. Without considering the airflow, a 1-D model is enough for coupling moisture transport within two domains. The left boundary condition of the MBR is fixed by the vapor density that is provided by the saturated salt solution. The transport equation in the MBR is solved simultaneously with moisture transport within the material for each iteration.

Fitted mass loss curves
The measured desorption isotherms for three studied cement pastes are fitted by the VG equation (see Eq. (7)). Results of parameters in VG model are provided in Table 2. The decrease of a with the increase of w/c indicates the shifting down of sorption isotherms at high RHs because of more large pores in materials with high w/c. Figure 2 shows that the fitted curves match the measured ones well at the high RHs, while at low RHs (in particular < 40%) the fitting accuracy is slightly low. The large drop in measured S from 40% to 30% RH is caused by the escape of dissolved air in the pore water. With the decreasing RH, the pressure difference between pore water and the ambient environment becomes larger and larger. Until the pressure difference reaches a certain level, the air bubbles grow and then suddenly collapse. [84,85] This process not only can create acoustic events but also results in more water loss from the porous material. Consequently, this sudden mass change cannot be captured by an empirical sorption isotherm equation.
To calibrate the moisture transport model with different boundary conditions, the thickness of the MBR x e in the coupling model must be fixed. Here, we take x e ¼ 10 cm which is roughly same as the distance between the specimen surface to the saturated salt solution in experiments. Results in Figure 3 show that the mass loss curves from the Dirichlet boundary condition are slightly better than the flux boundary condition for some drying cases but in general they are very similar. However, the coupling model is not able to provide accurate curves to the measured ones. The large discrepancy is observed at the early age of drying as the measured curves increase more rapidly than the calculated ones.
The calibrated intrinsic permeability K l values in Table 3 display that the maximum difference in K l about 36% provided by these boundary conditions is found for P06 drying at 33% RH. In fact, this difference is even smaller than that for the same material dried at different external RHs (comparing P05 and P06 dried at 77% and 33% RHs, respectively). The determined K l values for the case of drying at 33% RH are slightly higher than that for 71% RH, that is because the contribution of water vapor is included in the determined K l ; therefore, it was not suggested to determine K l using drying at low RHs. [11]

Validation with measured saturation profiles
The measured saturation profiles are used to validate the compared boundary conditions as the only unknown -intrinsic permeability K l -has already been calibrated by the measured mass loss curves. To compare the simulated and measured saturation profiles, we are able to tell which boundary condition is suitable for the studied cement-based materials. The general trend in Figure 4 is that the calculated profiles inside the material are lower than the measured profiles. However, profiles from the Dirichlet and flux boundary conditions are slightly higher than the coupling model, making them closer to measured ones. If only focusing on the part close to the surface, the coupling model provides the profiles higher than the measured ones, while the other two boundary conditions can show a good agreement with the experimental data.
When a specimen with high water content is moved to a constant RH chamber, the released moisture from the specimen can increase RH in the chamber. The increasing RH is dominated by vapor diffusion, meaning that it takes time to let the released moisture reach the surface of saturated salt solution which is able to absorb the extra moisture and then keep RH constant. This indicates that the actual RH at the specimen surface is always higher than the controlled RH by the saturated salt solution. This is the main reason that we see the saturation profiles of the coupling model is higher at the specimen surface in Figure 4. In fact, the flux boundary condition has the similar behavior as indicated by Eq. (5) that P m, 0 v is always higher than P e v to ensure a drying process. As a result, the saturation profiles of the flux boundary condition are slightly higher than measured ones, but this difference is very small and can be neglected.

Effect of MBR thickness
The above discussion implies that the value of x e has an effect on moisture transport in the porous material. To evaluate this effect, different values of x e are used for simulations (while keeping the same K l that was determined by using the flux boundary condition). One example for P05 paste drying at 71% RH is shown in Figure 5. Clearly, with the increasing x e , the mass loss curve gradually drops and becomes more straight. As indicated in Figure 5b, at the end of drying, RH at the specimen surface is higher for the wide MBR than the narrow ones. With the increase of x e , moisture needs more time to diffuse from the material to the saturated salt solution so that RH at the specimen surface decreases more slowly.
The question is that what the appropriate distance between the specimen and the saturated salt solution is for the coupling model to determine K l . This can be investigated by the relationship between x e and the inversely determined K l . As shown in Figure 6, the value of x e is adjusted between 0 and 20 cm and the corresponding K l is inversely determined by fitting measured mass loss curves. A very good correlation can be seen between x e and K l for all the studied materials as K l increases with x e . The K l values at x e ¼ 20 cm are a few times larger than values at x e ¼ 0 cm, but they are still in the same order of magnitude. The highlighted symbols in Figure 6 are K l determined by the flux boundary condition, which are located in x e % 10 cm, agreeing the value in the previous section. However, if taking K l values determined by the Dirichlet boundary condition, the thickness of the MBR x e becomes zero, which is in    agreement with the definition of the Dirichlet boundary condition. The thickness of MBR x e can in principle be experimentally measured by monitoring RH at different distances from the material surface. The measured x e can be then used to determine K l and further calculate the coefficient E in the flux boundary condition. This will provide a new indirect method to determine K l and E. However, to the best of our knowledge, the measurement of x e is rarely carried out in the literature.
The maximum difference of K l in these boundary conditions is about 36% according to values in Table  3. After reviewing the measured K l in the literature, studies [10,11] found that K l values show a large scatter even for the same cement-based material. Therefore, it was concluded that it is acceptable if the determined K l values hold within the same order of magnitude. [83] Based on this, we would say that three compared boundary conditions show a good agreement in terms of the fitted K l values. If considering the fitting accuracy (see Figure 3), the Dirichlet boundary condition is better than others and is suggested as the first choice for the inverse analysis method. The flux boundary condition is acceptable compared with the coupled modeling. Besides, these two methods are easier to implement in numerical modeling than the coupled modeling. However, this conclusion are drawn based on the laboratory observations, where air in the humidity chamber is stationary.

Effect of airflow on moisture drying kinetics
In the natural conditions, airflow is commonly non-negligible at the surface of concrete structures. It is generally agreed that the evaporation rate is enhanced with the increase of airflow velocity when drying the porous materials, such as soils, [42] gypsum boards, [43,86] woods, [43,44] and fresh concrete. [62,87,88] Nevertheless, this considerable effect is more pronounced for the first stage evaporation (the environment-controlled evaporation) at the low airflow velocity. At high airflow velocity, the evaporation rate becomes less dependent and finally independent on the airflow velocity. [42] The second stage of evaporation (diffusion-controlled evaporation) of course is always independent on the airflow velocity.
To consider the airflow in a numerical model, two boundary conditions, namely, flux boundary condition and the coupling model, can be used. In Eq. (3), the airflow velocity v e can be directly implemented in a 1 D model. However, for the coupling model (see Section 2.1.3), the direction of airflow velocity in 1 D is the same as moisture transport (perpendicular to the porous material surface). The question here is that the perpendicular airflow in the reality is negligible because the porous material can block the airflow and change the direction of airflow. Eventually, only airflow along with the material surface can affect moisture evaporation from the porous material. Therefore, to simulate airflow in both directions (perpendicular to and along with the material surface), a 2 D model is needed. The calculated velocities will be then implemented in the vapor transport model (see Eq. (4)) for both directions.

CFD model
The airflow in the environment is generally simulated by the Navier-Stokes equations. If considering the constant air density, the continuity equation for the case of 2 D in Cartesian coordinates is written as: The conservation of momentum is described by the Cauchy equations, which are written in 2 D as: q a u @v e @x þ v e @v e @y ¼ À @P @y þ g a @ 2 v e @x 2 þ where q a (kg/m 3 ) is the air density, g a (PaÁ s) represents the dynamic viscosity of air, u and v e (m/s) are the airflow velocities in x and y directions, respectively. The pressure P (Pa) is defined as P ¼ p þ p 0 : The isotropic pressure p represents the thermodynamic pressure at the quasi-equilibrium state (a measure of the local squeezing effect due to the fluids flow). The gravitational force g is balanced by the gradient of the reference pressure, rp 0 ¼ g, where the reference pressure p 0 is generally taken from the atmospheric pressure for the airflow. The model setup and boundary conditions for airflow and moisture transport are provided in Figure 7. The porous material on the right-hand side is sealed with only one side being exposed to the environment for the moisture exchange. The thickness of the MBR is set as 0.1 m which is consistent with the above simulations with the coupled model. The moisture transport in the MBR is described by the convectiondiffusion equation. The inlet boundary for the airflow (at the bottom of the MBR) is set as the constant velocity in y direction (zero in x direction). The outlet boundary condition (the upper boundary of the MBR) for the airflow is controlled by the constant pressure (atmospheric pressure). The slip boundary on the lefthand side of the MBR is used for the airflow. On the right-hand side, the Naiver slip condition is used, which creates a discontinuity in the tangential velocity at the material surface, that is, moisture velocity at the surface is higher than that in the porous material. This velocity jump at the surface is reasonable because for the low permeable material the moisture velocity is much smaller than the airflow velocity so that the transition region between the porous material and the environment is extremely thin. In addition, to simplify simulations of airflow, this study considers the airflow in the MBR as the laminar flow.

Simulation results
To analyze the effect of airflow velocity, we take P05 paste drying at 71% RH as an example. Simulation results for two boundary conditions are shown in Figures 8 and 9. If using K l in Table 3 (calibrated with experimental data), the comparison of Figures 8a  and 9a shows negligible effects of airflow velocity on moisture evaporation. The small difference is only observed at the very early age of drying. This indicates that the environment-controlled evaporation stage for the studied material is very short. When the drying starts, it rapidly enters the second stage of evaporation which is known as the diffusion-dominant stage (controlled by moisture transport in the porous material). This is different to the experimental and numerical simulation results in the literature for the most porous materials (e.g., soils, [42] gypsum boards [86] ) which show the evaporation rate increases with the airflow velocity. The main reason of this difference is that the mature cement-based materials are much less permeable than these porous materials studied in the literature.
Evaporation is the competition between the ambient condition and moisture transport in the porous material. If moisture transport in the material is extremely slow, moisture evaporation can only remove moisture in the part adjacent to the surface, while moisture within the material needs longer time to be transported to the surface. For this case, the transport property of the porous material becomes the determinant factor for moisture evaporation. On the contrary, if moisture transport in the material is high, we can expect that moisture in the porous material is able to quickly move the surface and to supply evaporation. To check this case, simulations with higher values of K l (simply Â 10 and 100 of K l in Table 3, representing more permeable porous materials than the studied cement pastes) were performed and results are shown in Figures 8b,c and 9b,c for the cases of flux boundary condition and the coupling model, respectively.
With increasing K l , the duration of the first stage evaporation increases as well (comparing the linear range of curves for three K l in Figures 8 or 9). In each case, the increase of the airflow velocity shortens the first stage evaporation and decreases the transition time between two evaporative stages. However, if the airflow velocity is high than a certain value, the evaporation rate becomes less dependent on it, because of reaching the diffusion-dominant evaporation stage, agreeing with results in the literature. [42] The differences between two boundary conditions for the high K l are obvious in Figures 8 and 9. Figure  8 shows that the mass loss curve gradually increases with the airflow velocity, which is understandable as the airflow basically increases the coefficient E in Eq. (11) and then enhances moisture evaporation at the specimen surface. Nevertheless, Figure 9 indicates that moisture evaporation is largely promoted with a small increase of the airflow velocity, meaning that moisture in the environment can be rapidly removed by the airflow with low velocity. This can be well demonstrated by simulation results in Figure 10 which compares two cases with different airflow velocities (v e ¼0 and v e ¼0.1 m/s). If v e ¼0, the moisture from the porous material can increase RH in the environment (leading to a thick MBR), while if v e is higher, airflow can reduce moisture content in the environment close to the porous material. In other words, the presence of airflow results in a thinner MBR than the case without airflow and thus maintains low RH at the surface of the porous material.
The next question investigated here is whether the consideration of airflow is able to improve the simulation accuracy of drying of cement-based materials. Figure 11 compares results of different airflow velocities for two boundary conditions. By considering the airflow, the model can slightly improve the fitted results, in particular for the initial drying stage, compared with the case when airflow is zero. This is due to the fact that the airflow can keep RH at the material surface as close as to the environmental RH. Thus, the initial moisture content gradient at the surface is much higher than the flux boundary condition and the case without airflow. Note that the coupling model in this study does not consider the boundary layer which we believe is very thin in our experiments. While this is true for concrete with smooth surface, if considering the boundary layer, there will be more room for the model to fit the experimental results.

Implications
In a laboratory measurement, it is generally accepted that a large chamber is better to keep the constant RH. However, if air does not flow, moisture Figure 8. Effects of airflow velocity on mass loss for the flux boundary condition (taking P05 paste drying at 71% as an example). Cases for three different K l are compared. Normal K l represents the case of taking K l in Table 3. Figure 9. Effects of airflow velocity on mass loss for the coupling model (taking P05 paste drying at 71% as an example). Cases for three different K l are compared. Normal K l represents the case of taking K l in Table 3. movement in the chamber is a pure diffusion process, which needs time to let moisture released from the specimen being captured by the saturated salt solution (see Figure 10a-c). This causes RH at the specimen surface higher than the controlled RH. Two ways are able to experimentally maintain RH at the specimen surface close to the controlled RH. The first method is to keep the specimen as close as possible to the saturated salt solution, which is well demonstrated by our simulations. The second way is to let air flow in the chamber. Airflow can largely reduce the thickness of the MBR (see Figure 10d-f). ). This figure only shows the case of the K l ten times higher than the normal K l (corresponding to Figure 9b). Meanwhile, the airflow does not enhance moisture evaporation of the low permeable materials, such as cement-based materials.
Nevertheless, for materials with permeability higher than cement-based materials, it becomes even difficult to maintain RH at the specimen surface close to the controlled RH. Because of the higher permeability, there is always enough moisture to support evaporation at the surface and the air circulation in the chamber is able to accelerate moisture loss. For these materials, the airflow method may not work anymore. The MBR becomes thicker than the low permeable materials and can be significantly affected by moisture transport in the porous materials. Therefore, it is important to couple the moisture transport in the MBR for the high permeable materials, instead of using of simplified boundary conditions (the Dirichlet or the flux boundary condition) in a moisture transport model which may lead to inaccurate results in determined K l .

Conclusion
This paper studies the effects of boundary conditions on moisture transport in cement-based materials regarding the inversely determined liquid water permeability and simulated mass loss kinetics. Three boundary conditions, namely, the Dirichlet (imposed) boundary condition, the flux boundary condition, and the boundary supplied by a coupling model (a convection-diffusion moisture transport model for the ambient environment), are investigated.
Calibration with the lab measured data (without airflow) shows that, the imposed (Dirichlet) boundary condition provides slightly better results in terms of fitting measured mass loss curves and saturation profiles than the other boundary conditions. The largest difference of determined permeability for a given RH is about 36%. Compared with the variations of data reported in the literature, the difference found in this paper is actually not significant. Therefore, we conclude that all studied boundary conditions are acceptable for the low permeable materials (e.g., cementbased materials) to inversely determine liquid water permeability. This may result from the fact that moisture evaporation at the surface of cement-based materials is dominated by the moisture transport in the porous material but not the type of boundary condition.
The flux boundary condition and the coupling model are able to consider the airflow in the environment. They show that the airflow velocity does not affect moisture evaporation for the studied materials. This is due to the fact that the very low permeability of the cement-based materials results in a very short initial drying stage (the environment-controlled evaporation stage). The presence of airflow can largely reduce the thickness of the MBR and keep RH at the surface of the porous material as close to the controlled one as possible. However, for materials with high permeability, the drying mass loss increases with the airflow velocity and the initial drying stage becomes much longer than a low permeable material because the material can supply enough moisture for evaporation. The Dirichlet boundary condition is not able to take into account the complicated moisture interactions between the ambient environment and the material. Therefore, if one wants to improve the simulation accuracy of drying kinetics and inversely determine K l for these materials, a coupling model that is able to simulate moisture transport in the environment becomes necessary.
The present study employed a conventional moisture transport model to simulate moisture state in cementbased materials. As mentioned in the Introduction, to consider the microstructural alteration, the conventional models need to be further developed, such as including the time-dependent permeability. The use of such models to couple with the environmental effect is able to further increase the simulation accuracy. This will be investigated in the future study.

Disclosure statement
No potential competing interest was reported by the authors.