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

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.

1 Introduction

Frozen soils have been reported to undergo significant changes under climate warming (Cheng and Wu, 2007; Hinzman et al., 2013; Biskaborn et al., 2019; Zhao et al., 2019). Changes in the freezing and thawing process can alter soil hydrothermal regimes and water flow pathways and, thus, affect vegetation development (Walvoord and Kurylyk, 2016). Such changes will further considerably affect the spatial pattern, the seasonal to interannual variability, and long-term trends in land surface water, energy and carbon budgets, and then the land surface–atmosphere interactions (Subin et al., 2013; Iijima et al., 2014; Schuur et al., 2015; Walvoord and Kurylyk, 2016). Understanding the soil freezing and thawing processes appears to be the necessary path for better water resources management and ecosystem protection in cold regions.

When soil experiences the freezing and thawing process, there is a dynamic thermal equilibrium system of ice, liquid 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, and heat flow was first proposed by Philip and De Vries (1957; hereafter termed PdV57), considering the enhanced vapor transport. The PdV57 theory has been widely applied for a detailed understanding of soil evaporation during the drying process (De Vries, 1958, 1987; Milly, 1982; Saito et al., 2006; Novak, 2010). The attempts to simulate the coupled water and heat transport in frozen soils started in the 1970s (e.g., Harlan, 1973; Guymon and Luthin, 1974). Since then, numerical tools for simulating 1D frozen soil were gradually developed. Flerchinger and Saxton (1989) developed the simultaneous heat and water (SHAW) model with the capacity to simulate the coupled water and heat transport process. Hansson et al. (2004) accounted for the phase changes in the HYDRUS-1D model and verified its numerical stability with rapidly changing boundary 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 called MarsFlo. Aiming to efficiently deal with the water phase change between liquid and ice, the enthalpy-based frozen soil model (using enthalpy and total water mass instead of temperature and liquid water content as the prognostic variables) was developed and demonstrated its capability to stably and efficiently simulate the soil freezing and thawing process (Li et al., 2010; Bao et al., 2016; Wang et al., 2017). These models work together with other modifications and simplifications to generate a hierarchy of frozen soil models; see the detailed review by Li et al. (2010) and Kurylyk and Watanabe (2013).

Airflow has been reported as being important to the soil water and heat transfer process under certain conditions (Touma and Vauclin, 1986; Prunty and Bell, 2007). Zeng et al. (2011a, b) found that soil evaporation is enhanced after precipitation events, by considering airflow, and demonstrated that the air-pressure-induced advective fluxes inject the moisture into the surface soil layers and increase the hydraulic conductivity at the top layer. The diurnal variations in air pressure resulted in the vapor circulation between the atmosphere and the land surface. Wicky and Hauck (2017) reported that the temperature difference between the upper and the lower 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 in the presence of soil ice. The abovementioned studies demonstrate that the explicit consideration of airflow has the potential to affect the soil hydrothermal regime. However, to what extent and under what conditions airflow plays significant roles in the subsurface heat budgets have not been detailed.

Current land surface models (hereafter LSMs), however, usually adopted simplified frozen soil physics with relatively coarse vertical discretization (Koren et al., 1999; Viterbo et al., 1999; Niu et al., 2011; Swenson et al., 2012). 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 absent. This mostly leads to oversimplifications of physical representations of hydrothermal and ecohydrological dynamics in cold regions (Novak, 2010; Su et al., 2013; Wang et al., 2017; Cuntz and Haverd, 2018; Grenier et al., 2018; Wang and Yang, 2018; Qi et al., 2019). Specifically, Su et al. (2013) evaluated the European Centre for Medium-Range Weather Forecasts (ECMWF) soil moisture analyses over the Tibetan Plateau and found that HTESSEL cannot capture phase transitions of soil moisture (i.e., underestimation during the frozen period while overestimation during thawing). There are continuous efforts to improve parameterizations and representations of cold region dynamics, including frozen ground (Boone et al., 2000; Luo et al., 2003), vapor diffusion (Karra et al., 2014), thermal diffusion (Bao et al., 2016), coupling water and heat transfer (Wang and Yang, 2018), and three-layer snow physics (Wang et al., 2017; Qi et al., 2019). While, to our knowledge, few studies have investigated the role of increasing complexities in soil physical processes (from the basic coupled to the advanced coupled water and heat transfer processes and, then, to the explicit consideration of airflow) in simulating the thermohydrological states in cold regions. How, and to what extent, do the complex mutual dependent physics affect the soil mass and energy transfer in frozen soils? Is it necessary to consider a fully coupled physical process in the LSMs? These two questions frame the scope of this work.

