STEMMUS-UEB v1.0.0: integrated modeling of snowpack and soil water and energy transfer with three complexity levels of soil physical processes

A snowpack has a profound effect on the hydrology and surface energy conditions of an area through its effects on surface albedo and roughness and its insulating properties. The modeling of a snowpack, soil water dynamics, and the coupling of the snowpack and underlying soil layer has been widely reported. However, the coupled liquid–vapor–air flow mechanisms considering the snowpack effect have not been investigated in detail. In this study, we incorporated the snowpack effect (Utah energy balance snowpack model, UEB) into a common modeling framework (Simultaneous Transfer of Energy, Mass, and Momentum in Unsaturated Soils with Freeze-Thaw, STEMMUS-FT), i.e., STEMMUS-UEB. It considers soil water and energy transfer physics with three complexity levels (basic coupled, advanced coupled water and heat transfer, and finally explicit consideration of airflow, termed BCD, ACD, and ACD-air, respectively). We then utilized in situ observations and numerical experiments to investigate the effect of snowpack on soil moisture and heat transfer with the abovementioned model complexities. Results indicated that the proposed model with snowpack can reproduce the abrupt increase of surface albedo after precipitation events while this was not the case for the model without snowpack. The BCD model tended to overestimate the land surface latent heat flux (LE). Such overestimations were largely reduced by ACD and ACD-air models. Compared with the simulations considering snowpack, there is less LE from no-snow simulations due to the neglect of snow sublimation. The enhancement of LE was found after winter precipitation events, which is sourced from the surface ice sublimation, snow sublimation, and increased surface soil moisture. The relative role of the mentioned three sources depends on the timing and magnitude of precipitation and the pre-precipitation soil hydrothermal regimes. The simple BCD model cannot provide a realistic partition of mass transfer flux. The ACD model, with its physical consideration of vapor flow, thermal effect on water flow, and snowpack, can identify the relative contributions of different components (e.g., thermal or isothermal liquid and vapor flow) to the total mass transfer fluxes. With the ACD-air model, the relative contribution of each component (mainly the isothermal liquid and vapor flows) to the mass transfer was significantly altered during the soil thawing period. It was found that the snowpack affects not only the soil surface moisture conditions (surface ice and soil water content in the liquid phase) and energy-related states (albedo, LE) but also the transfer patterns of subsurface soil liquid and vapor flow.

1 Introduction

In cold regions, the snowpack has a profound effect on hydrology and surface energy through its change of surface albedo, roughness, and insulating properties (Boone and Etchevers, 2001; Zhang, 2005). In contrast to rainfall, the melted snowfall enters the soil with a significant lag in time, and a large and sudden outflow or runoff may be produced because of the snowmelt effect. The heat-insulating property of snow cover also provides a buffer layer to reduce the magnitude of the underlying subsurface temperature variations and thus markedly affects the thickness of the active layer in cold regions. The effect of snow cover on the subsurface soils has been studied and reviewed (e.g., Zhang, 2005; Hrbáček et al., 2016). For instance, snow cover can act as an insulator between atmosphere and soil with its low thermal conductivity (Zhang, 2005; Hrbáček et al., 2016). The snowmelt functions as the energy sink via the absorption of heat due to phase change (Zhang, 2005). Yi et al. (2015) investigated the seasonal snow cover effect on the soil freezing and thawing process and its related carbon implications. Such studies mainly focus on the thermal effect of snowpack on the frozen soils. However, the effect of snowpack on the soil water and vapor transfer process is rarely reported (Hagedorn et al., 2007; Iwata et al., 2010; Domine et al., 2019).

Table 1 Brief overview of current soil–snow modeling efforts.

The abbreviations used in the table are as follows: HT_cond, heat conduction; Advc, advection; LH_phas, latent heat due to phase change; HT_Convect, convective heat due to liquid; SHP, soil physical process; Albedo_SNW_1A, snow albedo 1A, a function of snow age; Albedo_SNW_1B, snow albedo 1B, empirical function considering dry and wet states; Albedo_SNW_1C, snow albedo 1C, a function of extinction coefficient, grain size, and solar zenith angle; Albedo_SNW_2, snow albedo 2, a two-stream radiative transfer solution considering snow aging, solar zenith angle, optical parameters, and impurity; Albedo_SNW_3A, snow albedo 3A, prognostic snow albedo considering aging effect; Albedo_SNW_3B, snow albedo 3B, prognostic snow albedo considering aging effect and vegetation type dependent; Albedo_SNW_3C, snow albedo 3C, prognostic snow albedo considering aging and optical diameter; Albedo_SNW_3D, snow albedo 3D, prognostic snow albedo considering age and microstructure; Albedo_SNW_3E, snow albedo 3E, prognostic snow albedo considering aging effect and dry and wet states; Albedo_SNW_3F, snow albedo 3F, prognostic snow albedo considering aging effect and solar zenith angle; Albedo_SNW_4, snow albedo 4, diagnostic snow albedo considering snow aging, sleet and snowfall fraction, grain diameter, cloud fraction, and solar elevation effect; Density_SNW_1, snow density 1 relying on in situ measurements; Density_SNW_2A, snow density 2A, a function of air temperature; Density_SNW_2B, snow density 2B, a function of extinction coefficient and grain-size; Density_SNW_2C, snow density 2C, a function of old (densification), newly fallen (air temperature) snow pack density, and snow depth; Density_SNW_3, snow density 3, diagnostic density considering wet-bulb temperature; Density_SNW_4A, snow density 4A, prognostic density considering temperature, wind effect, snow compaction, and water and ice states; Density_SNW_4B, snow density 4B, prognostic density considering overburden and thermal metamorphisms; Density_SNW_4C, snow density 4C, prognostic snow density considering snow compaction and settling; Density_SNW_4D, snow density 4D, prognostic snow density considering snow compaction and wind-induced densification; Density_SNW_4E, snow density 4E, prognostic snow density considering snow compaction, settling, and vapor transfer; Density_SNW_4F, snow density 4F, prognostic density, a function of wind speed and air temperature; Density_SNW_4G, Snow density 4G, prognostic density, a function of stress state and microstructure; Density_SNW_4H, Snow density 4H, prognostic density considering snow temperature.

A great amount of effort has been made to better reproduce the snowpack characteristic and its effects in models. Initially, snowpack dynamics were expressed as a simple function of temperature. Nevertheless, these empirical relations have limited applications in complex climate conditions (Pimentel et al., 2015). Many physically based models for the mass and energy balance in the snowpack have been developed for their coupling with hydrological models or atmospheric models. Boone and Etchevers (2001) divided these snow models into three main categories: (i) simple force-restore schemes with the snow modeled as the composite snow–soil layer (Pitman et al., 1991; Douville et al., 1995; Yang et al., 1997) or a single explicit snow layer (Verseghy, 1991; Tarboton and Luce, 1996; Slater et al., 1998; Sud and Mocko, 1999; Dutra et al., 2010); (ii) detailed internal snow process schemes with multiple snow layers of fine vertical resolution (Jordan, 1991; Lehning et al., 1999; Vionnet et al., 2012; Leroux and Pomeroy, 2017); and (iii) intermediate-complexity schemes with physics from the detailed schemes but with a limited number of layers, which are intended for coupling with atmospheric models (e.g., Sun et al., 1999; Boone and Etchevers, 2001). The intercomparison results of the abovementioned snow models at an alpine site indicated that all three types of schemes are capable of representing the basic features of the snow cover over the 2-year period but behaved differently on shorter timescales. Furthermore, the Snow Model Intercomparison Project (SnowMIP) at two mountainous alpine sites revealed that the albedo parameterization was the major factor influencing the simulation of net shortwave radiation. Though this parameterization is independent of model complexity (Etchevers et al., 2004) it directly affects the snow simulation. SnowMIP2 evaluated 33 snowpack models across a wide range of hydrometeorological and forest canopy conditions. It identified the shortcomings of different snow models and highlighted the necessity of studying the separate contribution of individual components to the mass and energy balance of snowpack (Rutter et al., 2009). With the majority of research focused on the intercomparison of the snowpack models with various physical complexities, little attention has been paid to the treatment of the underlying soil physical processes (see the brief overview of current soil–snow modeling efforts in Table 1).

In current soil–snow modeling research, soil water and heat transfer are usually not fully coupled, and moreover the vapor flow and airflow are absent (Koren et al., 1999; Niu et al., 2011; Swenson et al., 2012). This may lead to the unrealistic interpretation of the underlying soil physical processes and the snowpack energy budgets (Su et al., 2013; Wang et al., 2017). Researchers have emphasized the need to consider the coupled soil water and heat transfer mechanisms (Scanlon and Milly, 1994; Bittelli et al., 2008; Zeng et al., 2009a, b; Yu et al., 2018a). As a consequence, dedicated efforts have been made to implement it in the recent updated models (e.g., Painter et al., 2016; Wang et al., 2017; Cuntz and Haverd, 2018). On the other hand, the role of the airflow has been reported as being important in many relevant studies, including retarding soil water infiltration (Touma and Vauclin, 1986; Prunty and Bell, 2007), enhancing surface evaporation after precipitation (Zeng et al., 2011a, b), enlarging the temperature difference between the upper and lower part of a permafrost talus slope (Wicky and Hauck, 2017), interacting with soil ice and vapor components, and enhancing the vapor transfer in frozen soils (Yu et al., 2018a, 2020c). However, to our knowledge, few soil–snow models have taken into account the soil–dry air transfer processes and moreover the multi-parameterization of the soil physical processes (from the basic coupled to the advanced coupled water and heat transfer processes and then to the explicit consideration of airflow), resulting in the lack of understanding on how and to what extent the complex soil physics affect the model interpretation of the snowpack effects.

