Understanding the mass, momentum, and energy transfer in the frozen soil with three levels of model complexities

Abstract. Frozen ground covers a vast area of the Earth's surface and it has important ecohydrological implications for cold regions under changing climate. However, it is challenging to characterize the simultaneous transfer of mass and energy in frozen soils. Within the modeling framework of Simultaneous Transfer of Mass, Momentum, and Energy in Unsaturated Soil (STEMMUS), the complexity of the soil heat and mass transfer model varies from the basic coupled model (termed BCM) to the advanced coupled heat and mass transfer model (ACM), and, furthermore, to the explicit consideration of airflow (ACM–AIR). The impact of different model complexities on understanding the mass, momentum, and energy transfer in frozen soil was investigated. The model performance in simulating water and heat transfer and surface latent heat flux was evaluated over a typical Tibetan plateau meadow site. Results indicate that the ACM considerably improved the simulation of soil moisture, temperature, and latent heat flux. The analysis of the heat budget reveals that the improvement of soil temperature simulations by ACM is attributed to its physical consideration of vapor flow and the thermal effect on water flow, with the former mainly functioning above the evaporative front and the latter dominating below the evaporative front. The contribution of airflow-induced water and heat transport (driven by the air pressure gradient) to the total mass and energy fluxes is negligible. Nevertheless, given the explicit consideration of airflow, vapor flow and its effects on heat transfer were enhanced during the freezing–thawing transition period.



Introduction
Frozen soils, have been reported with significant changes under climate change (Cheng and Wu, 2007;Hinzman et al., 2013;Zhao et al., 2019). Changes in freezing/thawing process can alter soil hydrothermal regimes, activate/close the water flow pathways and vegetation development (Walvoord and Kurylyk, 2016). Such changes will further considerably affect the spatial pattern, the seasonal to interannual 30 variability and long term trends in land surface water, energy and carbon budgets and then the land surface atmosphere interactions (Schuur et al., 2015;Subin et al., 2013;Walvoord and Kurylyk, 2016). Understanding the soil freeze/thaw processes appears to be the necessary path to the better water resources management and ecosystem protection in cold regions.
When soil experiences the freeze/thaw process, there is a dynamic thermal equilibrium system of ice, liquid 35 water, water vapor and dry air in soil pores. Water and heat flow are tightly coupled in frozen soils. Coupled water and heat physics, describing the concurrent flow of liquid, vapor as well as heat flow, was first proposed by Philip and De Vries (1957), (hereafter termed as PdV57) considering the enhanced vapor transport. The PdV57 theory has been widely applied for detailed understanding of soil evaporation during the drying process (De Vries, 1958;De Vries, 1987;Milly, 1982;Novak, 2010;Saito et al., 2006). The attempts to 40 simulate the coupled water and heat transport in frozen soils started in 1970s (e.g., Guymon and Luthin, 1974;Harlan, 1973). Since then, numerical tools able to and subjected to simulate one dimensional frozen soil were increasingly developed. Flerchinger and Saxton (1989) developed the SHAW model with the capacity of simulating the coupled water and heat transport process. Hansson et al. (2004) accounted for the phase changes in HYDRUS-1D model and verified its numerical stability with rapidly changing boundary 45 conditions. Considering the two components (water and gas) and three water phases (liquid, vapor, and solid), Painter (2011) developed a fully coupled water and heat transport model MarsFlo. These works together with other modifications, simplifications, generate a series of hierarchy of frozen soil models (detail reviewed by Li et al., (2010) and Kurylyk and Watanabe (2013)).
Air flow has been reported important to the soil water and heat transfer process under certain conditions 50 (Prunty and Bell, 2007;Touma and Vauclin, 1986). Zeng et al., (2011a, b) found that soil evaporation is enhanced after precipitation events by considering air flow and demonstrated that the air pressure induced advective fluxes inject the moisture into the surface soil layers and increase the hydraulic conductivity at top layer. The diurnal variations of air pressure resulted in the vapor circulation between the atmosphere and land surface. Wicky and Hauck (2017) reported that the temperature difference between the upper and the lower 55 part of a permafrost talus slope was significant and attributed it to the airflow induced convective heat flux. Yu et al., (2018) analyzed the spatial and temporal dynamics of air pressure induced fluxes and found an interactive effect as the presence of soil ice. The abovementioned studies demonstrate that the explicit consideration of air flow has the potential to affect the soil hydrothermal regime. https://doi.org/10.5194/hess-2020-253 Preprint. Discussion started: 15 June 2020 c Author(s) 2020. CC BY 4.0 License. Current land surface models (hereafter LSMs), however, usually adopted a simplified frozen soil physics 60 with relative coarse vertical discretization (Koren et al., 1999;Niu et al., 2011;Swenson et al., 2012;Viterbo et al., 1999). In their parameterizations, soil water and heat interactions can only be indirectly activated by the phase change processes, the mutual dependence of liquid water, water vapor, ice and dry air in soil pores is of course absent. This mostly lead to unrealistic physical interpretations and worse performance regarding to the hydrothermal, ecohydrological dynamics (Cuntz and Haverd, 2018;Novak, 2010;Wang and Yang, 65 2018). Specifically, Su et al. (2013)  Is it necessary to consider such a fully physical mechanism in LSMs? These two questions frame the scope 70 of this work.
In this paper, we incorporated the various complexity of water and heat transport mechanisms into a common modeling framework (STEMMUS-FT, Simultaneous Transfer of Energy, Momentum and Mass in Unsaturated Soils with Freeze-Thaw). With the aid of in situ measurements collected from a typical Tibetan meadow site, the pros and cons of different model complexities were investigated. Subsurface energy budgets 75 and latent heat flux density analyses were further conducted to illustrate the underlying mechanisms considering different coupled water-heat physics. Section 2 describes the experimental site and the implementation of increasing complexity of subsurface physics into STEMMUS framework. Performance of different models is presented in Section 3 together with the subsurface heat budgets and latent heat flux density analyses. Section 4 discusses the effects of considering coupled water-heat transport and air flow in 80 frozen soils. Conclusion is made in Section 5.