In this paper, we incorporated the various complexities of soil water and heat transport mechanisms into a common modeling framework (namely, the simultaneous transfer of energy, momentum, and mass in unsaturated soils with freeze–thaw – STEMMUS–FT). 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 and latent heat flux density analyses were further carried out to illustrate the underlying mechanisms of different coupled soil water–heat physics. Section 2 describes the experimental site and three different complexities of subsurface physics within the STEMMUS framework. The performance of different models is presented in Sect. 3, together with the subsurface heat budgets and latent heat flux density analyses. Section 4 discusses the effects of considering coupled soil water–heat transfer and airflow in frozen soils. The conclusion is drawn in Sect. 5.

Table 1 Governing equations for different complexities of water and heat coupling physics (see the Appendix for notations).

2 Methodology

2.1 Experimental site

Maqu station, equipped with a catchment-scale soil moisture and soil temperature (SMST) monitoring network and micrometeorological observing system (Su et al., 2011; Dente et al., 2012; Zeng et al., 2016), is situated on the northeastern 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 a dry winter and warm summer (Dwb). The mean annual air temperature is 1.2 ^∘C, and the mean air temperatures of the coldest month (January) and warmest month (July) are about −10.0 and 11.7 ^∘C, respectively. Alpine meadows (e.g., Cyperaceae and Gramineae), with heights varying from 5 to 15 cm throughout the growing season, are the dominant land cover in this region. In situ soil sampling determined that the soil has a mixture of sandy loam, silt loam, and organic soil, with a maximum of 18.3 % organic matter for the upper soil layers (Dente et al., 2012; Zheng et al., 2015a; Zhao et al., 2018).

The Maqu SMST monitoring network spans an area of approximately 40 km × 80 km, with the elevation ranging from 3200 to 4200 m above sea level (a.s.l.). SMST profiles are automatically measured by 5TM ECH₂O probes (Meter Group, Inc., USA) installed at the soil depths of 5, 10, 20, 40, and 80 cm. The micrometeorological observing system includes a 20 m planetary boundary layer (PBL) tower providing the meteorological measurements at five heights above ground (i.e., wind speed and direction, air temperature, and relative humidity), and an eddy covariance (EC150; Campbell Scientific, Inc., USA) system installed for measuring the turbulent, sensible, latent heat fluxes, and carbon fluxes. The four-component down and upwelling solar and thermal radiation (NR01-L; Campbell Scientific, Inc., USA) and liquid precipitation (T200B; Geonor, Inc., USA) are also monitored.

2.2 Mass and energy transport in unsaturated soils

On the basis of the STEMMUS modeling framework, the increasing complexity of vadose zone physics in frozen soils was implemented as three alternative models (Table 1). First, STEMMUS enabled isothermal water and heat transfer physics (Eqs. 1 and 2). The 1D Richards equation was utilized to solve the isothermal water transport in variably saturated soils. The heat conservation equation took into account the freezing and thawing process and the latent heat due to the water phase change. The effect of soil ice on soil hydraulic and thermal properties was considered. This is termed the basic coupled water and heat transfer model (BCM).

Second, the fully coupled water and heat physics, i.e., water vapor flow and thermal effect on water flow, was explicitly considered in STEMMUS and termed the advanced coupled model (ACM). For the ACM physics, the extended version of the Richards (1931) equation, with modifications made by Milly (1982), was used as the water conservation equation (Eq. 3). Water flow can be expressed as liquid and vapor fluxes driven by both temperature gradients and matric potential gradients. The heat transport in frozen soils mainly includes heat conduction (CHF; ${\mathit{\lambda }}_{\mathrm{eff}}\frac{\partial T}{\partial z}$), convective heat transferred by liquid flux (HFL; $-C_{\mathrm{L}}{q}_{\mathrm{L}}(T-T_{\mathrm{r}})$, $-C_{\mathrm{L}}S(T-T_{\mathrm{r}})$), vapor flux (HFV; $-C_{\mathrm{V}}{q}_{\mathrm{V}}(T-T_{\mathrm{r}})$), the latent heat of vaporization (LHF; −q_VL₀), the latent heat of freezing or thawing (−ρ_iθ_iL_f), and a source term associated with the exothermic process of the wetting of a porous medium (integral heat of wetting) ($-{\mathit{\rho }}_{\mathrm{L}}W\frac{\partial {\mathit{\theta }}_{\mathrm{L}}}{\partial t}$). It can be expressed as Eq. (4) (De Vries, 1958; Hansson et al., 2004).

