Models for liquid relative permeability of cementitious porous media at elevated temperature: comparisons and discussions

1.   Introduction

Cementitious materials such as concrete subjected to fire loading will experience a substantial pore pressure buildup, which represents the main factor resulting in potential explosive spalling, i.e., the splitting of concrete pieces from the fire-loaded surface, thereby restricting further reinforcement of the fire and greatly jeopardizing the load-carrying capacities of the structures, including tunnels, underground stations and bridges ^[2,3,4]. The buildup of pore pressure in heated concrete is closely related to the moisture transport and phase change processes of water, based on which the different regions in heated concrete are classified. Figure 1 illustrates the three regions in heated concrete: ⅰ) the hot region, ⅱ) the intermediate region and ⅲ) the cold region. Herein, the intermediate region is defined as the region between the depth with the maximum pore pressure and the depth with the maximum saturation. In the hot region, the liquid water vaporizes. In the intermediate region, the vapor condenses into liquid, which is accompanied by a decrease in the gas pressure and an increase in the degree of saturation; this process results in a so-called moisture clog with a very low gas permeability and consequently accelerates the accumulation of pore pressure. In the cold region, the characteristics of the concrete are generally similar to those in the initial state. With an increase in the temperature, the hot and intermediate regions will penetrate to greater depths in the concrete. These phase change processes and the formation of moisture clogs have been proven by experimental as well as numerical investigations ^[5,6,7,8,9]. However, the development and application of a coupled thermo-hydro-chemical multifield framework ^[10,11] is essential to numerically capture the coupled mass-energy transport and phase change processes of water.

The relative permeabilities of concrete with respect to gas $k^{rg}$ and liquid water $k^{rl}$ play important roles in the associated mass-transport processes, which are functions of the saturation degree. Many experimental investigations have been conducted to directly observe the evolution of $k^{rg}$ (for example, see ^{[12,13,14,15]}). In contrast, no reliable experimental tests have been presented to directly obtain $k^{rl}$ due to the difficulties associated with simultaneously testing the liquid permeability of a non-saturated medium and determining its degree of saturation. $k^{rl}$ is commonly determined by semi-analytical models, such as Mualem's model ^[16,17], albeit with questionable reliability. Consequently, to the best of our knowledge, the influence of $k^{rl}$ on the buildup of the pore pressure in heated concrete has not been studied thoroughly. Therefore, the present study was conducted to examine the importance of this effect on the formation of moisture clogs, to estimate $k^{rl}$ and to determine its influence on the pore pressure.

In this work, based on the results and achievements of previous studies, we will show that Mualem's model potentially overestimates the values of $k^{rl}$ for cementitious materials. Considering other models of $k^{rl}$, parametric studies are conducted with a fully coupled thermo-hydro-chemical multifield framework ^[18,19] to study the role of $k^{rl}$ in heated concrete with respect to different initial saturation degrees and fire loadings for high-performance concrete (HPC) as well as normal strength concrete (NSC), thereby indicating that the influences of $k^{rl}$ cannot be neglected. Finally, the water content experimentally obtained through a nuclear magnetic resonance (NMR) test ^[1] is used to disprove Mualem's model for $k^{rl}$, as this model cannot reproduce the steep rises in the moisture content with depth observed in the experiments.

We want to emphasize here that in the model the macro crack (fracture) is not taken into account, which is partly caused by unpredictable initial imperfection and requires sophisticated continuous-discontinuous numerical methods such as ^{[20,21,22,23,24,25,26]}. However, most of the material parameters we used are temperature dependent, based on experimental results. Especially, the temperature-dependent intrinsic permeability is directly taken from published papers which applied well-established experimental set-up, see ^{[27,28,29,30,31]}. In these tests the samples naturally include micro cracks such as interface cracks between cement paste and agammaegate. From this point of view, the influence of thermal induced micro crack on the mass-transport processes of heated concrete is already taken into account.