Experimental site
Maqu station, equipped with a catchment scale soil moisture and soil temperature (SMST) monitoring network and micro-meteorological observing system (Dente et al., 2012;Su et al., 2011;Zeng et al., 2016), is 85 situated on the north-eastern fringe of the Tibetan Plateau (33°30'-34°15'N, 101°38'-102°45'E). According to the updated Köppen-Geiger climate classification system, it can be characterized as a cold climate with dry winter and warm summer (Dwb). The mean annual air temperature is 1.2 ℃, and the mean air temperatures of the coldest month (January) and warmest month (July) are about -10.0 ℃ and 11.7 ℃, respectively. Alpine meadows (e.g., Cyperaceae and Gramineae), with heights varying from 5 cm to 15 cm 90 throughout the growing season, are the dominant land cover in this region. The sandy loam and silt loam are found by in situ soil sampling and organic soil with a maximum of 18.3 % organic matter for the upper soil layers (Dente et al., 2012;Zhao et al., 2018;Zheng et al., 2015).

Mass and energy transport in unsaturated soils
On the basis of STEMMUS modelling framework, the increasing complexity of vadose zone physics in frozen soils was implemented as three alternative models (Table 1) To accommodate the specific conditions of a Tibetan meadow, the total depth of soil column was set as 1.6 m. The vertical soil discretization was designed finer in the upper soil layers (0.1-2.5cm for 0-40cm) than that in the lower soil layers (5-20cm for 40-160cm). Three aforementioned models adopted the same set of soil parameters, shown as Table 2.