Lastly, STEMMUS expressed the freezing soil porous medium as the mutually dependent system of liquid water, water vapor, ice water, dry air, and soil grains in which, other than airflow, all other components are kept the same as in ACM (termed the ACM–AIR model; Eqs. 5–7; Zeng et al., 2011a, b; Zeng and Su, 2013). The effect of airflow on soil water and heat transfer can be two-fold. First, the airflow-induced water and vapor fluxes (q_La, q_Va) and the corresponding convective heat flow (HFa; $-q_{\mathrm{a}}{C}_{\mathrm{a}}(T-T_{\mathrm{r}})$) were considered. Second, the presence of airflow alters the vapor transfer processes, which can considerably affect the water and heat transfer in an indirect manner.

STEMMUS utilized the adaptive time step strategy, with maximum time steps ranging from 1 to 1800 s (e.g., with 1800 s as the time step under stable conditions). The maximum desirable change in soil moisture and soil temperature within one time step was set as 0.02 cm³ cm⁻³ and 2 ^∘C, respectively, to prevent too large a change in state variables that may cause numerical instabilities. If the changes between two adjacent soil moisture and temperature states are less than the maximum desirable change, STEMMUS continues without changing the length of current time step (e.g., 1800 s). Otherwise, STEMMUS will adjust the time step with a deduction factor, which is proportional to the difference between the changes that are too large and the maximum that are desirable in state variables. Within one single time step, the Picard iteration was used to solve the numerical problem, and the numerical convergence criteria is set as 0.001 for both soil matric potential (in centimeters) and soil temperature (in degrees Celsius).

To accommodate the specific conditions of a Tibetan meadow, the total depth of the soil column was set as 1.6 m (Fig. 1). The vertical soil discretization was designed to be finer for the upper soil layers (0.1–2.5 cm for 0–40 cm; 27 layers) than for the lower soil layers (5–20 cm for 40–160 cm; 10 layers). The surface boundary for the water transport was set as the flux-type boundary controlled by the atmospheric forcing (i.e., evaporation and precipitation), while the specific soil temperature was assigned as the surface boundary of the energy conservation equation. The free drainage (zero matric potential gradient) and measured soil temperature were set as the bottom boundary conditions for the water transport and heat transport, respectively. For the airflow, the surface boundary was set as the atmospheric pressure, and soil air was allowed to escape from the bottom of the soil column. Surface evapotranspiration was calculated using the Penman–Monteith method. Soil evaporation and transpiration can be separately estimated. The available radiation energy is partitioned into the canopy and soil component via the leaf area index (LAI); the canopy minimum surface resistance and soil surface resistance are then utilized to calculate the potential transpiration and soil evaporation. Actual transpiration is calculated as the function of potential transpiration and the root-length-density-weighted available soil liquid water (which is assumed to be zero if soil temperature falls below 0 ^∘C; Kroes et al., 2009; Orgogozo et al., 2019). For our simulation period, grassland stepped into the dormancy period as the soil froze. The accumulative positive temperature during the thawing period was not enough to break the dormancy of the vegetation. The contribution of plant transpiration to the land surface heat flux is negligible during the dormancy period. The effect of soil moisture on the actual soil evaporation is taken into account via the soil surface resistance (Eq. A6). All three aforementioned models adopted the same adaptive time step strategy and numerical solution and the same soil discretization, soil parameters (shown in Table 2), and boundary conditions. It indicated that the truncation errors, due to the numerical solution, among the three models were comparable. The differences among the models were mainly restricted to the various representations of soil physical processes (e.g., the inclusion of vapor flow and airflow or not).

Figure 1 (a) Conceptual illustration of the model setup, the surface and bottom boundary conditions, driving forces, and vertical discretization. (b) Half-hourly measurements of meteorological forcing, including air temperature (T_atm, ^∘C), relative humidity (HR_atm, %), net radiation (R_n, W m⁻²), wind speed (U_wind, m s⁻¹), and atmospheric pressure (P_atm, kPa) during the simulation period. Note that the dimensions are not drawn to scale; models were run at a 1D scale.

Table 2 The adopted average values of soil texture and hydraulic properties at different depths (see the Appendix for the notations).

3 Results