The remaining parts of this paper are organized as follows. In section 2, Mualem's model for $k^{rl}$ is introduced, and published studies that question the reliability of this model for cementitious materials are presented, indicating the necessity of this work. In section 3, the fully coupled thermo-hydro-chemical framework is presented with the corresponding input parameters. A comparison and discussion of the results are presented in section 3.3, demonstrating different evolutions of the pore pressure buildup, including the reproduction of an NMR test. The paper closes with concluding remarks in Section 4.

2.   Models of $k^{rl}$ for cementitious materials

2.1.   Mualem's model for multiphase flows in concrete

Considering multiphase gas-liquid flow in concrete ^[32,33], the gas and liquid flows denoted respectively as $\mathbf{q}^g$ and $\mathbf{q}^l$ are calculated by the following:

and

where $\mathbf{k}_{int}$ is the intrinsic permeability vector with $\mathbf{k}_{int} = \mbox{k}_{int}\ \mathbf{1}$ for an isotropic material. $\rho^g$ includes the vapor density $\rho^{gw}$ and dry air density $\rho^{ga}$ as $\rho^g = \rho^{gw}+\rho^{ga}$. $k^{rg}$ and $k^{rl}$ are the relative permeabilities of gas and liquid, respectively, as functions of the liquid saturation degree $S_l$. Theoretically, when $S_l$ changes from 0 to 1, $k^{rg}$ will drop from 1 to 0, and $k^{rl}$ will increase from 0 to 1.

Mualem's model ^[16], which is widely used to obtain both $k^{rg}$ and $k^{rl}$ in non-saturated porous media, assumes that the pore configuration is replaced by a pair of capillary elements (i.e., tubes) with $l_1/l_2 = d_1/d_2$ (see Figure 2). In Figure 2, both capillaries obey the Poiseuille equation. When using an equivalent tube with a diameter $D$ and a length $L$ to replace this pair of capillary elements, $D = \sqrt{d_1\ d_2}$ and $D^2\ L = d^2_1\ l_1+d^2_2\ l_2$ will be obtained. As mentioned in ^[16], large pores have an influence that is more important than is generally assumed in this model, since both the length and the cross section of the pore are considered.

Considering van Genuchten's sorption model ^[34,17], the following equations for $k^{rg}$ and $k^{rl}$ can be obtained:

and

where $m$ is a fitting parameter provided in the sorption model such that $m = 0.42$ for desorption and $m = 0.35$ for adsorption (see A).

2.2.   Potential overestimation of $k^{rl}$ by Mualem's model

$k^{rg}$ can be directly tested experimentally. The experimental results for cementitious materials indicate that using Eq 2.3 with the original parameters will greatly overestimate $k^{rg}$ (see Figure 3). Consequently, by fitting the results provided in ^[13,35] with an equation similar to Eq 2.3, we suggest the use of the following:

which is a fitting function mimicking Eq 2.3 that is not model-based but is instead based on experimental results. Zhang et al. provided a similar fitting equation in ^[36].

The overestimation of $k^{rg}$ using Mualem's model raises questions regarding $k^{rl}$. Since Eqs 2.3 and 2.4 are both based on the same model, when Eq 2.3 overestimates the values of $k^{rg}$, it is also possible that Eq 2.4 will overestimate the values of $k^{rl}$.

Mualem's model uses the Poiseuille equation, which is designed for laminar flow and is appropriate for porous media with large pores, such as soils. For cementitious materials, however, the large number of capillary pores renders this assumption unsuitable. In the framework of a macro model, liquids flow along thin films (due to a high capillary pressure) but not through pores directly, resulting in longer and narrower flow paths constituting a type of detour effect (see Figure 4).

Herein, we would like to emphasize that we use a macro model ^[37,36]. Despite their differences at the microscopic level (e.g., differences in pore size distributions), different types of concrete, such as NSC and HPC, are characterized by different macroscopic material properties, including the sorption, permeability, porosity and strength. In this model, all pores are assumed to be connected. Figure 4 presents a simplified schematic showing the influence of capillary forces on the liquid flow in the framework of this macro model. In reality, the capillary pressure cannot exceed the tensile strength of the liquid ^[38]. Accordingly, when $S_l$ drops below some limit value, parts of the liquid will become disconnected, isolated and trapped inside the pores, thereby stopping the flow of the liquid (and gas in the similar way). This indicates the role of the microstructure on $k^{rl}$. In contrast, micromechanical-based models are very commonly used in particle systems involving materials such as sand and powder ^[39,40]. Despite their sophistication, however, such models need precise descriptions or models of the microstructure, thereby necessitating the inclusion of additional variables within the model that are difficult to determine, which is beyond the scope of this paper.