Soil hydrothermal profile simulations
The performance of model with various soil physics in simulating the soil thermal profile information is illustrated in Fig. 1. Both CPLD and CPLD-AIR model well reproduced the time series of soil temperature at different soil depth except for the 40cm, which is probably due to the inappropriate measurements (e.g., improper placement of sensors). However, there are significant discrepancies of soil temperature simulated 140 by the unCPLD model. Compared to the observations, a stronger diurnal behavior of soil temperature in response to the fluctuating atmospheric forcing was found and the earlier stepping in/stepping out of the frozen period was reproduced by the unCPLD model. Such differences enlarged at deeper soil layers with large BIAS and RMSE values (Table 3).

145
During the rapid freezing period, a noticeable overestimation of diurnal fluctuations and early and fast decreasing of soil liquid water content was simulated by unCPLD model. Moreover, the strong diurnal fluctuations and early increase of liquid water content were also found during the thawing period. The early thawing of soil water even lead to an unrealistic refreezing process at 80 cm (from 88 th to 92 th day after December 2015), which is due to the simulated early warming of soil by unCPLD model (Fig. 1). Such 150 discrepancies were significantly ameliorated from CPLD and CPLD-AIR simulations. Nevertheless, all three models can well capture the diurnal variations and magnitude of liquid water content during the frozen period.
Note that there is an observable difference between CPLD and CPLD-AIR simulated soil liquid water content at shallower soil layers during the thawing process (e.g., Fig. 2, 5cm).

155
The time series of freezing front propagation derived from the measured and model simulated soil temperature was reproduced as Figure 3. Initialized from the soil surface, the freezing front quickly develops downwards till the maximum freezing depth. The thawing process starts from both the top and bottom, mainly driven by the atmospheric heat and geothermal heat source, respectively. Such characteristics were well https://doi.org/10.5194/hess-2020-253 Preprint. Discussion started: 15 June 2020 c Author(s) 2020. CC BY 4.0 License. captured by both the CPLD and CPLD-AIR model in terms of freezing rate, maximum freezing depth and 160 surface thawing process. While the unCPLD model tended to present a more fluctuated and rapid freezing front propagation and a deeper maximum freezing depth which is early reached. The effect of atmospheric heat source on soil temperature was overestimated by the unCPLD model as shown by the stronger diurnal early onset of the thawing process.

165
The performance of model with different soil physics in reproducing the latent heat flux dynamics is shown in Fig. 4. Compared to the observed LE, there is a significant overestimation of half-hourly latent heat flux, which significantly degraded the overall performance using unCPLD model. The occurrence of such overestimation was notably reduced using CPLD and CPLD-AIR model. While the general underestimation of latent heat flux by the CPLD and CPLD-AIR model was found mostly during the freezing-thawing 170 transition period (Fig. 5b), when the soil hydrothermal states are not well captured ( Fig. 1 &2).
The overestimation of surface evapotranspiration by unCPLD model was significant during the initial freezing and freezing-thawing transition period (Fig. 5a, December & February). During the rapid freezing period (January), unCPLD model presented a good match in the diurnal variation compared to the observations. The monthly average diurnal variations were found to be well captured by CPLD and CPLD-