Given the same atmospheric forcing and the same set of parameters, the performance of models with varying complexities of soil water and heat physics was illustrated in Sect. 3.1–3.3. Section 3.4 and 3.5 further analyzed the variations in heat budgets and subsurface latent heat flux density, illustrating differences in the underlying mechanisms among various models.

Figure 2 Comparison of measured (Obs) and estimated time series of soil temperature at various soil layers using the basic coupled model (BCM), advanced coupled model (ACM), and advanced coupled model with airflow (ACM–AIR).

3.1 Soil hydrothermal profile simulations

The performance of the model, with various soil physics, simulating the soil thermal profile information is illustrated in Fig. 2. Both ACM and ACM–AIR reproduced the time series of the soil temperature at different soil depths well, except at 40 cm which is most probably due to inappropriate measurements (e.g., improper placement of sensors). However, there are significant discrepancies in soil temperature simulated by the BCM. Compared to the observations, a stronger diurnal behavior of soil temperature, in response to the fluctuating atmospheric forcing, was found, and an earlier stepping in and stepping out of the frozen period was simulated by the BCM. Such differences enlarged at deeper soil layers, with large bias and root mean square error (RMSE) values (Table 3).

Table 3 Comparative statistical values of observed and simulated soil temperature and moisture with three models, with the bold font indicating the best statistical performance.

Note: $\mathrm{Bias}=\frac{\sum _{i=\mathrm{1}}^n\left(y_i-\widehat{y_i}\right)}{n}$, $\mathrm{RMSE}=\sqrt{\frac{\sum _{i=\mathrm{1}}^n{\left(y_i-\widehat{y_i}\right)}^{\mathrm{2}}}{n}}$, where y_i, $\widehat{y_i}$ are the measured and model-simulated soil temperature or moisture, and n is the number of data points.

Figure 3 Comparison of measured (Obs) and model-simulated time series of soil moisture at various soil layers using the basic coupled model (BCM), advanced coupled model (ACM), and advanced coupled model with airflow (ACM–AIR).

Figure 3 presents the time series of observed and simulated soil liquid water content at five soil layers. During the rapid freezing period, a noticeable overestimation of the diurnal fluctuations and an early and fast decrease in soil liquid water content was simulated by BCM. Moreover, stronger diurnal fluctuations and an early increase in liquid water content were also found during the thawing period. The early thawing of soil water even led to an unrealistic refreezing process at 80 cm (from the 88th to 92nd day after December 2015), which is due to the simulated early warming of soil by BCM (Fig. 2). Such discrepancies were significantly ameliorated by ACM and ACM–AIR simulations. Nevertheless, all three models can capture the diurnal variations in and magnitude of liquid water content during the frozen period well. Note that there is an observable difference between ACM- and ACM–AIR-simulated soil liquid water content at shallower soil layers during the thawing process (e.g., Fig. 3; 5 cm).

3.2 Freezing front propagation

The time series of freezing front propagation derived from the measured and simulated soil temperature was reproduced in Fig. 4. Initialized from the soil surface, the freezing front quickly develops downwards until the maximum freezing depth is reached. The thawing process starts from both the top and bottom and is mainly driven by the atmospheric heat and bottom soil temperature, respectively. Such characteristics were captured well by both the ACM and ACM–AIR models in terms of freezing rate, maximum freezing depth, and surface thawing process, while the BCM tended to present a more fluctuating and rapidly freezing front propagation and a deeper maximum freezing depth that was reached early. The effect of atmospheric heat sources on the soil was overestimated by the BCM, as shown by the stronger diurnal early onset of the thawing process.

Figure 4 Comparison of measured (Obs) and model-simulated freezing front propagation (FFP) using the basic coupled model (BCM), advanced coupled model (ACM), and advanced coupled model with airflow (ACM–AIR). Note that the measured FFP was seen as the development of zero degree isothermal lines from the measured soil temperature field.

3.3 Surface evapotranspiration

The performance of the model with different soil physics in reproducing the latent heat flux dynamics is shown in Fig. 5. Compared to the observed LE, there is a significant overestimation of the half-hourly latent heat flux, which significantly degraded the overall performance when using BCM. The occurrence of such an overestimation was notably reduced when using ACM and ACM–AIR. The general underestimation of the latent heat flux by the ACM and ACM–AIR was found mostly during the freezing–thawing transition period (Fig. 6b) when the soil hydrothermal states are not well captured (Figs. 2 and 3).