2.3.   Baroghel-Bouny's and Poyet's results

The potential overestimation of $k^{rl}$ with Mualem's model was partly proven experimentally by Baroghel-Bouny in ^[35,14] following tests on the $\mbox{k}_{int}$ values of cementitious materials. As follows, Baroghel-Bouny demonstrated that there are three general methods for experimentally testing $\mbox{k}_{int}$ (as shown in Figure 5:

1. Testing the air flow in a completely dry material ($S_l = 0$);

2. Testing the liquid flow in a completely wet material ($S_l = 1$);

3. Testing the weight change of a wet material within a chamber with a controlled relative humidity (variable $S_l$).

In these methods, methods 1 and 2 avoid the influence of relative permeability, which can be considered as direct testing methods. However, it is actually very hard to make the sample completely saturated at room temperature and normal pressure, questioning the reliability of method 2 ^*. Only method 1 can be considered as a practical and direct testing method^†. Method 3 is an indirect testing method with varying $k^{rg}$ and $k^{rl}$.

^*We still show the results obtained with method 2 for comparisons.

^†In method 1, firstly the samples are slowly (heating rate smaller than 1 $^\circ$C/min) heated from room temperature to around 105 $^\circ$C and kept in the heating chamber for several hours, assuring the samples are completely dry without micro cracks. Then the samples are slowly cooled and tested with dry air or nitrogen.

Regarding method 3, when using Eqs 2.3 and 2.4 with $m$ values ranging from 0.4191 to 0.4854, Baroghel-Bouny ^[14,35] back-calculated the intrinsic permeabilities of cementitious materials with different water-cement ratios; the resulting permeabilities were much smaller than those obtained by method 1 and method 2 for all types of cementitious materials tested, see Table 1. Similar results were also obtained in other papers; $\mbox{k}_{int, 1}$ is approximately in the range of $10^{-16}\sim 10^{-18}\mbox{m}^2$ ^{[27,28,29,30,31]} when $\mbox{k}_{int, 3}$ is in the range of $10^{-20}\sim 10^{-23}\mbox{m}^2$ ^{[41,42,43,44]}. This fact is astonishing when $\mbox{k}_{int}$ is intrinsic, as its values should primarily depend on the microscopic structure of the porous medium and not the physical properties of the flow.

Furthermore, in the drying test (method 3), Mainguy et al. ^[45] proved that the flow of gas (vapor) has a minor influence on the mass transport of water inside a material compared with the flow of liquid. Hence, the small value of $\mbox{k}_{int, 3}$ can be attributed to the overestimation of $k^{rl}$ with Eq 2.4. Poyet ^[42] back-calculated $k^{rl}$ with the results obtained using method 3. Based on the results of the drying test on a concrete with WCR = 0.47, Poyet used Pickett's equation ^[46] in addition to linear and polynomial functions for sorption curves; the obtained values of $k^{rl}$ were smaller than those obtained with Mualem's model (see Figure 6). In summary, the differences between $\mbox{k}_{int, 1}$ and $\mbox{k}_{int, 3}$ are potentially caused by overestimated values of $k^{rl}$ obtained via Mualem's model. A similar conclusion was also drawn by Waldeh and Perrin ^[47].

2.4.   Cases of $k^{rl}$ for numerical studies