In this paper, one of the widely used snowpack models (Utah energy balance snowpack model, UEB, Tarboton and Luce, 1996) was incorporated into a common soil modeling framework (Simultaneous Transfer of Energy, Mass and Momentum in Unsaturated Soils with Freeze-Thaw, STEMMUS-FT, Zeng et al., 2011a, b; Zeng and Su, 2013; Yu et al., 2018a). The new model is named STEMMUS-UEB and is configured with various levels of model complexity in terms of mass and energy transport physics. We utilized in situ observations and numerical experiments with STEMMUS-UEB to investigate the effect of snowpack on the underlying soil mass and energy transfer with different complexities of soil models. The description of the coupled soil–snow modeling framework STEMMUS-UEB and the model setup for this study are presented in Sect. 2. Section 3 verifies the proposed model and identifies the effect of snowpack on soil liquid–vapor fluxes. The uncertainties and limitations of this study and the applicability of the proposed model are discussed in Sect. 4.

2 Description of the coupled soil–snow modeling framework and model setup

This section first presents the coupling procedure of STEMMUS-FT and UEB model, followed by the detailed description of the two models and their successful applications. Then the used model configurations and two tested experimental sites in the Tibetan Plateau were elaborated. The Maqu case is for investigating the effect of snowpack on the underlying soil hydrothermal regimes. The Yakou case is for demonstrating the validity of the developed STEMMUS-UEB model in reproducing the snowpack dynamics (results were presented in Appendix B). In addition, the relationship between the snow cover properties and albedo was presented in Appendix B4, which confirmed the validity of using the albedo to identify the presence of snowpack and its lasting time.

2.1 Coupling procedure

The coupled process between the snowpack model (UEB) and the soil water model (STEMMUS-FT) was illustrated in Fig. 1. The sequential coupling is employed to couple the soil model with the current snowpack model. The role of the snowpack is explicitly considered by altering the water and heat flow of the underlying soil. The snowpack model takes the atmospheric forcing as the input (precipitation, air temperature, wind speed and direction, relative humidity, shortwave and longwave radiation) and solves the snowpack energy and mass balance (Eqs. A8 and A9; subroutines: ALBEDO, PARTSNOW, PREDICORR), which provides the melt water flux and heat flux as the surface boundary conditions for the soil model STEMMUS-FT (subroutines: h_sub and Enrgy_sub for the advanced coupled models and Diff_Moisture_Heat for the basic coupled model). The soil–snow coupling variables are the snowmelt water flux M_r, the convective heat flux due to snowmelt water Q_m, and the heat conduction flux Q_g. STEMMUS-FT then solves the energy and mass balance equations of soil layers in one time step. To highlight the effect of the snowpack on the soil water and vapor transfer process, we constrained the soil surface energy boundary as the Dirichlet-type condition (take the specific soil temperature as the surface boundary condition). Surface soil temperature was derived from the soil profile measurements and was not permitted to be higher than zero when there is snowpack. In such way, the reliability of the soil surface energy boundary condition is maintained and the snow thermal effect is implicitly considered. The snowmelt water flux, in addition to the rainfall, was added to the topsoil boundary for solving soil water transfer. To ensure numerical convergence, the adapted time step strategy was used. Half-hourly meteorological forcing measurements were linearly interpolated to the running time steps (Subroutine Forcing_PARM). The precipitation rate (validated at 3 h time intervals) was regarded uniformly within the 3 h duration (refer to Table S6.1 in the Supplement for details). The general description of the primary subroutines in STEMMUS-UEB was presented in Table 2. It includes the main functions, input and output, and their connections with other subroutines (linked with Tables S6.1 and S6.2 in the Supplement for the description of model input parameters and outputs for this study; see the detailed description in Tarboton and Luce, 1996; Zeng and Su, 2013; Yu et al., 2018a).

Figure 1 The overview of the coupled STEMMUS-FT and UEB model framework and model structure. SFCC is soil freezing characteristic curve, θ_L and θ_i are soil liquid water and ice content, K_Lh is soil hydraulic conductivity, and λ_eff is thermal conductivity. $\mathit{\psi },T,P_g$ are the state variables for soil module STEMMUS-FT (matric potential, temperature, and air pressure, respectively). U, SWE, and τ are the state variables for snow module UEB (snow energy content, snow water equivalent, and snow age, respectively). UEB is the Utah energy balance module. Precip, Ta, HRa, Rn, and u are the meteorological inputs (precipitation, air temperature, relative humidity, radiation, and wind speed, respectively). M_r is the snowmelt water flux, Q_m is the convective heat flux due to snowmelt water, and Q_g is the heat conduction flux. Model subroutines are in red.

Table 2 Main subroutines in STEMMUS-UEB.

Note that —> means the relevant subroutines that are incoming to the current one, and –> means the relevant subroutines for which the current subroutine is output to. agesn² and ALBEDO² mean the use of subroutines agesn and ALBEDO after solving the snowpack energy and mass conservation equations to update the snow age and albedo, respectively.

2.2 Soil mass and heat transfer module

The detailed physically based two-phase flow soil model (STEMMUS) was first developed to investigate the underlying physics of soil water, vapor, and dry air transfer mechanisms and their interaction with the atmosphere (Zeng et al., 2011a, b; Zeng and Su, 2013). It is achieved by simultaneously solving the balance equations of soil mass, energy, and dry air in a fully coupled way. The mediation effect of vegetation on such interactions was recently incorporated via the root water uptake sub-module (Yu et al., 2016) and by coupling with the detailed soil and vegetation biogeochemical process (Wang et al., 2021; Yu et al., 2020a). It facilitates our understanding of the hydrothermal dynamics of respective components in the frozen soil medium (i.e., soil liquid water, water vapor, dry air, and ice) by implementing the freeze–thaw process (hereafter STEMMUS-FT, for applications in cold regions, Yu et al., 2018a, 2020c).