Figure 5 Scatterplot of observed and model-estimated half-hourly latent heat flux using the (a) basic coupled model (BCM), (b) advanced coupled model (ACM), and (c) advanced coupled model with airflow (ACM–AIR). The color indicates the data composite of surface latent heat flux.

Figure 6 Comparison of observed and model-simulated (a) mean diurnal variations in surface evapotranspiration and (b) cumulative evapotranspiration (ET) by the basic coupled model (BCM), advanced coupled model (ACM), and advanced coupled model with airflow (ACM–AIR).

The overestimation of surface evapotranspiration by BCM was significant during the initial freezing and transition period (Fig. 6a; December and February). During the rapid freezing period (January), BCM presented a good match in the diurnal variation as compared to the observations. The monthly average diurnal variations were found to be well captured by ACM and ACM–AIR. Figure 6b shows the comparison of observed and simulated cumulative surface evapotranspiration. The overall overestimation of surface evapotranspiration by BCM can be clearly seen in Fig. 6b. Days in 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 ACM and ACM–AIR. A slight underestimation of cumulative surface evapotranspiration was simulated by ACM and ACM–AIR, with values of 3.98 % and 4.78 %, respectively.

3.4 Heat budgets

Figure 7 shows the time series of simulated energy budget components at 5 cm using BCM, ACM, and ACM–AIR during the freezing period (5th–11th day after 1 December) and the freezing–thawing transition period (83rd–89th day after 1 December). For the BCM, only the change rate of heat content (HC) and conductive heat flux divergence (CHF) are considered in the left-hand side (LHS) and right-hand side (RHS) of Eq. (2); see Table 1. Three additional terms, namely convective heat flux divergence of liquid flow (HFL) and vapor flow (HFV) and latent heat flux divergence, were included for the ACM, while for the ACM–AIR the convective heat flux divergence of airflow (HFa) was also added.

Figure 7 Time series of model-simulated heat budget components at the soil depth of 5 cm, using the (a, d) basic coupled model (BCM), (b, e) advanced coupled model (ACM), and (c, f) advanced coupled model with airflow (ACM–AIR) simulations during the typical 6 d freezing (a–c) and freezing–thawing transition (d–f) periods. HC – change rate of heat content; CHF – conductive heat flux divergence; HFL – convective heat flux divergence due to liquid water flow; HFV – convective heat flux divergence due to water vapor flow; HFa – convective heat flux divergence due to airflow; and LHF – latent heat flux divergence. Note that, for graphical purposes, HFL, HFV, HFa, and LHF were enhanced by a factor of 10 during the freezing period.

There is a strong diurnal variation in heat budget components (HC, CHF, and LHF; Table 1) corresponding to the diurnal fluctuation in soil temperature. For the BCM, the change rate of the heat content was almost completely balanced by the conductive heat flux divergence (CHF; Fig. 7a). Compared to the BCM, a stronger diurnal fluctuation of HC and CHF was found in ACM results. As inferred from the results in Fig. 2, the time series of soil temperature change ($\partial T/\partial t$) simulated by BCM was larger than that simulated by ACM. This indicates that BCM produced less fluctuation in apparent heat capacity ($C_{\mathrm{app}}={\mathit{\rho }}_{\mathrm{soil}}{C}_{\mathrm{soil}}+{\mathit{\rho }}_{\mathrm{i}}\frac{L_{\mathrm{f}}^{\mathrm{2}}}{gT}\frac{\mathrm{d}{\mathit{\theta }}_{\mathrm{L}}}{\mathrm{d}\mathit{\psi }}$) than ACM. During the freezing period, the latent heat flux divergence (LHF) was lower than the conductive heat flux divergence (CHF) by 1–2 orders of magnitude (Fig. 7b). The positive value of the LHF term during the daytime indicates that condensation happens at 5 cm as the water vapor moves downward. The convective heat fluxes of liquid flow and vapor flow were even smaller compared to conductive heat flux (Fig. 7b). There is no significant difference in heat budget components between ACM and ACM–AIR in terms of diurnal variation and magnitude. The convective heat flux divergence of airflow played a negligible role in the change in the thermal state (HC) (Fig. 7c).

The dynamics of heat balance components at the 5 cm soil layer were simulated for the freezing–thawing transition period (Fig. 7d–f). Both HC and CHF underwent strong diurnal variations with increasing fluctuation magnitude, indicating soil warming at 5 cm. For the ACM, CHF outnumbered HC during daytime, and the difference increased with time. Negative values were found for LHF and developed further over time. The sum of CHF and LHF nearly balanced the HC term. Such behavior was similarly reproduced by ACM–AIR, with a slightly larger difference between the HC and CHF terms. This means a larger amount of water vapor was evaporated from the 5 cm soil layer (with a more negative LHF term) from ACM–AIR simulations than that from ACM simulations, which explains the lower liquid water content for ACM–AIR (Fig. 3; 5 cm).