The combination of Eq 2.4 with $\mbox{k}_{int, 3}$ has been widely used as a set of input parameters for performing durability analyses on concrete structures with environmental temperatures of less than 100 ℃, in which the overestimation of $k^{rl}$ has already been mitigated by small values of $\mbox{k}_{int, 3}$. However, method 3 cannot be used to ascertain the evolution of $\mbox{k}_{int, 3}$ (even the residual values) with elevated temperatures, because wetting a dehydrated cementitious sample will cause rehydration. When using $\mbox{k}_{int, 3}$ to analyze concrete specimens subjected to fire loading, extra assumptions must be introduced to account for the evolution of $\mbox{k}_{int, 3}$ with the temperature (for example, see ^[48,49]). To simulate heated concrete, $\mbox{k}_{int, 1}$ is more widely used as the input parameters than $\mbox{k}_{int, 3}$ ^[50,51,52] because the overestimation of $k^{rl}$ with Mualem's model must be considered when using the latter.

Other models and curves exist in addition to Mualem's model for estimating the values of $k^{rl}$. Ranaivomanana et al. ^[53] presented a pore network model to obtain $k^{rl}$ numerically in consideration of the achievements of preceding research, such as that conducted by Ishida et al. ^[54]. The authors assumed that three modes of pores exist in cementitious materials that can be categorized according to their sizes and that the cross-sectional shape of the pore is an oblong (see Figures 1 and 3 of ^[53]). The following relationship was given:

where $r_{cr}$ is a calibration parameter. Comparisons with the experimental results of the isotherm drying tests in ^[55,41] prove the reliability of this model insomuch that the obtained $\mbox{k}_{int, 3}$ values are much more reasonable than those of $\mbox{k}_{int, 1}$. With this model, Ranaivomanana et al. ^[56] numerically obtained the $k^{rl}$ value of an NSC with WCR = 0.4872 and that of a cement paste with WCR = 0.2667. Based on their results, we fit the data by mimicking Eq 2.4 and obtain the following fitting equation:

Furthermore, based on the experimental results of drying tests, Chaparro et al. ^[57] derived the following empirical equation from back-calculations:

where $k^{rl}\left(S_l = 1\right) = 0.01\neq 1$.

Herein, we consider three cases of $k^{rl}$ for the numerical simulations in the next section as follows (see Figure 7):

1. Case 1: using Eq 2.4 with $m = 0.35$ as Mualem's model-based function;

2. Case 2: using Eq 2.7 as a fitting function from the model-based numerical results in ^[56];

3. Case 3: using Eq 2.8 as an empirical function.

In all cases, the values of $k^{rl}$ drop below $10^{-3}$ when $S_l < 0.4$.

3.   Coupled thermo-hydro-chemical analysis

3.1.   Control equations

For the numerical simulations, we use the fully coupled thermo-hydro-chemical (THC) model presented in ^[58,59] and applied previously by ^[60,61,36]. Detailed discussions regarding this model and the other numerical models are given in ^[62,63]. The following governing field equations are constructed by obeying mass and energy conservation with $p^c$, $p^g$ and $T$ as the unknowns:

Mass-balance equation for the water phase:

Mass-balance equation for the dry-air phase:

Energy-balance equation:

where $\dot{m}_{vap}$ is determined by

In Eqs 3.1 through 3.4, the symbol $D^s\ / \ Dt$ denotes the time derivative of the corresponding parameter relative to the solid skeleton. The nonlinear equations are discretized by finite element methods in the space domain and the finite difference method in the time domain, where $D^s\ / \ Dt\approx\Delta\ / \ \Delta t$. The obtained global stiffness matrix is non-symmetric, Newton's method is used to solve the equations with small time steps for enhancing numerical stability. The program is coded in hybrid Fortran-C language with Gmsh for pre-processing and Gnuplot for post-processing. Interested readers are referred to ^[64,65,66].

3.2.   Input parameters

Regarding the size of normal structures potentially suffered from fire loading such as tunnel segment and walls of building. A model with depth of 50 cm is considered here. One side of the model is prescribed with a heated boundary, and the other side is isolated. The following equation is used for the heat flux through the heated boundary $q_T$ [W m$^{-2}$]:

where $\beta_T = 25$ W m$^{-2}$ K$^{-1}$ is the temperature convection coefficient, $e = 0.56$ is the emissivity of concrete, and $s = 5.67\cdot10^{-8}$ W m$^{-2}$ K$^{-4}$ is the Stefan-Boltzmann constant. Two types of heat sources $T_\infty$ are considered: RABT and ISO 834 (see Figure 8). The fire temperature of the RABT source increases linearly from 20 ℃ at 0 s to 1200 ℃ at 300 s. Although the RABT fire temperature will drop linearly from 1200 ℃ at 3600 s to 20 ℃ at 10200 s, the maximum pore pressure will be reached within 3600 s. Hence, we consider only the parts of the RABT temperature curve that remain, including the increasing components of the curve. The capillary boundary is taken as suggested in ^[63] while assuming a boundary vapor pressure of $p^{gw}_\infty = 1020$ Pa and a vapor convection coefficient of $\beta_v = 0.025$ m s$^{-1}$. The gas pressure of the boundary is kept constant at $p^{g}_\infty = 101325$ Pa.

We consider one type of HPC with WCR = 0.34 ^[28] and one type of NSC with WCR = 0.43 ^[27], the intrinsic permeabilities of which are shown in Figure 9. Considering the work presented in ^[63,36,30], we adopt and list the main input parameters in Table 2 with regard to the following empirical equation for the evolution of $n$ with $T$:

where $n_0$ is the initial porosity and $T_0$ is the initial temperature.

3.3.   Parameter studies and discussion

3.3.1.   Maximum pore pressure built-up

We first analyze the evolution of the maximum pore pressure $p^g_{max}$ in heated concrete. The time-$p^g_{max}$ curves of the RABT fire loading and different values of $S_{l, 0}$ are shown in Figure 10. The results indicate that the historical maximum value of $p^g_{max}$ will be reached within 400 s after fire loading. Once the historical maximum value of $p^g_{max}$ is reached, $p^g_{max}$ will change only slightly, as was also reported in ^[67,18,49]. The historical maximum value of $p^g_{max}$ for the HPC is almost twice that for the NSC, see for example Figure 10(c) and (d), indicating that the intrinsic permeability has a substantial effect on the pore pressure buildup, which was summarized as a "log-linear influence" in ^[36].

The values of $p^g_{max}$ obtained in $k^{rl}$ Cases 2 and 3 are higher than those obtained in $k^{rl}$ Case 1, in which $p^g_{max}$ also begins to rise earlier (see the circle in Figure 10(d)), suggesting the faster formation of moisture clogs in Cases 2 and 3. We would like to emphasize here that although the differences in $k^{rl}$ are not tremendous when $S_l$ = 0.5 (see Figure 7), a difference in $p^g_{max}$ of approximately 20% can be obtained when $S_{l, 0}$ = 0.5, see Figure 10(a) and (b). This is partly induced by the high gradient of $p^c$ in non-saturated concrete. Similar characteristics can also be found in the results with the ISO 834 fire loading, in which $p^g_{max}$ becomes almost constant after 600 s (see Figure 11).

3.3.2.   Distributions of $p^g$, $T$ and $S_l$

The general evolutionary trends of $p^g_{max}$ obtained with the RABT and ISO 834 fire loadings are similar; the value of $p^g_{max}$ will increase until a historical maximum value is reached, after which it will change only slightly. Hence, we consider only the results of the RABT fire loading in this section.

The distributions of $p^g$ and $T$ at 200 s and 300 s shown in Figure 12 indicate the following:

● The temperature distributions are very similar. For the same type of concrete, higher values of $S_{l, 0}$ result in slightly lower temperatures when more energy is consumed by a greater initial liquid water content. See the arrow marked in Figure 12 (temperature results at 300 s). Similar results can be found in ^[67];

● For concrete with higher $S_{l, 0}$ values, the influence of $k^{rl}$ on the pore pressure is more conspicuous; see Figure 12(c) and (d);

● For the NSC with higher $S_{l, 0}$ values, the influence of $k^{rl}$ on the pore pressure is considerable, especially at the early stage of fire loading; see the maximum pore pressure regarding Cases 1 and 2 at 200 s in Figure 12(d);