The frozen soil physics considered in STEMMUS-FT includes three parts: (i) the ice blocking effect on soil hydraulic conductivities (see Sect. S2.2.2 in the Supplement), (ii) the inclusion of ice effect in the calculation of soil thermal capacity and conductivity (see Sect. S2.2.8), and (iii) the exchange of latent heat flux during phase change periods. With the aid of Clausius–Clapeyron relation, which characterizes the phase transition between liquid and solid phase in the thermal equilibrium system, the soil water characteristic curve (e.g., van Genuchten, 1980) is then extended to consider the freezing temperature dependence, i.e., soil freezing characteristic curve (Hansson et al., 2004; Dall'Amico et al., 2011). The fraction of soil liquid–solid water at a given temperature was then calculated prognostically with the soil freezing characteristic curve. Soil hydraulic parameters were further used in the Mualem (1976) model to compute the soil hydraulic conductivity. The ice effect is considered by reducing the soil saturated hydraulic conductivity as a function of ice content (Yu et al., 2018a).

In response to minimize the potential model-comparison uncertainties from various model structures and to figure out which process matters, three levels of complexity of mass and heat transfer physics are made available in the current STEMMUS-FT modeling framework (Yu et al., 2020c). First, the 1D Richards equation and heat conduction were deployed in STEMMUS-FT to describe the isothermal water flow and heat flow (termed BCD). The BCD model considers the interaction of soil water and heat transfer implicitly via the parameterization of heat capacity, thermal conductivity, and the water phase change effect. The water flow is fully affected by soil temperature regimes in the advanced coupled water and heat transfer model (termed ACD model). The movement of water vapor, as the primary linkage between soil water and heat flow, is explicitly characterized. STEMMUS-FT further enables the simulation of temporal dynamics of three water phases (liquid, vapor, and ice), together with the soil dry air component (termed ACD-air model). The governing equations of liquid water flow, vapor flow, airflow, and heat flow were listed in Appendix A1 (see the more detailed model description in Zeng et al., 2011a, b; Zeng and Su, 2013; Yu et al., 2018a, 2020c).

2.3 Snowpack module UEB

The Utah energy balance snowpack model (UEB; Tarboton and Luce, 1996) is a single-layer physically based snow accumulation and melt model. Two precipitation types, i.e., rainfall and snowfall, are discriminated by their dependence on air temperature. The snowpack is characterized using two primary state variables, snow water equivalent (SWE), and the internal energy U. Snowpack temperature is expressed diagnostically as the function of SWE and U together with the states of the snowpack (i.e., solid, solid and liquid mixture, and liquid). Given the insulation effect of the snowpack, snow surface temperature differs from the snowpack bulk temperature, which is mathematically considered using the equilibrium method (i.e., balances energy fluxes at the snow surface). The age of the snow surface, as the auxiliary state variable, is utilized to calculate the snow albedo (see Appendix A3). When the snowpack is shallow, the albedo is the weighting function of the snow albedo and the bare-ground albedo. The solar radiation penetration in the shallow snowpack is exponentially attenuated and expressed in the weighting factor. The melt outflow is calculated using Darcy's law with the liquid fraction as inputs. The conservation of mass and energy forms the physical basis of UEB (Tarboton and Luce, 1996, as presented in Appendix A2).

UEB is recognized as one simple yet physically based snowmelt model. It captures the snow process well (e.g., diurnal variation of meltwater outflow rate, snow accumulation, and ablation; see the general overview of UEB model development and applications in Table S6.3). It requires little effort in parameter calibration and can be easily transferable and applicable to various locations (e.g., Gardiner et al., 1998; Schulz and de Jong, 2004; Watson et al., 2006; Sultana et al., 2014; Pimentel et al., 2015; Gichamo and Tarboton, 2019), especially for data-scarce regions like the Tibetan Plateau. We thus selected the original parsimonious UEB (Tarboton and Luce, 1996) as the snow module to be coupled with the soil module (STEMMUS-FT).

2.4 Configurations of numerical experiments

On the basis of the aforementioned STEMMUS-UEB coupling framework, the various complexities of vadose zone physics were further implemented as three alternative model versions. First, the soil ice effect on soil hydraulic and thermal properties and the heat flow due to the water phase change were taken into account, while the water and heat transfer is not coupled in STEMMUS-FT and is termed the BCD model. Second, the STEMMUS-FT with the fully coupled water and heat transfer physics (i.e., water vapor flow and thermal effect on water flow) was applied and termed the ACD model. Lastly, on top of the ACD model, the air pressure was independently considered as a state variable (therefore, the airflow) and termed the ACD-air model. With the abovementioned model versions (STEMMUS-FT_Snow), taking into account the no-snow scenarios (STEMMUS-FT_No-Snow), Table 3 lists the configurations of all six designed numerical experiments. The model parameters used for all simulations for the tested experimental site are listed in Table S6.2.

Table 3 Numerical experiments with various mass and energy transfer schemes with and without explicit consideration of snow cover (Eqs. A1–A7 are listed in Appendix  A1; Eqs. (A8)–(A9) are listed in Appendix A2).

2.5 Description of the tested experimental sites

Maqu station, equipped with a catchment-scale soil moisture and soil temperature (SMST) monitoring network and micro-meteorological observing system, is situated on the northeastern edge of the Tibetan Plateau (Su et al., 2011; Dente et al., 2012; Zeng et al., 2016). According to the updated Köppen–Geiger climate classification system, it can be characterized as a cold climate with dry winter and warm summer. The annual mean precipitation is about 620 mm, and the annual average potential evaporation is about 1353.4 mm. Precipitation in Maqu is uneven over the year, with most of the precipitation events occurring from May to October and little precipitation or snowfall during the wintertime. The average annual air temperature is 1.2 ^∘C, and the mean air temperatures of the coldest month (January) and the warmest month (July) are about −10.0 and 11.7 ^∘C, respectively. Alpine meadows (e.g., Cyperaceae and Gramineae), with a height varying from 5 to 15 cm throughout the growing season, are the dominant land cover in this region. This site is seasonally snow covered, with temporal snow in the non-growing season, which is due to the intermittent snowfall and the rapid snow melting and sublimation caused by the high air temperature and strong solar radiation in the daytime. The general soil types are sandy loam, silt loam, and organic soil for the upper soil layers (Dente et al., 2012; Zheng et al., 2015; Zhao et al., 2018a). The soil texture and hydraulic properties were listed in Table S6.2, and how these were used in STEMMUS-UEB is illustrated in Fig. 1 and Table 2.

The Maqu SMST monitoring network spans an area of approximately 40 km × 80 km with an elevation ranging from 3200 m to 4200 m a.s.l. ($\mathrm{33}{}^{\circ }{\mathrm{30}}^{\prime }$–$\mathrm{34}{}^{\circ }{\mathrm{15}}^{\prime }$ N, $\mathrm{101}{}^{\circ }{\mathrm{38}}^{\prime }$–$\mathrm{102}{}^{\circ }{\mathrm{45}}^{\prime }$ E). SMST profiles are automatically measured by 5TM ECH₂O probes (METER Group, Inc., USA) installed at different soil depths, i.e., 5, 10, 20, 40, and 80 cm. The micro-meteorological observing system consists of 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 system (EC150, Campbell Scientific, Inc., USA) equipped for measuring the turbulent sensible and latent heat fluxes and carbon fluxes. The equipment for four-component downwelling and upwelling solar and thermal radiation (NR01-L, Campbell Scientific, Inc., USA) and liquid precipitation (T200B, Geonor, Inc., USA) are also deployed. A dataset from 1 December 2015 to 15 March 2016 was utilized in this study. Independent precipitation data (3 h time interval) during the same testing period from an adjacent meteorological station were used as the mutual validation data.

Yakou super snow station ($\mathrm{38}{}^{\circ }{\mathrm{00}}^{\prime }\mathrm{36}$ N, $\mathrm{100}{}^{\circ }{\mathrm{14}}^{\prime }\mathrm{24}$ E, 4145 m) is located in the upstream Heihe basin in the northeastern Tibetan Plateau. It is a high-elevation snow-covered site with the wet summers and dry winters. The dominant land type is tundra with frozen ground below. There is a unique seasonal variation of snow depth with the maximum snow depth usually being in the springtime (32 cm during the period 2014–2017). Loam is the main soil type with the silt loam near the surface and sandy soil for the deeper soil layers.

The integrated hydrometeorological, snow cover, and frozen ground data were published and available from the Cold and Arid Regions Science Data Center at Lanzhou (Che et al., 2019; Li et al., 2019; Li, 2019). The meteorological data (air temperature, wind speed, precipitation, downward shortwave and longwave radiation, and relative humidity) were recorded by the automatic meteorological station (AMS). In situ measurements of snow cover properties (snow depth and snow water equivalent) were obtained using the state-of-the-art instruments (SR50A and GammaMONitor, Campbell Scientific, USA). Soil moisture profiled at 4, 10, 20, 40, 80, 120, and 160 cm soil depth was measured using ECH2O-5 probes (METER Group, Inc., USA). In addition to the seven soil depths, the surface soil temperature (0 cm) was also recorded using the Avalon AV-10T sensors (Avalon Scientific, Inc., USA). The eddy covariance system was equipped at the Yakou site for measuring land surface turbulent fluxes. The dataset from 1 September to 31 December 2016 was used to validate the model performance in mimicking the dynamics of snow water equivalent, soil hydrothermal regimes, and land surface evaporation. The calibrated soil hydraulic and snow cover properties were listed in the Supplement in Table S6.2.

3 Results: comparison of simulation results of surface variables with and without snowpack effect

3.1 Albedo

The time series of surface albedo, calculated as the ratio of upwelling shortwave radiation to the downwelling shortwave radiation and estimated using BCD, ACD, and ACD-air models, is shown in Fig. 2 together with precipitation. As the snowpack has a higher albedo than the underlying surface (e.g., soil, vegetation) compared to the observations, models without snow module presented a relatively flat variation of daily average surface albedo and lacked the response to the winter precipitation events (Fig. 2, Table 4). With the snow module, STEMMUS-UEB models can mostly capture the abrupt increase of surface albedo after winter precipitation events. The mismatches in terms of the magnitude or absence of increased albedo after precipitation events indicated that the model tended to underestimate the albedo dynamics. The shallow snowfall events might be not well captured by the model (see Sect. 4.1). Three model versions (BCD-Snow, ACD-Snow, and ACD-air-Snow) produced similar fluctuations regarding the presence of snow cover with slight differences in terms of the magnitude of albedo.

Figure 2 Time series of observed and model simulated daily average albedo using (a) BCD, (b) ACD, and (c) ACD-air soil models with and without consideration of the snow module (including precipitation).

Table 4 Comparative statistics values of various model versions for snow albedo, LE, soil temperature, and soil moisture. The best statistical performance is highlighted in bold font​​​​​​​, while the values with poor statistical model performance are highlighted in italic font.

Note that BIAS $=\frac{\sum _{i=\mathrm{1}}^n\left(y_{\mathrm{i}}-\widehat{y_{\mathrm{i}}}\right)}{n}$, $R^{\mathrm{2}}=\mathrm{1}-\frac{{\sum }_{i=\mathrm{1}}^n{\left(y_{\mathrm{i}}-\widehat{y_{\mathrm{i}}}\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{\left(y_{\mathrm{i}}-\overline{y}\right)}^{\mathrm{2}}}$, and RMSE $=\sqrt{\frac{\sum _{i=\mathrm{1}}^n{\left(y_{\mathrm{i}}-\widehat{y_{\mathrm{i}}}\right)}^{\mathrm{2}}}{n}}$, where y_i and $\widehat{y_{\mathrm{i}}}$ are the measured and model simulated values of the selected variable (snow albedo, latent heat flux or LE, and soil temperature and moisture), $\overline{y}$ is the mean value of the measurements of the selected variable (snow albedo, LE, and soil temperature and moisture), and n is the number of data points. The correlations are all significant at the 0.01 level excluding values marked with “–”, which indicates that the correlation is not significant.

3.2 Soil temperature and moisture dynamics

The observed spatial and temporal dynamics of soil temperature from five soil layers were used to verify the performance of different models (Fig. 3). The initial soil temperature state can be characterized as the warm bottom and cool surface soil layers (based on in situ observations). The freezing front (indicated by the zero-degree isothermal line, ZDIL) developed downwards rapidly until the 70th day after 1 December 2015, when it reached its maximum depth. Following this, the freezing front stabilized as an offset effect of latent heat release (termed the zero-curtain effect). Such influence can be sustained until all the available water to that layer is frozen, at which point the latent heat effect is negligible compared to the heat conduction. At shallower layers, the atmospheric forcing dominates the fluctuation of thermal states. The isothermal lines (e.g., −2 ^∘C) had a larger variation than that of ZDIL. At deeper soil layers, the temporal dynamics of isothermal lines were smoother than that of ZDIL, indicating that the effect of fluctuated atmospheric force on soil temperature was damped with the increase of soil depth. Compared to the observations, BCD-Snow model presented an earlier development of the freezing front and arrival of the maximum freezing depth (60th day after 1 December 2015). The deeper and more fluctuated freezing front indicates that a stronger control of atmospheric forcing on soil thermal states was produced by BCD-Snow model. The ACD models can capture the propagation characteristic of the freezing front well in terms of the variation magnitude and maximum freezing depth. There is no significant difference in soil thermal dynamics between the model with and without the snow module, except at the surface soil layers (Table 4).

Figure 3 The spatial and temporal dynamics of observed (a) and simulated soil temperature using BCD, ACD, and ACD-air soil models both with and without consideration of the snow module (snow: b, e, h; no snow: c, f, i) and the difference between them (d, g, j) (simulations with snow minus simulations without snow). The red line indicates the zero-degree isothermal line (ZDIL) from the measured soil temperature. The observed soil freezing stage and stabilization stage is marked in Fig. 3a.

Figure 4 The spatial and temporal dynamics of observed (a) and simulated soil volumetric water content using BCD, ACD, and ACD-air soil models both with and without consideration of the snow module (snow: b, e, h; no snow: c, f, i) and the difference between them (d, g, j) (simulations with snow minus simulations without snow). The red line indicates the ZDIL from the measured soil temperature. The observed wet zone, dry zone, and rewet zone of soil moisture is indicated in Fig. 4a.

Figure 4 shows the spatial and temporal dynamics of observed and simulated soil water content in the liquid phase (SWCL). The SWCL of active layers depends to a large extent on the soil freezing and thawing status. Soil is relatively wet at soil layers of 10–60 cm for the starting period. Its temporal development was disrupted by the presence of soil ice and tended to increase wetness during the thawing period. A relatively dry zone (θ_L<0.06 m³ m⁻³) above the freezing front was found, indicating the nearly completely frozen soil during the stabilization stage. The initial wet zone of soil moisture was narrowed down and the rewetting zone tended to enlarge from BCD-Snow simulation due to its early freezing and thawing of soil (Fig. 4b). The position of the dry zone occurred earlier due to the early reaching of the stabilization period in the BCD-Snow model (Fig. 3b). For the ACD models, the position and development of initial wet zone, rewetting zone, and dry zone are similar to those from the observations, indicating the soil moisture dynamics can be captured well by the ACD models. Compared to the STEMMUS-FT_Snow model, there was no observable difference in the SWCL dynamics at deeper soil layers from STEMMUS-FT_No-Snow simulations. The surface SWCL was found affected from STEMMUS-FT_Snow simulations (Table 4).

3.3 Surface latent heat flux

Figure 5 shows the comparison of time series of observed and model-simulated surface cumulative latent heat flux using three models with and without consideration of the snow module. Considerable overestimation of latent heat flux was produced by the BCD-Snow model: 121.79 % more than was observed. Such overestimations were largely reduced by ACD and ACD-air models. There is a slight underestimation of cumulative latent heat flux in the ACD-Snow and ACD-air-Snow models, with values of −8.33 % and −7.05 %, respectively. Compared with STEMMUS-FT_Snow simulations, there is less latent heat flux produced by the STEMMUS-FT_No-snow simulations. This is mainly due to the sublimation of snow cover, which cannot be simulated by the STEMMUS-FT_No-snow models. The difference in cumulative latent heat flux between STEMMUS-FT with and without snow module increases from BCD to ACD-air schemes, with the values of 2.02 %, 7.69 %, and 8.97 % for BCD, ACD, and ACD-air schemes, respectively.

Figure 5 Time series of observed and model simulated surface cumulative latent heat flux (LE) using (a) BCD, (b) ACD, and (c) ACD-air soil models with and without consideration of the snow module (including precipitation). The top row shows the comparisons, and the bottom row shows the model bias of the cumulative surface LE.

3.4 Liquid and vapor fluxes

To further elaborate the effect of snowpack on LE, we presented the diurnal variations of LE and its components at two typical episodes with precipitation events (freezing and thawing period, respectively). The relative contribution of liquid and vapor flow to the total mass transfer after precipitation events was separately presented in Figs. 8 and 9, i.e., the liquid water flux driven by temperature q_LT, matric potential q_Lh and air pressure q_La, water vapor flux driven by temperature q_VT, matric potential q_Vh, and air pressure q_Va.

3.4.1  LE

Diurnal dynamics of the observed and simulated latent heat flux during the rapid freezing period with the occurrence of precipitation events, from 10th to 14th days after 1 December 2015, are shown as Fig. 6a, b, and c. Compared to the observations, the diurnal variations of latent heat flux were captured by the proposed model with various levels of complexities. Performance of BCD, ACD, and ACD-air models in simulating LE differed mainly regarding the magnitude and response to precipitation events. For the BCD-Snow model, the overestimation of LE was found at the 10th and 11th day after 1 December due to relatively high surface soil moisture simulation (Fig. S6.1b). A certain amount of enhanced surface evaporation was produced shortly after precipitation, which is most probably due to the snow sublimation. Snow sublimation does not appear to intuitively match with observations. The mismatch in the LE enhancement after precipitation events can be attributed to the fact that the partition process of precipitation into various components (rainfall, snowfall, canopy interception) might not be captured by the model well. Such a response to the winter precipitation events was absent from the BCD-No-Snow simulations.

Figure 6 Observed and model simulated latent heat flux using (a) BCD, (b) ACD, and (c) ACD-air soil models with and without the snow module of a typical 5 d freezing period (from the 10th to 14th day after 1 December 2015). P is the precipitation and Ps is the snowfall. All precipitation is in the form of snowfall.

The overestimation of LE was reduced by ACD and ACD-air models (Fig. 6b and c). Compared to the ACD-Snow model simulations, the ACD-No-snow model produced a stronger diurnal variation of LE after the precipitation and is approaching the measured LE. Lower diurnal variation of LE for the ACD-Snow model can be ascribed to the lower surface SWCL (see Fig. S6.1d and g). For the ACD-Snow model, precipitation was partitioned into rainfall and snowfall, part of which was directly evaporated as sublimation. The sum of rainfall and the melting part of snowfall reached the soil surface as the incoming water flux, which is less than that for the ACD-No-snow model (taking all the precipitation as the incoming water flux). There is no significant difference in the dynamics of LE between simulations by ACD models and ACD-air models.

During the thawing period, the diurnal variations of LE were simulated well by the models (Fig. 7). There are some discrepancies regarding the peak values of LE. For the BCD-Snow model, overestimations were found in 100th, 101st, and 102nd day after 1 December 2015. The high LE values on 100th and 101st day are probably due to the high surface soil moisture by the thawing water (Fig. S6.2b), whereas on the 102nd day it is due to the snow sublimation (Fig. 7a). The peak values were reproduced but shifted by BCD-No-Snow simulation, which occurred on 100th day and at the end of 102nd day, indicating the shift of surface soil moisture states (Fig. S6.2b).

Figure 7 Observed and model simulated latent heat flux using (a) BCD, (b) ACD, and (c) ACD-air soil models with and without the snow module of a typical 5 d thawing period (from the 100th to 104th day after 1 December 2015). P is the precipitation and Ps is the snowfall.

For the ACD model, the difference in latent heat flux between snow and no-snow simulations was noticeable 2 d after precipitation. The larger values of LE from the ACD-No-snow model occurred earlier than those from the ACD-Snow model due to the earlier response of surface soil moisture to the precipitation event (Fig. S6.2). Compared to the observations, the enhancement of LE advanced from the ACD-Snow simulations (Fig. 7b). This enhanced evaporation can be attributed to the snow sublimation and increased surface soil moisture content. Similar lag behavior of precipitation-enhanced evaporation was produced by the ACD-air-Snow models (Fig. 7c). There are mismatches in the time and magnitude of LE enhancement between ACD-Snow model simulations and observations (Fig. 7b). This discrepancy lies in the uncertainties of snowpack simulations, which can be attributed to either the inaccurate precipitation measurements (Barrere et al., 2017; Günther et al., 2019) or to the fact that the precipitation partition process is not described well by the model (Harder and Pomeroy, 2014; Ding et al., 2017).

3.4.2 LE and decomposition of surface mass transfer

During the freezing period, the soil water vapor rather than the liquid water flux dominated the surface mass transfer process. Missing the description of the vapor diffusion process hindered the BCD models ability to realistically depict the decomposition of surface mass transfer dynamics (Fig. 8a and b).

Figure 8 Model-simulated latent heat flux and surface soil (0.1 cm) thermal and isothermal liquid water and vapor fluxes (LE, q_VT, q_Vh, q_LT, q_Lh, q_La, q_Va) with and without the snow module of a typical 5 d freezing period (from the 10th to 14th day after 1 December 2015). Panels (a), (c), and (e) are the surface soil thermal and isothermal liquid water and vapor fluxes simulated by the BCD-Snow, ACD-Snow, and ACD-air-Snow models, respectively. Panels (b), (d), and (f) are the surface soil thermal and isothermal liquid water and vapor fluxes simulated by the BCD-No-Snow, ACD-No-Snow, and ACD-air-No-Snow models, respectively. LE is the latent heat flux, q_VT and q_Vh are the water vapor fluxes driven by temperature and matric potential gradients, respectively, q_LT and q_Lh are the liquid water fluxes driven by temperature and matric potential gradients, respectively, and q_La and q_Va are the liquid and vapor water fluxes driven by air pressure gradients, respectively. Positive and negative values indicate upward and downward fluxes, respectively. Note that the surface LE fluxes without snow sublimation are presented here. P is the precipitation, and Ps is the snowfall. All precipitation is in the form of snowfall.

There is a visible diurnal variation of thermal vapor flux q_VT from the ACD model simulation (Fig. 8c and d). The isothermal vapor flux q_Vh contributed to most of the mass transfer during the freezing period. It should be noted that the sum of water–vapor fluxes at the 0.1 cm soil layer cannot balance the surface evaporation, especially after the precipitation events (Fig. 8c). We assumed this premise and attributed it to the surface ice sublimation process. Precipitation water was frozen on the soil surface, and only vapor fluxes are active in the topsoil layers. Sublimation of surface ice may contribute to the gaps between liquid–vapor fluxes and LE (Yu et al., 2018a). As more precipitation water was frozen on the soil surface from the ACD-No-Snow model (Fig. 8d), the difference between the sum of water–vapor fluxes at the top 0.1 cm soil layer and the surface evaporative water enlarged compared to ACD-Snow simulations. Thermal liquid water flux q_LT appears negligible to the total mass flux during the whole simulation period. There is no significant difference recognized in the mass transfer between the ACD-air and ACD during the freezing period.

During the thawing period, a certain amount of upward liquid water flux was produced by the BCD model, supplying the water to the topsoil and evaporate into the atmosphere (Fig. 9a and b). Compared to the isothermal liquid flux q_Lh, the thermal liquid flux q_LT was negligible to the total mass flux.

Figure 9 Model-simulated latent heat flux and surface soil (0.1 cm) thermal and isothermal liquid water and vapor fluxes (LE, q_VT, q_Vh, q_LT, q_Lh, q_La, q_Va) using BCD (a, b), ACD (c, d), and ACD-air (e, f) simulations with and without the snow module, respectively, during the typical 5 d thawing periods (from the 100th to 104th day after 1 December 2015). Panel (a), (c), and (e) are the surface soil thermal and isothermal liquid water and vapor fluxes simulated by BCD-Snow, ACD-Snow, and ACD-air-Snow model, respectively. Panels (b), (d), and (f) are the surface soil thermal and isothermal liquid water and vapor fluxes simulated by BCD-No-Snow, ACD-No-Snow, and ACD-air-No-Snow model, respectively. LE is the latent heat flux, q_VT and q_Vh are the water vapor fluxes driven by temperature and matric potential gradients, respectively, q_LT and q_Lh are the liquid water fluxes driven by temperature and matric potential gradients, respectively, and q_La and q_Va are the liquid and vapor water fluxes driven by air pressure gradients, respectively. Positive and negative values indicate upward and downward fluxes, respectively. Note that the surface LE fluxes without snow sublimation are presented here. P is the precipitation, and Ps is the snowfall.

For the ACD model, the diurnal variation of thermal vapor flux q_VT was enhanced after precipitation, producing a larger amount of upward and downward vapor flux during the nighttime and daytime, respectively (e.g., Fig. 9c). As the surface soil is relatively dry, the isothermal vapor flux q_Vh contributes nearly all of the mass flux during the selected thawing period. Driven by the matric potential gradient, a large amount of isothermal water vapor flux q_Vh, accompanied by downward liquid water flux q_Lh, can be found after the nighttime precipitation event (Fig. 9c, d, e, f). These precipitation-induced isothermal liquid–vapor fluxes were lagged and less intense from the ACD-Snow model than that from the ACD-No-Snow model simulation (e.g., Fig. 9c vs. Fig. 9d). The snowpack reduces the instant precipitation infiltration process and enables the snowmelt afterwards, which led to the lagged and weaker response of surface SWCL to the precipitation (Fig. S6.2). This breaks the balance between isothermal vapor flux and evaporative LE (around the 103rd day after 1 December 2015). Compared to the ACD-No-Snow model, the imbalance was enlarged for the ACD-Snow model during the thawing period (Fig. 9c and d).

Compared to the ACD-No-Snow simulations, the upward thermal vapor flux q_VT was enhanced after precipitation for the ACD-air-No-Snow model (Fig. 9f). This enhanced upward vapor flux reduced the soil liquid water content at 0.1 cm (Fig. S6.2f) and decreased the soil hydraulic conductivity and then the downward isothermal liquid–vapor flux (q_Lh, q_Vh). Other than that there is no significant difference between the ACD-air model and the ACD model during the thawing period.

4 Discussion

4.1 Uncertainties in simulations of surface albedo and limitations

After a winter precipitation event, land surface albedo increases considerably (Fig. 2), indicating the presence of the snowpack. However, such snowfall events were episodic with small magnitudes (similar to those in Li et al., 2017), which means that they are difficult to capture well. Such difficulties can be partially attributed to the inherent uncertainties in precipitation measurements (both the precipitation amount and types). Due to the spatial variability of precipitation, the accurate observation of winter precipitation has proven to be a challenge, especially during windy winters (Barrere et al., 2017; Pan et al., 2017). It is necessary to have more snowpack-relevant measurements (e.g., the high-resolution measurements of the spatiotemporal field of wind speed, precipitation, and snowpack variations) to understand the dynamics of snowpack and its effect on energy and water fluxes. Furthermore, the temporal resolution of precipitation measurements adopted in this study is relatively coarse (3 h). In the current precipitation partition parameterization, the amount of snowfall was determined as a function of precipitation and air temperature thresholds. Given the coarse temporal resolution of precipitation measurements, the model may produce a time shift of snowfall events or even the misidentification of snowfall. The simple relation between the air temperature and precipitation types may be not suitable for this region because air temperature is not the best indicator of precipitation types, as argued by Ding et al. (2014). Other factors, i.e., relative humidity, surface elevation, and wet-bulb temperature, are also very relevant and should be taken into account for the discrimination of precipitation types. The other uncertainty lies in the representation of the snow process. For example, the wind blow effect and canopy snow interception, which have been recognized as important to the accurate simulation of snowpack dynamics (Mahat and Tarboton, 2014), are not taken into account in detail. Last but not least, the interpretation of surface albedo dynamics needs to be adapted to the specific site, especially regarding the shallow snow situations (Ueno et al., 2007, 2012; Ding et al., 2017; Wang et al., 2017). The albedo of the underlying surface should also be properly accommodated to this Tibetan meadow system. Regardless of the aforementioned uncertainties, our proposed model was capable of capturing the surface albedo variations with precipitation (Fig. 2) and can be seen as acceptable for analyzing snow cover effects in such a harsh environment.

4.2 Snow-cover-induced evaporation enhancement

In contrast to precipitation water from rainfall, precipitation water from snowfall enters the soil considerably lagged in time due to the water storage by snow cover (You et al., 2019). With the snow module, precipitation was partitioned into rainfall and snowfall. Part of the snowfall evaporated into the atmosphere as sublimation, and the other part, together with the rainfall, infiltrated into the underlying soil. It resulted in the delay of incoming water to the soil with a lower amount compared to that without consideration of the snow module. This amount of incoming water increased the evaporation after precipitation (Figs. 6 and 7). The other source for the enhanced evaporation flux after precipitation is snow sublimation, which is absent from the model without the snow module. Sublimation occurs readily under certain weather conditions (e.g., with freezing temperatures, enough energy). It can be more active in regions with low relative humidity, low air pressure, and dry winds. Such an amount of sublimation has been reported as being important from the perspective of climate and hydrology (e.g., Strasser et al., 2008; Jambon-Puillet et al., 2018), especially in high-altitude regions with low air pressure. During the freezing period, the evaporation enhancement can be also sourced from the sublimation of surface ice. The amount of the ice sublimation appeared to decrease during the freezing period in the presence of a transient snowpack (e.g., Fig. 8c vs. Fig. 8d). This is consistent with the results of Hagedorn et al. (2007), who investigated the effect of snow cover on the mass balance of ground ice with an artificially continuous annual snow cover. According to their results, the snow cover enhanced the vapor transfer into the soil and thus reduced the long-term ice sublimation. The relative contribution of increased surface soil moisture, snow sublimation, and surface ice sublimation to the enhanced evaporation is dependent on the pre-precipitation soil moisture and temperature states, air temperature, and the time and magnitude of precipitation events. Under the conditions of the low pre-precipitation SWCL with a freezing soil temperature (e.g., Fig. 8e, 11th vs. 12th day after 1 December), the precipitation falls on the surface as snowfall and rainfall (most freezes as ice). The sublimation from surface ice can contribute to most of the total mass transfer (e.g., Fig. 8e, 11th day after 1 December). If the soil temperature rises above the freezing temperature, there will be no sublimation of surface ice, in terms of contributing to the enhanced evaporation (e.g., Fig. 9e, 102nd day after 1 December).

4.3 Snow cover impacts with different soil model complexities

The model with different complexities of soil mass and energy transfer physics behaves differently in response to the winter precipitation events. During the freezing period, there is no significant difference in soil moisture simulated using the BCD models with and without the snow module. The precipitation water freezes at the soil surface, which cannot be transferred downwards with the BCD model physics. The sublimation, from either the snow or the surface ice, contributes to the precipitation-enhanced evaporation for the BCD model. As with vapor flow, the surface ice increases the soil moisture at lower layers via the downward isothermal vapor flux (Fig. 8). The surface ice sublimation and increased moisture-induced soil evaporation enhancement can be identified from the ACD model simulation. The role of airflow was negligible for the mass transfer during the freezing period.

When it comes to the thawing period, the BCD model produced a certain amount of liquid water flow, contributing considerably to the mass transfer. The obvious fluctuation of SWCL was noticed due to the thawing water and precipitation event. The main source for the increased evaporation was interpreted as isothermal liquid water flow, while for the ACD model the situation becomes more complex. Thawing surface ice and snowmelt water may coexist at the soil surface, resulting in different soil moisture response to precipitation events. The ice sublimation, snow sublimation, and increased soil moisture contribute to the evaporation enhancement after precipitation. When considering airflow, dry air interacts with soil ice and liquid and vapor water in soil pores (Yu et al., 2018a) and alters the soil moisture state. It thus considerably changes the relative contribution of each component to the mass transfer (Fig. 9).

5 Conclusions

With the aim to investigate the hydrothermal effect of the snowpack on the underlying soil system, we developed the integrated process-based soil–snow–atmosphere model, STEMMUS-UEB v1.0.0, which is based on the easily transferable and physically based description of the snowpack process and the detailed interpretation of the soil physical process with various complexities. From STEMMUS-UEB simulations, snowpack affects not only the soil surface conditions (surface ice and SWCL) and energy-related states (albedo, latent heat flux) but also the transfer patterns of subsurface soil liquid and vapor flow. STEMMUS-FT model can mostly capture the abrupt increase of surface albedo after winter precipitation events with consideration of the snow module. There is a significant overestimation of cumulative surface latent heat flux by the BCD model. The ACD and ACD-air models produce a slight underestimation of cumulative LE compared to the observations. Without sublimation from snowpack, there is less latent heat flux produced by STEMMUS-FT_No-snow simulations compared with snow simulations. The presence of snowpack alters the partition process of precipitation and thus the surface SWCL. BCD models with and without snowpack produced similar surface SWCL during the freezing period while resulting in an abrupt increase of soil moisture in response to the precipitation during the thawing period. The ACD-Snow model simulated a less intensive and lagged soil moisture variation in response to precipitation compared to the ACD-No-Snow model during both the freezing and thawing period, respectively. The ACD-air model affected the intensity of increased surface soil moisture, especially during the thawing period.

Three mechanisms, surface ice sublimation, snow sublimation, and increased soil moisture, can contribute to enhanced latent heat flux after winter precipitation events. The relative role of each mechanism in the total mass transfer can be affected by the time and magnitude of precipitation and pre-precipitation soil moisture and temperature states (see Sect. 4.3). The simple BCD model cannot provide a realistic partitioning of mass transfer. The ACD model, which takes into consideration vapor diffusion and thermal effect on water flow and snowpack, can produce a reasonable analysis of the relative contributions of different water flux components. When considering airflow, the relative contribution of each component to the mass transfer was substantially altered during the thawing period. Further work will take into account the thermal interactive effects between snowpack and the underlying soil, which explicitly considers the convective and conductive heat fluxes and the solar radiation attenuation due to the snowpack. Such work will inevitably enhance our confidence in interpreting the underlying mechanisms and physically elaborating on the role of snowpack in cold regions.

Appendix A

A1 STEMMUS-FT model with three levels of complexity

A1.1 Uncoupled soil water and heat transfer physics

The Richard equation, which describes the water flow under gravity and capillary forces in isothermal conditions, is solved for variably saturated soils.

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}=-{\displaystyle \frac{\partial q}{\partial z}}-S={\mathit{\rho }}_{\mathrm{L}}{\displaystyle \frac{\partial }{\partial z}}\left[K\left({\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}+\mathrm{1}\right)\right]-S,\end{array}$$

where θ (m³ m⁻³) is the volumetric water content, q (kg m⁻² s⁻¹) is the water flux, z (m) is the vertical-direction coordinate (positive upwards), S (s⁻¹) is the sink term for root water uptake, ρ_L (kg m⁻³) is the soil liquid water density, K (m s⁻¹) is the soil hydraulic conductivity, ψ (m) is the soil water potential, and t (s) is the time.

The heat conservation equation, considering the latent heat due to water phase change, can be expressed as follows:

$$\begin{array}{}\text{(A2)}& C_{\mathrm{soil}}{\displaystyle \frac{\partial T}{\partial t}}-{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{f}}{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{i}}}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\lambda }}_{\mathrm{eff}}{\displaystyle \frac{\partial T}{\partial z}}\right),\end{array}$$