3.5 Subsurface latent heat flux density

To give more context to the results, the spatial and temporal distributions of the simulated latent heat flux density (S_h), namely $-{\mathit{\rho }}_{\mathrm{w}}L\partial q_{\mathrm{v}}/\partial z$, during the freezing and freezing–thawing transition period are shown in Fig. 8. For the BCM, the latent heat flux density (S_h) is not available as it neglects the vapor flow.

Figure 8 The spatial and temporal distributions of model-estimated soil latent heat flux density using the (a, d) advanced coupled model (ACM), (b, e) advanced coupled model with airflow (ACM–AIR), and (c, f) the difference between ACM and ACM–AIR simulations ($S_{\mathrm{h},\mathrm{ACM}\text{-}\mathrm{AIR}}-S_{\mathrm{h},\mathrm{ACM}}$) during the typical 6 d freezing and freezing–thawing transition periods. The left (a–c) and right (d–f) columns are for the freezing and freezing–thawing transition periods, respectively. Note that figures for the basic coupled model (BCM) are absent as the model cannot simulate the subsurface soil latent heat flux density.

Figure 8a shows that there is a strong diurnal variation of S_h in the upper 0.1 cm soil layers. Such diurnal behavior along the soil profile was interrupted at 1 cm, at which point the water vapor consistently moved upwards as an evaporation source (termed as evaporative front). The path of this upward water vapor was disrupted at 20 cm from 6 December onward when the freezing front developed. Compared to the upper 0.1 cm soil, a weaker diurnal fluctuation of S_h was found at lower soil layers. For ACM–AIR, the vapor transfer patterns were similar to those of ACM (Fig. 8b). There were isolated connections of condensed water vapor between the upper 1 cm of soil and the lower soil layers (S_h>0; e.g., 6, 7, 9, and 10 December), possibly associated with the downward airflow (see Fig. 12 in Yu et al., 2018). The large difference in magnitude of latent heat flux density between ACM and ACM–AIR appeared to be mainly isolated in the upper soil layers (Fig. 8c). At soil layers between 1 and 20 cm, ACM–AIR simulated less in the condensation vapor area (S_h>0) and more in the evaporation area (S_h<0), indicating that ACM–AIR produced an additional amount of condensation and evaporation water vapor compared with ACM (Fig. 8c).

Similar to during the freezing period, the S_h during the transition period can be characterized as having strong diurnal variations at upper soil layers, interrupting diurnal patterns with the constant upward evaporation of intermediate soil layers, and having weak diurnal variations at lower soil layers (Fig. 8d and e). While the maximum evaporation rate was less than that during the freezing period, the consistent evaporation zone developed to a depth of 5 cm. The path for the upward-moving water vapor tended to develop deeper than 30 cm with the absence of soil ice. The simulation by ACM–AIR produced more condensation and less evaporative water vapor than that of ACM (Fig. 8f). In addition, steadily more evaporative water vapor from 5 cm was simulated by ACM–AIR compared to ACM. This confirms the aforementioned point that, during the freezing–thawing transition period, large LHF values were simulated by ACM–AIR (Fig. 7).

4 Discussion

4.1 Coupled water and heat transfer processes

Vapor flow, which is dependent on soil matric potential and temperature, links soil water and heat transfer processes. The mutual dependence of soil water in different phases (liquid, water vapor, and ice) and heat transport on vapor flow enables the facilitation of our better understanding of the complex soil physical processes (e.g., Figs. 7–8). Specifically, the interdependence of soil moisture and soil temperature (SMST) profiles simulated by ACM was closer to the observation than that by BCM. In addition, a significant enhancement in the portrayal of the monthly average diurnal variations in surface evapotranspiration and cumulative evapotranspiration can be found in ACM simulations which constrain the hydrothermal regimes, especially during the freezing–thawing transition periods (Figs. 2, 3, and 6).