● The temperatures at the depths with $p^g_{max}$ are mostly in the range of 170 ℃~0 ℃, which was also reported in ^[36].

To investigate the formation of moisture clogs, the distributions of $S_l$ and $T$ are shown in Figure 13, indicating the following:

● $S_l$ begins to increase at the corresponding depth where $p^g$ begins to decrease;

● The $S_l$ curves for Case 1 are different for the HPC and NSC insomuch that the $S_l$ curves rise faster with depth in the HPC than in the NSC (see the arrows marked in Figure 13);

● For the HPC and NSC, the $S_l$ curves for Cases 2 and 3 are generally similar when the $S_l$ curves for Cases 2 are even steeper than those for Cases 3. Meanwhile, in Cases 2 and 3, the $S_l$ curves increase rapidly with increasing depth (i.e., a steep increase in $S_l$), and moisture clogs are formed in a narrow region;

● The different shapes of $S_l$ are primarily caused by the values of $k^{rl}$ at corresponding $S_l$. When the value of $k^{rl}$ is lower (see Figure 7(b)), it becomes more difficult to achieve further moisture transport once vapor condenses into liquid water, accelerating the formation of moisture clogs.

3.4.   Comparison with the NMR results from van der Heijden et al ^[1]

van der Heijden et al. ^[1] used an NMR system for the non-destructive real-time monitoring of moisture transport in heated concrete, in which the concrete samples were heated by an array of four 100-W halogen lamps. A maximum temperature of approximately 400 ℃ was reached in the specimens.

The test specimens in this study are concrete samples (WCR = 0.5) with different initial relative humidity $\mbox{RH}_0$. Two cases, namely, $\mbox{RH}_0 = 0.97$ and $\mbox{RH}_0 = 0.75$, are considered in our numerical simulation. All of the material parameters with values that are different from those in Table 2 are listed in Table 3. The evolution of the intrinsic permeability is assumed to be the same as that of the NSC shown in Figure 9.

Based on the monitored $S_l$, we assume a bubble pressure $p^c_b = 2.2\ \mbox{MPa}$. Based on the monitored temperature profile, we assume a heating boundary with $\beta_T = 0$ W m$^{-2}$ K$^{-1}$, $e = 0.56$, $T_\infty = 550\; ^\circ\mbox{C}$ for $\mbox{RH}_0 = 0.75$ and $T_\infty = 520\; ^\circ\mbox{C}$ for $\mbox{RH}_0 = 0.75$; see 14 for the distributions and evolutions of the numerically and experimentally obtained temperature profiles, which indicate the same rules provided in the last subsection insomuch that the temperatures from the different cases of $k^{rl}$ are very similar.

In ^[1], the chemically bound water was first determined by calibrating the NMR results. In the coupled THC model we applied, chemically bound water is considered as a "source" of free water in the equations, which is not accounted in the water volume content $\left(n\ S_l \right)$.

Regarding moisture profiles, for $\mbox{RH}_0 = 0.75$, the distributions and evolutions of the water volume content $\left(n\ S_l \right)$ are shown in Figure 15. Because of the low values of $S_{l, 0}$, $S_{l}$ does not have great potential for increasing, and the values of $k^{rl}$ are generally small for all cases. The values of $\left(n\ S_l \right)$ obtained with different cases of $k^{rl}$ do not exhibit substantial differences, whereas the peak values of $\left(n\ S_l \right)$ obtained in Cases 2 and 3 are slightly higher than those obtained in Case 1. However, compared with the experimental results, the peak values of numerically obtained $\left(n\ S_l \right)$ are much lower. This phenomenon is currently not well understood; we believe it is caused by the assumed capillary models. In cementitious materials, in addition to capillary pores, there exist some big pores with $p^c\approx0$ that provide channels for liquid water transport. However, these pores are not fully considered in the present macro model, indicating the necessity of developing a multiscale pore model for better capturing mass transport in heated concrete. Moreover, in the numerically obtained results, there is increasing moisture content in the shallow depth (for an example, see the arrow in Figure 15(b)). This increase is not caused by chemically bound water but rather by the sorption model; see A. When the temperature is below $T_{crit}$, $S_l$ will not become zero because the corresponding value of $p^c$ is infinity; see Eq A.1.