where C_soil (J kg^−1 ∘C⁻¹) is the specific heat capacity of bulk soil, T (^∘C) is the soil temperature, ρ_i (kg m⁻³) is the density of soil ice, L_f (J kg⁻¹) is the latent heat of fusion, θ_i (m³ m⁻³) is the soil ice volumetric water content, and λ_eff (W m^−1 ∘C⁻¹) is the effective thermal conductivity of the soil.

A1.2 Coupled water and heat transfer

For the coupled water and heat transfer physics, the liquid water flow is non-isothermal and affected by soil temperature regimes. The movement of water vapor, as the linkage between soil water and heat flow, is explicitly characterized. With modifications made by Milly (1982), the extended version of Richards (1931) equation with consideration of the liquid and vapor flow is written as follows:

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& {\displaystyle \frac{\partial }{\partial t}}\left({\mathit{\rho }}_{\mathrm{L}}{\mathit{\theta }}_{\mathrm{L}}+{\mathit{\rho }}_{\mathrm{V}}{\mathit{\theta }}_{\mathrm{V}}+{\mathit{\rho }}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{i}}\right)=-{\displaystyle \frac{\partial }{\partial z}}\left(q_{\mathrm{L}}+q_{\mathrm{V}}\right)-S\\ & =-{\displaystyle \frac{\partial }{\partial z}}\left(q_{\mathrm{Lh}}+q_{\mathrm{LT}}+q_{\mathrm{Vh}}+q_{\mathrm{VT}}\right)-S\\ & ={\mathit{\rho }}_{\mathrm{L}}{\displaystyle \frac{\partial }{\partial z}}\left[K_{\mathrm{Lh}}\left({\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}+\mathrm{1}\right)+K_{\mathrm{LT}}{\displaystyle \frac{\partial T}{\partial z}}\right]\\ & +{\displaystyle \frac{\partial }{\partial z}}\left[D_{\mathrm{Vh}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}+D_{\mathrm{VT}}{\displaystyle \frac{\partial T}{\partial z}}\right]-S,\end{array}\end{array}$$