During the freezing period, liquid water in the soil freezes, which is analogous to the soil drying process, and water vapor fluxes instead of liquid fluxes dominate the mass transfer process (Zhang et al., 2016). Neglecting such an important water flux component unavoidably results in different or unrealistic simulations of surface evapotranspiration and SMST profiles (Li et al., 2010; Karra et al., 2014; Wang and Yang, 2018). Li et al. (2010) reported that vapor fluxes were comparable to the liquid water fluxes and affected the freezing and melting processes. On the basis of long-term 1D soil column simulations, Karra et al. (2014) reported that the inclusion of the vapor diffusion effect significantly increased the thickness of the ice layer, as explained by the positive vapor cold trapping–thermal conductivity feedback mechanism. From the energy budget perspective, latent heat fluxes contribute more, due to the vapor phase change (LHF), to the heat balance budget at soil layers above the evaporative front than that below it (see LHF in Fig. 7e vs. Fig. 7b; corresponding evaporative front shown in Fig. 8d vs. Fig. 8a). This is consistent with findings by Zhang et al. (2016), who concluded that the latent heat of vapor, due to phase change, is 2 orders of magnitude less than the heat fluxes, due to conduction during wintertime. This corresponds to our results in Fig. 7b and c during the freezing period, while our results further showed that the latent heat fluxes, due to vapor phase change, can be considerable during the transition period (Fig. 7e and f). The downward latent heat flux from ACM makes the subsurface soil warmer, which reduces the temperature gradient ($\partial T/\partial z$) (Wang and Yang, 2018). This further results in weaker diurnal fluctuations of the HC term in ACM than in BCM (see HC in Fig. 7e vs. Fig. 7d). At the soil layers below the evaporative front, the heat flux source from the vapor transfer process (LHF) is negligible (e.g., Fig. 7b). The thermal retard effect occurs as the presence of soil ice, expressed as the apparent heat capacity term (C_app), dominates the heat transfer process in frozen soils. By considering the thermal effect on water flow, ACM usually has a larger water capacity value $\partial \mathit{\theta }/\partial \mathit{\psi }$ than BCM does. As such, the intense thermal impedance effect leads to the result that ACM produced a weaker diurnal fluctuation of soil temperature than BCM at subsurface soil layers (e.g., Fig. 2; 20 cm).

4.2 Airflow in the soil

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

5 Conclusions

On the basis of the STEMMUS modeling framework, with various representations of water and heat transfer physics (BCM, ACM, and ACM–AIR), the performance of each model in simulating water and heat transfer and surface evapotranspiration was evaluated over a typical Tibetan meadow ecosystem. Results indicated that, compared to in situ observations, the BCM tended to present earlier freezing and thawing dates, with a stronger diurnal variation in soil temperature and liquid water in response to the atmospheric forcing. Such discrepancies were considerably reduced by the model with the advanced coupled water–heat physics. Surface evapotranspiration was overestimated by BCM, mainly due to the mismatches during the initial freezing and freezing–thawing transition period. ACM 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 the frozen period. The analysis of the heat budget components and latent heat flux density revealed that the improvement in soil temperature simulations by ACM is ascribed to its physical consideration of vapor flow and the thermal effect on water flow, with the former mainly functioning at regions above the evaporative front and the latter dominating below the evaporative front. The nonconductive heat process (liquid or vapor or air-induced heat convection flux) contributed very minimally to the total energy fluxes during the frozen period, except for the latent heat flux divergence at the topsoil layers. 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 during the freezing–thawing transition period. This work highlighted the role of considering the vapor flow and the 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 the freeze–thaw mechanisms of frozen soils, which will subsequently contribute to the quantification of permafrost carbon feedback (Burke et al., 2013; Kevin et al., 2014; Schuur et al., 2015) if the STEMMUS–FT model is coupled with a biogeochemical model, as has been implemented lately (Yu et al., 2020a).

Appendix A: Calculation of surface evapotranspiration

The one-step calculation of actual soil evaporation (E_s) and potential transpiration (T_p) is achieved by incorporating canopy minimum surface resistance and actual soil resistance into the Penman–Monteith model (i.e., the ET_dir method in Yu et al., 2016). LAI is implicitly used to partition available radiation energy into the radiation reaching the canopy and soil surface.

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{p}}={\displaystyle \frac{\mathrm{\Delta }\left(R_{\mathrm{n}}^{\mathrm{c}}-G\right)+{\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{p}}\frac{\left(e_{\mathrm{s}}-e_{\mathrm{a}}\right)}{r_{\mathrm{a}}^{\mathrm{c}}}}{\mathit{\lambda }\left(\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\frac{r_{\mathrm{c},\min }}{r_{\mathrm{a}}^{\mathrm{c}}}\right)\right)}}\text{(A2)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{s}}={\displaystyle \frac{\mathrm{\Delta }\left(R_{\mathrm{n}}^{\mathrm{s}}-G\right)+{\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{p}}\frac{\left(e_{\mathrm{s}}-e_{\mathrm{a}}\right)}{r_{\mathrm{a}}^{\mathrm{s}}}}{\mathit{\lambda }\left(\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\frac{r_{\mathrm{s}}}{r_{\mathrm{a}}^{\mathrm{s}}}\right)\right)}},\end{array}$$