For $\mbox{RH}_0 = 0.97$, the distributions and evolutions of $\left(n\ S_l \right)$ are shown in Figure 16. Although none of the cases of $k^{rl}$ capture the high peak values of $\left(n\ S_l \right)$ observed experimentally, Cases 2 and 3 capture the steep increase in $S_l$; see the arrow in Figure 16(a), (c) and (d). In contrast, Case 1 of $k^{rl}$ is not capable of reproducing this steep increasing trend and overestimates the depth of the drying zone (see the circle in Figure 16(b)). From this perspective, Cases 2 and 3 provide better results than Case 1.

Finally we summarize the advantages and disadvantages of the models used in Cases 1 to 3 in Table 4. As mentioned before, the values of maximum pore pressure obtained when using Cases 2 and 3 are around 15% to 20% higher than the values obtained when using Case 1, while the temperature results are generally similar. On the other hand, none of these models completely reproduces the moisture increasing observed in NMR results, indicating the necessities for developing improved models.

4.   Conclusions

In this paper, we reviewed the relative permeabilities of cementitious materials both to air $k^{rg}$ and to liquid $k^{rl}$, and we investigated the suitability of Mualem's model for estimating $k^{rl}$. By showing the experimental and numerical results from other researchers, we indicate the possibility that Mualem's model potentially overestimates not only $k^{rg}$ but also $k^{rl}$, through which we derive the motivation of this investigation. By using a fully coupled thermo-hydro-chemical multifield model, we study the influence of $k^{rl}$ on the pore pressure buildup in heated concrete while considering Mualem's model as well as other models given in the literature. Furthermore, through a comparison with NMR experimental results, we prove that certain other models produce results that are more agreeable than those of Mualem's model, which should be considered by researchers for spalling analyses in the future.

Acknowledgments

The authors gratefully acknowledge financial support by the NSFC (41572260), (51208378), and (51809069) and the Fundamental Research Funds for the Central Universities.

Conflict of interest

The authors declare no conflict of interest in this paper.

Appendix

A.   Sorption model

In this paper, the sorption model is the same as the model used in ^[36], which was modified from the van-Genuchten model ^[34,55] and presented in ^[62]. A recent modified version was presented in ^[9], and the experimental investigations are found in ^[68]. The relationship between the capillary pressure $p^c$ and water saturation degree $S_l$ is described by

where

and

In Eqs (A.2) and (A.3), $B_s = 30$ MPa, $N_s = 1.2$, $Z_s = 0.5$ ℃, $T_{crit} = 374.15$ ℃, $T_0 = 23$ ℃, $E_{s, 0} = 2605.11$, and $T_b = 100$ ℃ ^[63]. In addition to this model, other models, such as that of ^[49], are also available for investigating elevated temperatures.

The parameters $p^c_b$ and $m$ in Eq. (A.1) are determined based on sorption isotherms obtained at $T\leq T_b$ using the Kelvin-Laplace equation represented in Eq. (A.4). By fitting the experimental results given in ^[14] (see Fig. 17) regarding both adsorption and desorption processes, the values of $p^c_b$ and $m$ are obtained and shown in Table 5 for different water-cement ratios (WCRs). With Eq. A.1, the relationship between the relative humidity ($\mbox{RH}$) and $S_l$ at high temperatures can also be obtained (see Figs. 18 and 19). Although the physical meaning of the capillary pressure at a high temperature is slightly different ^[69], the results indicate that the influence of $p^c_b$ on the sorption isotherms at high temperatures becomes negligible. Moreover, the differences between the desorption and adsorption curves become smaller with increasing temperature.

Accounting for both the desorption curve and the adsorption curve in the same numerical model when they correspond to different values of $S_l$ at the same $\mbox{RH}$ is a complicated task. In heated concrete, since the desorption and adsorption curves become similar at high temperatures and the adsorption process is dominant at lower temperatures, we use the adsorption curve (right column in Table 5) throughout this paper.