where ρ_V and ρ_i (kg m⁻³) are the density of water vapor and ice, respectively; θ_L and θ_V (m³ m⁻³) are the volumetric water content (liquid and vapor, respectively); q_L and q_V (kg m⁻² s⁻¹) are the soil water fluxes of liquid water and water vapor (positive upwards), respectively; K_Lh (m s⁻¹) and K_LT (m² s^−1 ∘C⁻¹) are the isothermal and thermal hydraulic conductivities, respectively; D_Vh (kg m⁻² s⁻¹) is the isothermal vapor conductivity; and D_VT (kg m⁻¹ s^−1 ∘C⁻¹) is the thermal vapor diffusion coefficient.

On the basis of the work of De Vries (1958) and Hansson et al. (2004), the heat transport function in frozen soils, considering the fully coupled water and heat transport physics, can be expressed as follows:

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& {\displaystyle \frac{\partial }{\partial t}}\left[\left({\mathit{\rho }}_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{s}}{C}_{\mathrm{s}}+{\mathit{\rho }}_{\mathrm{L}}{\mathit{\theta }}_{\mathrm{L}}{C}_{\mathrm{L}}+{\mathit{\rho }}_{\mathrm{V}}{\mathit{\theta }}_{\mathrm{V}}{C}_{\mathrm{V}}+{\mathit{\rho }}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{i}}{C}_{\mathrm{i}}\right)\left(T-T_{\mathrm{r}}\right)\right.\\ & \left.+{\mathit{\rho }}_{\mathrm{V}}{\mathit{\theta }}_{\mathrm{V}}{L}_{\mathrm{0}}-{\mathit{\rho }}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{i}}{L}_{\mathrm{f}}\right]-{\mathit{\rho }}_{\mathrm{L}}W{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{L}}}{\partial t}}\\ & ={\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\lambda }}_{\mathrm{eff}}{\displaystyle \frac{\partial T}{\partial z}}\right)-{\displaystyle \frac{\partial }{\partial z}}\left[q_{\mathrm{L}}{C}_{\mathrm{L}}\left(T-T_{\mathrm{r}}\right)\right.\\ & \left.+q_{\mathrm{V}}(L_{\mathrm{0}}+C_{\mathrm{V}}(T-T_{\mathrm{r}}\left)\right)\right]-C_{\mathrm{L}}S\left(T-T_{\mathrm{r}}\right),\end{array}\end{array}$$