175
AIR models. Figure 5b shows the comparison of observed and model simulated cumulative surface evapotranspiration. The overall overestimation of surface evapotranspiration by unCPLD model can be clearly seen in Fig. 5b. Days at the initial freezing periods, with high liquid water content simulations, accounted for more than 90% of the overestimation. The initial stage overestimation of surface evapotranspiration was significantly reduced by CPLD and CPLD-AIR simulations. Slight underestimation 180 of cumulative surface evapotranspiration was simulated by CPLD and CPLD-AIR model with values of 3.98% and 4.78%, respectively. The dynamics of heat balance components at 5 cm soil layer was tested for the freezing-thawing transition 205 period (Fig. 6 d,

210
means larger amount of water vapor was evaporated from 5 cm soil layer (with more negative LHF term) from CPLD-AIR simulations than that from CPLD simulations, which explains the lower liquid water content for CPLD-AIR model (Fig. 2, 5 cm).

Subsurface latent heat flux density
To give more context to the results, the spatial and temporal distribution of model simulated latent heat flux 215 density (Sh), − / , during the freezing and freezing-thawing transition period was shown as Fig. 7.
For the unCPLD model, the latent heat flux density (Sh) is not available due to its inability to depict the vapor flow process. Figure 7a shows that there is a strong diurnal variation of Sh at upper 0.1cm soil layers. Such diurnal behavior along the soil profile was interrupted by soil layer of 1cm, at which the water vapor consistently moved 220 upwards as evaporation source (termed as evaporative front). The path of this upward water vapor ended at soil depth of 20cm from the 6 th of December, where the freezing front developed. Compared to the upper 0.1cm soil, a weaker diurnal fluctuations of Sh was found at lower soil layers. For CPLD-AIR model, the vapor transfer patterns are similar to that of CPLD model (Fig. 6b). There were isolated connections of condensed water vapor between upper 1cm soil and the lower soil layers (Sh>0, e.g., 6 th , 7 th , 9 th , and 10 th of 225 December), possibly associated with the downward air flow (see Fig. 12 in Yu et al. (2018)). The large difference in magnitude of latent heat flux density between CPLD and CPLD-AIR model appeared mainly isolated at upper soil layers (Fig. 7c). At soil layers between 1cm and 20cm, CPLD-AIR model simulated less in condensation vapor area (Sh>0) and more in the evaporation area (Sh<0), indicating that CPLD-AIR https://doi.org/10.5194/hess-2020-253 Preprint. Discussion started: 15 June 2020 c Author(s) 2020. CC BY 4.0 License.
model produced an additional amount of condensation and evaporation water vapor compared with CPLD 230 model (Fig. 7c).
Similar to that during the freezing period, strong diurnal variations at upper soil layers, interruption of diurnal patterns by the constant upward evaporation of intermediate soil layers, and weak diurnal variation at lower soil layers of Sh can be clearly observed along soil profile during the freezing-thawing transition period (Fig.   7d, e). While the maximum evaporation rate was less than that during the freezing period. The consistent 235 evaporation zone developed to a depth of 5 cm. The path for the upwards water vapor tended to develop deeper than 30cm with the absence of soil ice. The simulation by CPLD-AIR model produced more condensation and less evaporation water vapor than that by CPLD model can be seen more clearly (Fig. 7f).
In addition, steadily more evaporation water vapor from soil depth of 5 cm was simulated by CPLD-AIR model compared to CPLD model. This confirms the aforementioned point that during the freezing-thawing 240 transition period, large LHF values were simulated by CPLD-AIR model (Fig. 6).

Coupled Water and Heat Transfer Processes
The coupled water and heat transfer is realized via considering the vapor transfer processes. The mutual dependence of soil water, in different phases (liquid, water vapor, and ice), and heat transport is enabled to 245 facilitate our better understanding of the complex soil physical processes (e.g., Fig. 6-7). Specifically, the coevolution of soil moisture and soil temperature (SMST) profiles simulated by CPLD model was closer to the observation than that by unCPLD model. In addition, significant enhancement in portraying the monthly average diurnal variations of surface evapotranspiration and cumulative evapotranspiration can be found from CPLD model simulations, which constraints the hydrothermal regimes especially during the freezing-250 thawing transition periods (Fig. 1, 2& 5). During the freezing period, liquid water in the soil freezes, which is analog to the soil drying process, and water vapor fluxes instead of liquid fluxes dominate the mass transfer process. Neglecting such important water flux component unavoidably results in unrealistic simulations of surface evapotranspiration and SMST profiles. From the energy budget perspective, the contribution of vapor fluxes to the heat balance budget is more evidenced at soil layers above the evaporative front than that below 255 it (e.g., Fig. 6a vs. Fig. 6d, corresponding evaporative front shown as Fig. 7a vs. Fig. 7d). The downward latent heat flux from CPLD model makes the subsurface soil warmer, which reduces the temperature gradient ( / ). This further results in the weaker diurnal fluctuations of conduction term for CPLD model than that for the unCPLD model (Fig. 6). At the soil layers below the evaporative front, the heat flux source from vapor diffusion process (LHF) is negligible. Thermal retard effect as the presence of soil ice, expressed as 260 the apparent heat capacity term (Capp), dominates the heat transfer process in frozen soils. CPLD model, considering the thermal effect on water flow, usually has a larger water capacity value / than unCPLD model. As such, the intense thermal impedance effect leads to the results that CPLD model has a weaker diurnal fluctuation of soil temperature than unCPLD model at subsurface soil layers.

265
Since soil pores are filled with liquid water, vapor and dry air, taking dry air as an independent state variable can facilitate better understanding of the relative contribution of each component in soil pores to the mass and heat transfer in soils. The results show that the dry air-induced water and heat flow is negligible to the total mass and energy transfer (Yu et al., 2018;Zeng et al., 2011a). Nevertheless, dry air can affect soil hydrothermal regimes significantly under certain circumstances. Wicky and Hauck (2017) reported that the 270 airflow-induced convective heat transfer resulted in a considerable temperature difference between the upper and lower part of a permafrost talus slope and thus have a remarkable effect on the thermal regime of the talus slope. Zeng et al. (2011a) demonstrated the airflow-induced surface evaporation enhanced after precipitation events, since the hydraulic conductivity of topsoil layers increased tremendously due to the airflow from the atmosphere into the soil. In this study, we found that the explicit consideration of airflow 275 have the model produce an additional amount of subsurface condensation and evaporative water vapor in the condensation region and evaporation region, respectively ( Fig. 7c & f). The effect of latent heat flux on heat transfer was enhanced by airflow during the freezing-thawing transition period (Fig. 6), which further affects the subsurface hydrothermal simulations (e.g., Fig. 2).

280
On the basis of STMMUS modeling framework with various complexity of water and heat transfer physics (unCPLD, CPLD and CPLD-AIR model), the performance of each model in simulating water and heat transfer and surface evapotranspiration was tested on a typical Tibetan meadow. Results indicate that compared to the in situ observations, the unCPLD model tended to present an earlier freezing and thawing date with a stronger diurnal variation of soil temperature/liquid water in response to the atmospheric forcing.

285
Such discrepancies were considerably reduced by model with the coupled water-heat physics. Surface evapotranspiration was overestimated by unCPLD model, mainly due to the mismatches during the initial freezing and freezing-thawing transition period. CPLD models, with the coupled constraints from the perspective of water and energy conservation, significantly improve the model performance in mimicking the surface evapotranspiration dynamics during frozen period. The analysis of heat budget components and 290 latent heat flux density revealed that the improvement of soil temperature simulations by CPLD model is ascribed to its physical consideration of vapor flow and thermal effect on water flow, with the former mainly functions at regions above the evaporative front, the latter dominates below the evaporative front. The contribution of airflow induced water and heat flow to the total mass and energy fluxes is negligible. However, given the explicit consideration of airflow, the latent heat flux and its effect on heat transfer were enhanced 295 during the freezing-thawing transition period. This work highlighted the role of considering the vapor flow, thermal effect on water flow, and airflow in portraying the subsurface soil hydrothermal dynamics especially during freezing-thawing transition periods. To sum up, this study can contribute to a better understanding of freeze-thaw mechanisms of permafrost, which will subsequently contribute to the quantification of https://doi.org/10.5194/hess-2020-253 Preprint. Discussion started: 15 June 2020 c Author(s) 2020. CC BY 4.0 License.

300
FT model is to be coupled with a biogeochemical model, as lately implemented (Yu et al., 2020).