where $R_{\mathrm{n}}^{\mathrm{c}}$ and $R_{\mathrm{n}}^{\mathrm{s}}$ (MJ m⁻² d⁻¹) are the net radiation at the canopy surface and soil surface, respectively; ρ_a (kg m⁻³) is the air density; c_p (J kg⁻¹ K⁻¹) is the specific heat capacity of air; $r_{\mathrm{a}}^{\mathrm{c}}$ and $r_{\mathrm{a}}^{\mathrm{s}}$ (s m⁻¹) are the aerodynamic resistance for canopy surface and soil surface, respectively; r_c,min (s m⁻¹) is the minimum canopy surface resistance; and r_s (s m⁻¹) is the soil surface resistance.

The net radiation reaching the soil surface can be calculated using Beer's law as follows:

$$\begin{array}{}\text{(A3)}& R_{\mathrm{n}}^{\mathrm{s}}=R_{\mathrm{n}}\exp (-\mathit{\tau }\mathrm{LAI}).\end{array}$$

And the net radiation intercepted by the canopy surface is the residual part of total net radiation, as follows:

$$\begin{array}{}\text{(A4)}& R_{\mathrm{n}}^{\mathrm{c}}=R_{\mathrm{n}}(\mathrm{1}-\exp (-\mathit{\tau }\mathrm{LAI}\left)\right).\end{array}$$

The minimum canopy surface resistance r_c,min is given by the following:

$$\begin{array}{}\text{(A5)}& r_{\mathrm{c},\min }=r_{\mathrm{l},\min }/{\mathrm{LAI}}_{\mathrm{eff}},\end{array}$$

where r_l,min is the minimum leaf stomatal resistance and LAI_eff is the effective leaf area index, which considers that, generally, the upper and sunlit leaves in the canopy actively contribute to the heat and vapor transfer.

The soil surface resistance can be estimated following van de Griend and Owe (1994):

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& r_{\mathrm{s}}=r_{\mathrm{sl}},{\mathit{\theta }}_{\mathrm{1}}>{\mathit{\theta }}_{\min },h_{\mathrm{1}}>-\mathrm{100}\mathrm{000}\mathrm{cm}\\ & r_{\mathrm{s}}=r_{\mathrm{sl}}{e}^{a\left({\mathit{\theta }}_{\min }-{\mathit{\theta }}_{\mathrm{1}}\right)},{\mathit{\theta }}_{\mathrm{1}}\le {\mathit{\theta }}_{\min },h_{\mathrm{1}}>-\mathrm{100}\mathrm{000}\mathrm{cm}\\ & r_{\mathrm{s}}=\mathrm{\infty },h_{\mathrm{1}}\le -\mathrm{100}\mathrm{000}\mathrm{cm},\end{array}\end{array}$$

where r_sl (10 s m⁻¹) is the resistance to molecular diffusion of the water surface, a (0.3565) is the fitted parameter, θ₁ is the topsoil water content, and θ_min is the minimum water content above which soil is able to deliver vapor at a potential rate.

The root water uptake term described by Feddes et al. (1978) is as follows:

$$\begin{array}{}\text{(A7)}& S\left(h\right)={\mathit{\rho }}_{\mathrm{L}}\mathit{\alpha }\left(h\right)S_{\mathrm{p}},\end{array}$$

where α(h) (dimensionless) is the reduction coefficient related to soil water potential h, and S_p (s⁻¹) is the potential water uptake rate.

$$\begin{array}{}\text{(A8)}& S_{\mathrm{p}}=b\left(z\right)T_{\mathrm{p}},\end{array}$$

where b(z) is the normalized water uptake distribution, which describes the vertical variation of the potential extraction term, S_p, over the root zone. Here the asymptotic function was used to characterize the root distribution as described in Gale and Grigal (1987), Jackson et al. (1996), Yang et al. (2009), and Zheng et al. (2015b). T_p is the potential transpiration in Eq. (A1).

Table A1 Notation.