where C_s, C_L, C_V, and C_i (J kg^−1 ∘C⁻¹) are the specific heat capacities of solids, liquid water, water vapor, and ice, respectively; ρ_s (kg m⁻³) is the density of solids; θ_s is the volumetric fraction of solids in the soil; T_r (^∘C) is the arbitrary reference temperature; L₀ (J kg⁻¹) is the latent heat of vaporization of water at the reference temperature T_r; and W (J kg⁻¹) is the differential heat of wetting (the amount of heat released when a small amount of free water is added to the soil matrix).

A1.3 Coupled mass and heat physics with airflow

In STEMMUS-FT, the temporal dynamics of three phases of water (liquid, vapor and ice), together with the soil dry air component are explicitly presented and simultaneously solved by spatially discretizing the corresponding governing equations of liquid water flow, vapor flow, and airflow.

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& {\displaystyle \frac{\partial }{\partial t}}\left({\mathit{\rho }}_{\mathrm{L}}{\mathit{\theta }}_{\mathrm{L}}+{\mathit{\rho }}_{\mathrm{V}}{\mathit{\theta }}_{\mathrm{V}}+{\mathit{\rho }}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{ice}}\right)\\ & =-{\displaystyle \frac{\partial }{\partial z}}\left(q_{\mathrm{Lh}}+q_{\mathrm{LT}}+q_{\mathrm{La}}+q_{\mathrm{Vh}}+q_{\mathrm{VT}}+q_{\mathrm{Va}}\right)-S\\ & ={\mathit{\rho }}_{\mathrm{L}}{\displaystyle \frac{\partial }{\partial z}}\left[K\left({\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}+\mathrm{1}\right)+D_{\mathrm{TD}}{\displaystyle \frac{\partial T}{\partial z}}+{\displaystyle \frac{K}{{\mathit{\gamma }}_w}}{\displaystyle \frac{\partial P_g}{\partial z}}\right]\\ & +{\displaystyle \frac{\partial }{\partial z}}\left[D_{\mathrm{Vh}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial z}}+D_{\mathrm{VT}}{\displaystyle \frac{\partial T}{\partial z}}+D_{\mathrm{Va}}{\displaystyle \frac{\partial P_g}{\partial z}}\right]-S,\end{array}\end{array}$$

where q_Lh, q_LT, and q_La (kg m⁻² s⁻¹) are the liquid water fluxes driven by the gradient of matric potential $\frac{\partial \mathit{\psi }}{\partial z}$, temperature $\frac{\partial T}{\partial z}$, and air pressure $\frac{\partial P_g}{\partial z}$, respectively. q_Vh, q_VT, and q_Va (kg m⁻² s⁻¹) are the water vapor fluxes driven by the gradient of matric potential $\frac{\partial \mathit{\psi }}{\partial z}$, temperature $\frac{\partial T}{\partial z}$, and air pressure $\frac{\partial P_g}{\partial z}$, respectively. P_g (Pa) is the mixed pore air pressure. γ_W (kg m⁻² s⁻²) is the specific weight of water; D_TD (kg m⁻¹ s^−1 ∘C⁻¹) is the transport coefficient for adsorbed liquid flow due to temperature gradient; D_Vh (kg m⁻² s⁻¹) is the isothermal vapor conductivity; D_VT (kg m⁻¹ s^−1 ∘C⁻¹) is the thermal vapor diffusion coefficient; and D_Va is the advective vapor transfer coefficient (Zeng et al., 2011a, b).

STEMMUS-FT takes into account different heat transfer mechanisms, including heat conduction (${\mathit{\lambda }}_{\mathrm{eff}}\frac{\partial T}{\partial z}$ ), convective heat transferred by liquid flux ($-C_{\mathrm{L}}{q}_{\mathrm{L}}(T-T_{\mathrm{r}})$, $-C_{\mathrm{L}}S(T-T_{\mathrm{r}})$), vapor flux ($-[L_{\mathrm{0}}{q}_{\mathrm{V}}+C_{\mathrm{V}}{q}_{\mathrm{V}}(T-T_{\mathrm{r}}\left)\right]$), and airflow (q_aC_a(T−T_r)). The latent heat of vaporization (ρ_Vθ_VL₀), the latent heat of freezing and thawing (−ρ_iθ_iL_f), and a source term associated with the exothermic process of wetting of a porous medium (integral heat of wetting)($-{\mathit{\rho }}_{\mathrm{L}}W\frac{\partial {\mathit{\theta }}_{\mathrm{L}}}{\partial t}$) are all considered here.

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& {\displaystyle \frac{\partial }{\partial t}}\left[\left({\mathit{\rho }}_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{s}}{C}_{\mathrm{s}}+{\mathit{\rho }}_{\mathrm{L}}{\mathit{\theta }}_{\mathrm{L}}{C}_{\mathrm{L}}+{\mathit{\rho }}_{\mathrm{V}}{\mathit{\theta }}_{\mathrm{V}}{C}_{\mathrm{V}}+{\mathit{\rho }}_{\mathrm{da}}{\mathit{\theta }}_{\mathrm{a}}{C}_{\mathrm{a}}\right.\right.\\ & \left.\left.+{\mathit{\rho }}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{i}}{C}_{\mathrm{i}}\right)\left(T-T_{\mathrm{r}}\right)+{\mathit{\rho }}_{\mathrm{V}}{\mathit{\theta }}_{\mathrm{V}}{L}_{\mathrm{0}}-{\mathit{\rho }}_{\mathrm{i}}{\mathit{\theta }}_{\mathrm{i}}{L}_{\mathrm{f}}\right]\\ & -{\mathit{\rho }}_{\mathrm{L}}W{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{L}}}{\partial t}}\\ & ={\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\lambda }}_{\mathrm{eff}}{\displaystyle \frac{\partial T}{\partial z}}\right)-{\displaystyle \frac{\partial }{\partial z}}\left[q_{\mathrm{L}}{C}_{\mathrm{L}}\left(T-T_{\mathrm{r}}\right)\right.\\ & \left.+q_{\mathrm{V}}(L_{\mathrm{0}}+C_{\mathrm{V}}(T-T_{\mathrm{r}}\left)\right)+q_{\mathrm{a}}{C}_{\mathrm{a}}\left(T-T_{\mathrm{r}}\right)\right]\\ & -C_{\mathrm{L}}S\left(T-T_{\mathrm{r}}\right),\end{array}\end{array}$$

where ρ_da (kg m⁻³) is the density of dry air, C_a (J kg^−1 ∘C⁻¹) is the specific heat capacity of dry air, and q_a (kg m⁻² s⁻¹) is the air flux. The airflow balance equation for solving the coupled water and heat equations is written as in Zeng et al. (2011a, b) and Zeng and Su (2013):

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}& {\displaystyle \frac{\partial }{\partial t}}\left[\mathit{\epsilon }{\mathit{\rho }}_{\mathrm{da}}\left(S_{\mathrm{a}}+H_{\mathrm{c}}{S}_{\mathrm{L}}\right)\right]={\displaystyle \frac{\partial }{\partial z}}\left[D_e{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{da}}}{\partial z}}\right.\\ & \left.+{\mathit{\rho }}_{\mathrm{da}}{\displaystyle \frac{S_{\mathrm{a}}{K}_g}{{\mathit{\mu }}_{\mathrm{a}}}}{\displaystyle \frac{\partial P_g}{\partial z}}-H_{\mathrm{c}}{\mathit{\rho }}_{\mathrm{da}}{\displaystyle \frac{q_{\mathrm{L}}}{{\mathit{\rho }}_{\mathrm{L}}}}+\left({\mathit{\theta }}_{\mathrm{a}}{D}_{\mathrm{Vg}}\right){\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{da}}}{\partial z}}\right],\end{array}\end{array}$$

where ε is the porosity, S_a ($=\mathrm{1}-S_{\mathrm{L}}$) is the degree of air saturation in the soil, S_L ($={\mathit{\theta }}_{\mathrm{L}}/\mathit{\epsilon }$) is the degree of saturation in the soil, H_c is Henry's constant, D_e (m² s⁻¹) is the molecular diffusivity of water vapor in soil, K_g (m²) is the intrinsic air permeability, μ_a (kg m⁻² s⁻¹) is the air viscosity, θ_a(=θ_V) is the volumetric fraction of dry air in the soil, and D_Vg (m² s⁻¹) is the gas-phase longitudinal dispersion coefficient.

A2 Snowpack module UEB

A2.1 Mass balance equation

The increase or decrease of snow water equivalence with time equals the difference of income and outgoing water flux:

$$\begin{array}{}\text{(A8)}& {\displaystyle \frac{\mathrm{dSWE}}{\mathrm{d}t}}=P_{\mathrm{r}}+P_{\mathrm{s}}-M_{\mathrm{r}}-E,\end{array}$$

where SWE (m) is the snow water equivalent, P_r (m s⁻¹) is the rainfall rate, P_s (m s⁻¹) is the snowfall rate, M_r (m s⁻¹) is the meltwater outflow from the snowpack, and E is the sublimation from the snowpack.

A2.2 Energy balance equation

The energy balance of snowpack can be expressed as follows:

$$\begin{array}{}\text{(A9)}& {\displaystyle \frac{\mathrm{d}U}{\mathrm{d}t}}=Q_{\mathrm{sn}}+Q_{\mathrm{li}}+Q_{\mathrm{p}}+Q_{\mathrm{g}}-Q_{\mathrm{le}}+Q_h+Q_e-Q_{\mathrm{m}},\end{array}$$

where Q_sn (W m⁻²) is the net shortwave radiation, Q_li (W m⁻²) is the incoming longwave radiation, Q_p (W m⁻²) is the advected heat from precipitation, Q_g (W m⁻²) is the ground heat flux, Q_le (W m⁻²) is the outgoing longwave radiation, Q_h (W m⁻²) is the sensible heat flux, Q_e (W m⁻²) is the latent heat flux due to sublimation and condensation, and Q_m (W m⁻²) is the advected heat removed by meltwater.

Equations (8) and (9) form a coupled set of first-order, nonlinear ordinary differential equations. The Euler predictor–corrector approach was employed in the UEB model to solve the initial value problems of these equations (Tarboton and Luce, 1996).

A3 Albedo calculation

A3.1 Ground albedo

Instead of the constant bare soil albedo in the original UEB model, the bare soil albedo is expressed as a decreasing linear function of soil moisture in STEMMUS-UEB.

$$\begin{array}{}\text{(A10)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{g},\mathrm{v}}={\mathit{\alpha }}_{\mathrm{sat}}+\min {\mathit{\alpha }}_{\mathrm{sat}},\max [\left(\mathrm{0.11}-\mathrm{0.4}\mathit{\theta }\right),\mathrm{0}],\text{(A11)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{g},\mathrm{ir}}=\mathrm{2}{\mathit{\alpha }}_{\mathrm{g},\mathrm{v}},\end{array}$$

where α_g,v and α_g,ir are the bare soil and ground albedo for the visible and infrared band, respectively. α_sat is the saturated soil albedo, depending on local soil color. θ is the surface volumetric soil moisture.

A3.2 Vegetation albedo

The calculation of vegetation albedo is developed to capture the essential features of a two-stream approximation model using an asymptotic equation. It approaches the underlying surface albedo α_g,λ or the thick canopy albedo α_c,λ when the L_SAI is close to zero or infinity.

$$\begin{array}{}\text{(A12)}& {\displaystyle }\begin{array}{rl}{\mathit{\alpha }}_{\mathrm{Veg},\mathrm{b},\mathit{\lambda }}& ={\mathit{\alpha }}_{\mathrm{c},\mathit{\lambda }}\left[\mathrm{1}-\exp \left(-{\displaystyle \frac{{\mathit{\omega }}_{\mathit{\lambda }}\mathit{\beta }{L}_{\mathrm{SAI}}}{\mathit{\mu }{\mathit{\alpha }}_{c,\mathit{\lambda }}}}\right)\right]\\ & +{\mathit{\alpha }}_{\mathrm{g},\mathit{\lambda }}\exp [-\left(\mathrm{1}+{\displaystyle \frac{\mathrm{0.5}}{\mathit{\mu }}}\right)L_{\mathrm{SAI}}],\end{array}\text{(A13)}& {\displaystyle }\begin{array}{rl}{\mathit{\alpha }}_{\mathrm{Veg},\mathrm{d},\mathit{\lambda }}& ={\mathit{\alpha }}_{c,\mathit{\lambda }}\left[\mathrm{1}-\exp \left(-{\displaystyle \frac{\mathrm{2}{\mathit{\omega }}_{\mathit{\lambda }}\mathit{\beta }{L}_{\mathrm{SAI}}}{{\mathit{\alpha }}_{\mathrm{c},\mathit{\lambda }}}}\right)\right]\\ & +{\mathit{\alpha }}_{\mathrm{g},\mathit{\lambda }}\exp [-\mathrm{2}{L}_{\mathrm{SAI}}],\end{array}\end{array}$$

where subscripts Veg, b, d, c, g, and λ represent vegetation, direct beam, diffuse radiation, thick canopy, ground, and spectrum bands of either visible or infrared bands. μ is the cosine of solar zenith angle; ω_λ is the single-scattering albedo, amounting to 0.15 for the visible band and 0.85 for the infrared band, respectively; β is assigned as 0.5; L_SAI is the sum of leaf area index LAI and stem area index (SAI); and α_c,λ is the thick canopy albedo, which is dependent on vegetation type.

The bulk snow-free surface albedo, averaged between bare-ground albedo and vegetation albedo, is written as follows:

$$\begin{array}{}\text{(A14)}& {\mathit{\alpha }}_{\mathit{\eta },\mathit{\lambda }}={\mathit{\alpha }}_{\mathrm{Veg},\mathit{\lambda }}{f}_{\mathrm{Veg}}+{\mathit{\alpha }}_{\mathrm{g},\mathit{\lambda }}(\mathrm{1}-f_{\mathrm{Veg}}),\end{array}$$

where α_η,λ is the averaged bulk snow-free surface albedo and f_Veg is the fraction of vegetation cover.

A3.3 Snow albedo

According to Dickinson et al. (1993), snow albedo can be expressed as a function of snow surface age and solar illumination angle. The snow surface age, which is dependent on snow surface temperature and snowfall, is updated with each time step in UEB. Visible and near-infrared bands are separately treated when calculating reflectance and are further averaged as the albedo with modifications of illumination angle and snow age. The reflectance in the visible and near-infrared bands can be written as follows:

$$\begin{array}{}\text{(A15)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{vd}}=\left(\mathrm{1}-C_{\mathrm{v}}{S}_{\mathrm{age}}\right){\mathit{\alpha }}_{\mathrm{vo}},\text{(A16)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{ird}}=\left(\mathrm{1}-C_{\mathrm{ir}}{S}_{\mathrm{age}}\right){\mathit{\alpha }}_{\mathrm{iro}},\end{array}$$

where α_vd and α_ird represent diffuse reflectance in the visible and near-infrared bands, respectively. C_v (=0.2) and C_ir (=0.5) are parameters that quantify the sensitivity of the visible and infrared band albedo to snow surface aging (grain size growth), and α_vo (=0.85) and α_iro (=0.65) are fresh snow reflectance in visible and infrared bands, respectively. S_age is a function to account for aging of the snow surface and is given by

$$\begin{array}{}\text{(A17)}& S_{\mathrm{age}}={\displaystyle \frac{\mathit{\tau }}{\mathrm{1}+\mathit{\tau }}},\end{array}$$

where τ is the non-dimensional snow surface age that is incremented at each time step by the quantity designed to emulate the effect of the growth of surface grain sizes.

$$\begin{array}{}\text{(A18)}& \mathrm{\Delta }\mathit{\tau }={\displaystyle \frac{r_{\mathrm{1}}+r_{\mathrm{2}}+r_{\mathrm{3}}}{{\mathit{\tau }}_o}}\mathrm{\Delta }t,\end{array}$$

where Δt is the time step in seconds with τ_o=10⁶ s. r₁ is the parameter to represent the effect of grain growth due to vapor diffusion and is dependent on snow surface temperature:

$$\begin{array}{}\text{(A19)}& r_{\mathrm{1}}=\exp \left[\mathrm{5000}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{273.16}}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{s}}}}\right)\right].\end{array}$$

r₂ describes the additional effect near and at the freezing point due to melt and refreeze:

$$\begin{array}{}\text{(A20)}& r_{\mathrm{2}}=\min (r_{\mathrm{1}}^{\mathrm{10}},\mathrm{1}).\end{array}$$

r₃=0.03 (0.01 in Antarctica) represents the effect of dirt and soot.

The reflectance of radiation with illumination angle (measured relative to the surface normal) is computed as follows:

$$\begin{array}{}\text{(A21)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{v}}={\mathit{\alpha }}_{\mathrm{vd}}+\mathrm{0.4}f\left(\mathit{\phi }\right)(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{vd}}),\text{(A22)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{ir}}={\mathit{\alpha }}_{\mathrm{ird}}+\mathrm{0.4}f\left(\mathit{\phi }\right)(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{ird}}),\end{array}$$

where

$$f\left(\mathit{\phi }\right)=\left\{\begin{array}{ll}\frac{\mathrm{1}}{b}\left[\frac{b+\mathrm{1}}{\mathrm{1}+\mathrm{2}b\cos \left(\mathit{\phi }\right)}-\mathrm{1}\right],& \text{for }\cos \left(\mathit{\phi }\right)<\mathrm{0.5}\\ \mathrm{0},& \text{otherwise},\end{array}\right.$$

where b is a parameter set at 2 as in Dickinson et al. (1993).

When the snowpack is shallow (depth $z<h=\mathrm{0}$.01 m), the albedo is calculated by interpolating between the snow albedo and bare-ground albedo with the exponential term approximating the exponential extinction of radiation penetration of snow.

$$\begin{array}{}\text{(A23)}& A_{\mathrm{v}/\mathrm{ir}}=r{\mathit{\alpha }}_{\mathrm{g},\mathrm{v}/\mathrm{ir}}+(\mathrm{1}-r){\mathit{\alpha }}_{\mathrm{v}/\mathrm{ir}},\end{array}$$

where $r=\left(\mathrm{1}-\frac{z}{h}\right)e^{-z/\mathrm{2}h}$.

Appendix B

B1 Snow water equivalent

STEMMUS-UEB can reproduce the dynamics of snow water equivalent (Fig. B1). The discrepancies mainly happened under conditions with lower snow water equivalent. These intermittent shallow snowpack processes are difficult to capture well due to the drifting snow effect and temporal and complex ground heat conditions, and they require both high-quality observations and advanced snowpack models.

B2 Daily surface evaporation

Compared to the observations, surface evaporation was underestimated by the model with no snow module during the snowfall periods (Fig. B2). Models with snow module, however, produced a generally good agreement but with overestimations and underestimations, which corresponds to the mismatches in the snow water equivalent results. When the snow water equivalent is overestimated, snowpack sublimation and surface evaporation were overestimated.

Compared to the model without the snow module, the model with the snow module produced a better correlation with the measured daily surface evaporation (Fig. B3). Surface evaporation was underestimated by the model without the snow module and slightly overestimated by the model with snow module.

B3 Soil moisture and temperature

Models both with and without the snow module can reproduce the soil moisture dynamics in terms of their response to precipitation events (Fig. B4). Soil moisture was underestimated by the model without the snow module due to the lower amount of incoming water flux. Such underestimation was damped as the soil depth increases. Models with the snow module gain more incoming water (snowmelt water), and thus the underestimation of soil moisture was alleviated.

The dynamics of soil temperature were reproduced well by models both with and without the snow module (Fig. B5). There is no significant difference between soil temperature simulations of models with and without the snow module.

B4 Snow cover properties and albedo

There is a good correlation between the snow depth and surface albedo (Fig. B6). Figure B7 shows that surface albedo variations correspond well to the dynamics of the snow cover properties. This demonstrated that surface albedo is a reliable indicator to identify the presence of the snowpack and its influencing periods. Three example periods were selected to illustrate the validity of using the indirect method (albedo variation and ancillary meteorological data, i.e., air temperature, and precipitation) to define the presence and lasting time of the snowpack. Results indicated that the snowpack duration was successfully characterized using the indirect method (results were shown in Table S6.4 in the Supplement).

Figure B1 Time series of the observed and estimated snow water equivalent using the developed STEMMUS-UEB model.

Figure B2 Intercomparison of the observed and estimated surface evaporation using the model with and without the snow module.

Figure B3 Measured and estimated daily surface evaporation using the model with and without snow module (a and b, respectively).

Figure B4 Observed and estimated soil moisture at various soil layers using the model with and without the snow module.

Figure B5 Observed and estimated soil temperature at various soil layers using the model with and without the snow module.

Figure B6 Scatterplot of snow depth and albedo (Yakou station, 2014–2017).

Figure B7 Time series of the snow depth, snow water equivalent (SWE), and albedo (Yakou station).

Table B1 Notation.