STEMMUS-UEB v1.0.0: Integrated Modelling of Snowpack and Soil Mass and Energy Transfer with Three Levels of Soil Physical Process Complexities

Snowpack, as the indispensable component in cold regions, has a profound effect on the hydrology and surface energy conditions through its modification of the surface albedo, roughness, and insulating property. Although the modelling of the snowpack, soil water dynamics, and the coupling of the snowpack and underlying soil layer has been widely reported, the analysis of coupled liquid-vapor-air flow mechanisms 15 considering the snowpack effect was not yet investigated in detail. In this study, we incorporated the snowpack effect (Utah Energy Balance model, UEB) into a common modeling framework (Simultaneous Transfer of Energy, Mass, and Momentum in Unsaturated Soils with Freeze-Thaw, STEMMUS-FT), with various complexities of mass and energy transfer physics (from the basic coupled, to advanced coupled water and heat transfer, and further to the explicit consideration of airflow, termed BCD, ACD, and ACD-air, 20 respectively). We then utilized the in situ observations and numerical experiments to investigate the effect of snowpack on soil moisture and heat transfer with the above-mentioned model complexities. Results indicated that the abrupt increase of surface albedo after precipitation events can be only reproduced by models considering snowpack. The BCD model tended to overestimate the land surface latent heat flux. Such overestimations were largely reduced by ACD and ACD-air models. Compared with the simulation 25 considering snowpack, there is less surface latent heat flux from no-snow simulations due to the neglect of snow sublimation. With coupled models, the enhanced latent heat flux after precipitation events can be sourced from the surface ice sublimation, snow sublimation, and increased surface soil moisture, while the simple BCD model cannot provide the realistic partition of surface latent heat flux. The ACD model, with its physical consideration of vapor flow, thermal effect on water flow, and snowpack, can identify relative 30 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. https://doi.org/10.5194/gmd-2020-416 Preprint. Discussion started: 17 February 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
In cold regions, the snowpack has a profound effect on hydrology and surface energy through its modification 35 of the surface albedo, roughness and insulating property (Boone and Etchevers, 2001;Zhang, 2005).
Different than rainfall, precipitation water of snowfall enters the soil with a significant lag in time. However, a large and sudden outflow or runoff may be produced because of the snowmelt effect. The heat insulating effect of snow cover also provides a buffer layer to reduce the magnitude of the underlying subsurface temperature and thus markedly affect the thickness of the active layer in cold regions. The effect of snow 40 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 with the absorption of heat due to phase change (Zhang, 2005). Yi et al. (2015) investigated the seasonal snow cover effect on soil freezing/thawing process and its related carbon implications. Such studies mainly focus on the thermal effect 45 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).
Great amounts of modeling efforts have been made to better reproduce the snowpack characteristic and its effects. Initially, snowpack dynamics was expressed as a simple empirical function of temperature.
Nevertheless, these empirical relations have limited applications in complex climate conditions (Pimentel et 50 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; 55 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); iii) intermediate-complexity schemes with physics from the detailed schemes but with a limited amount 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 60 types of schemes are capable of representing the basic features of the snow cover over the 2-yr period but behaved differently on shorter timescales. Furthermore, 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, which was independent of model complexity (Etchevers et al., 2004) but affect directly snow simulations. SnowMIP2 evaluated thirty-three snowpack models across a wide range 90 facilitates our understanding of the hydrothermal dynamics of respective components in frozen soil medium (i.e., soil liquid water, water vapor, dry air, and ice) (Yu et al., 2018;Yu et al., 2020c).
In response to minimizing the potential model-comparison uncertainties rising from various model structures (Clark et al., 2015) 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 modelling framework (Yu et al., 2020c).

95
First, the 1-D Richards equation and heat conduction were deployed in STEMMUS-FT to describe the isothermal water flow and heat flow (termed BCD). In the BCD model, the interaction of soil water and heat transfer is only implicitly via the parameterization of heat capacity, thermal conductivity and the water phase change effect. For the advanced coupled water and heat transfer (ACD model), the water flow is affected by soil temperature regimes. The movement of water vapor, as the linkage between soil water and heat flow, is (see the more detailed model description in Zeng et al., 2011a, b;Zeng and Su, 2013;Yu et al., 2018;Yu et al., 2020c).

Snowpack module UEB
The Utah energy balance (UEB) snowpack model (Tarboton and Luce, 1996) is a single-layer physicallybased snow accumulation and melt model. 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 snowpack (i.e., solid, solid and liquid mixture, and 110 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 snow albedo (see Appendix A.3). The melt outflow is calculated using Darcy's law with the liquid fraction as inputs. The conservation of mass and energy, as presented in Appendix A.2, forms the physical basis of 115 UEB (Tarboton and Luce, 1996).
UEB is recognized as one simple yet physically-based snowmelt model, which can capture the first order snow process (e.g., diurnal variation of meltwater outflow rate, snow accumulation, and ablation, see a general overview of UEB model development and applications in Table S3). It requires little effort in parameter calibration and can be easily transportable and applicable to various locations (e.g., Gardiner et al., 120 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 as for example Tibetan Plateau. We thus selected the original parsimonious UEB (Tarboton and Luce, 1996) as the snow module to be coupled with soil module (STEMMUS-FT).

125
The coupled process between the snowpack model (UEB) and the soil water model (STEMMUS-FT) was illustrated in Figure 1. The one-way 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 130 snowpack energy and mass balance (Eq. A.8 & A.9, Subroutines: ALBEDO, PARTSNOW, PREDICORR), 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 ACD models; Diff_Moisture_Heat for BCD model).
STEMMUS-FT then solves the energy and mass balance equations of soil layers in one timestep. To highlight the effect of snowpack on the soil water and vapor transfer process, we constrained the soil surface energy 135 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 not permitted to be higher than https://doi.org/10.5194/gmd-2020-416 Preprint. Discussion started: 17 February 2021 c Author(s) 2021. CC BY 4.0 License. zero when there is snowpack. To ensure the numerical convergence, the adapted timestep strategy was used.
The half-hourly meteorological forcing measurements were linearly interpolated to the running timesteps (Subroutine Forcing_PARM). The precipitation rate (validated at 3-hour time intervals) was regarded 140 uniformly within the 3-hour duration (see Table S1 for detail). The general description of the main subroutines in STEMMUS-UEB, including the main functions, input/output, and its connection with other subroutines, was presented in Table 2 (linked with Table S1 and S2 for the description of model input parameters and outputs for this study, see the detailed general description in Tarboton and Luce, 1996;Zeng and Su, 2013;Yu et al., 2018).

Configurations of numerical experiments
On the basis of the aforementioned STEMMUS-UEB coupling framework, the various complexity of vadose zone physics was 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 termed BCD model. Second, the 150 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 ACD model. Lastly, on top of ACD model, the air pressure was independently considered as a state variable (therefore, the airflow) and termed ACD-air model. With the abovementioned model versions (STEMMUS-FT_Snow) and taking into account the no-snow scenarios (STEMMUS-FT_No-Snow), Table 3 lists the configurations of all six designed numerical experiments. The 155 model parameters used for all simulations for the tested experimental site are listed in Table S2.

Description of the Tested Experimental site
Maqu station, equipped with a catchment scale soil moisture and soil temperature (SMST) monitoring network and micro-meteorological observing system, is situated on the north-eastern edge of the Tibetan Plateau (Su et al., 2011;Dente et al., 2012;Zeng et al., 2016). According to the updated Köppen-Geiger 160 climate Classification System, it can be characterized as a cold climate with dry winter and warm summer.
The average annual air temperature is 1.2 ℃, and the mean air temperatures of the coldest month (January) and warmest month (July) are about -10.0 ℃ and 11.7 ℃, respectively. Alpine meadows (e.g., Cyperaceae and Gramineae), with a height varying from 5 cm to 15 cm throughout the growing season, are the dominant land cover in this region. The general soil types are sandy loam, silt loam and organic soil for the upper soil 165 layers (Dente et al., 2012;Zheng et al., 2015;Zhao et al., 2018). The soil texture and hydraulic properties were listed in Table S2 and how it used in STEMMUS-UEB is illustrated in Figure 1 and Table 2.  i.e., 5 cm, 10 cm, 20 cm, 40 cm, 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 https://doi.org/10.5194/gmd-2020-416 Preprint. Discussion started: 17 February 2021 c Author(s) 2021. CC BY 4.0 License. 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 down and upwelling solar and thermal radiation 175 (NR01-L, Campbell Scientific, Inc., USA), and liquid precipitation (T200B, Geonor, Inc., USA) are also deployed. The dataset from December 1, 2015 to March 15, 2016 was utilized in this study. An independent precipitation data (3-hour time interval) during the same testing period from an adjacent meteorological station was used as the mutual validation data.

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, was shown in Figure 2 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 185 average surface albedo, and lacked the response to the winter precipitation events ( Figure 2, Table 4). With the snow module, STEMMUS-UEB models can capture mostly 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 dynamics of albedo and the shallow snowfall events might be not well captured by the model (see the Sect. 4.1). Three model versions  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.

Soil Temperature and Moisture Dynamics
The observed spatial and temporal dynamics of soil temperature from five soil layers was used to verify the performance of different models (Fig. 3). The initial soil temperature state can be characterized as the warm 195 bottom and cool surface soil layers. The freezing front (indicated by the zero degree isothermal line, ZDIL) developed downwards rapidly till 70 th day after December 1, 2015, reaching its maximum depth. Then the freezing front stabilized as the offset effect of latent heat release (termed as zero-curtain effect). Such effect can sustain 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 200 thermal states. The isothermal lines (e.g., -2 o 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. model. The ACD models can well capture the propagation characteristic of freezing front 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 snow module, except at the surface soil layers (Table 4). Figure 4 shows the spatial and temporal dynamics of observed and simulated soil water content in the liquid 210 phase (SWCL). The SWCL of active layers is largely dependent on the soil freezing/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 ( < 0.06 3 −3 ) 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 215 to enlarge from BCD-Snow simulations due to its early freezing and thawing of soil (Fig. 4b). The position of the dry zone occurred earlier as the early reaching of the stabilization period by BCD-Snow model (Fig.   3b). For the ACD models, the position and development of initial wet zone, rewetting zone and the dry zone is similar to that from the observations, indicating the soil moisture dynamics can be well captured by the ACD models. Compared to STEMMUS-FT_Snow model, there was no observable difference in the SWCL 220 dynamics at deeper soil layers from STEMMUS-FT_No-Snow simulations. The surface SWCL was found affected from STEMMUS-FT_Snow simulations (Table 4). Figure 5 shows the comparison of time series of observed and model simulated surface cumulative latent heat flux using three models with/without consideration of snow module. Considerable overestimation of 225 latent heat flux was produced by the BCD-Snow model, with 121.79% more than observed. Such overestimations were largely reduced by ACD and ACD-air models. There is a slight underestimation of cumulative latent heat flux by ACD-Snow and ACD-air-Snow models, with -8.33% and -7.05%, respectively.

Surface Latent Heat Flux
Compared with STEMMUS-FT_Snow simulations, there is less latent heat flux produced by STEMMUS-FT_No-snow simulations. It is mainly due to the sublimation of snow cover, which cannot be simulated by 230 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.

Liquid/vapor fluxes
To further elaborate the effect of snowpack on LE, we presented the diurnal variations of LE and its 235 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 Figure 6 & 7, i.e., the liquid water flux driven by temperature qLT, matric potential qLh and air pressure qLa, water vapor flux driven by temperature qVT, matric potential qVh and air pressure qVa. https://doi.org/10.5194/gmd-2020-416 Preprint. Discussion started: 17 February 2021 c Author(s) 2021. CC BY 4.0 License.

240
Diurnal dynamics of observed and simulated latent heat flux during the rapid freezing period with the occurrence of precipitation events, from 10 th to 14 th Days after Dec. 1. 2015, is shown as Fig. 6a, d &g.
Compared to the observations, the diurnal variations of latent heat flux was 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 245 overestimation of LE was found at 10 th and 11 th day after December 1 due to relatively high surface soil moisture simulations (Fig. S1b). A certain amount of enhanced surface evaporation was produced shortly after precipitation, which is most probably due to the snow sublimation, which presents in the model simulations while not intuitively in observations. The mismatch in the LE enhancement after precipitation events can be attributed to that the partition process of precipitation into various components (rainfall, 250 snowfall, canopy interception) might not be well captured by the model. Such a response to the winter precipitation events was absent from the BCD-No-Snow simulations.
The overestimation of LE was reduced by ACD and ACD-air models ( Fig. 6d & g). Compared to ACD-Snow model simulations, ACD-No-snow model produced a stronger diurnal variation of LE after the precipitation, and is more approaching to the measured LE. The lower diurnal variation of LE for ACD-Snow model can 255 be ascribed to the lower surface SWCL (see Fig. S1d & g). For 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 (took 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.

260
During the freezing period, the soil water vapor, instead of liquid water flux, dominated the surface mass transfer process. Missing the description of the vapor diffusion process hindered the BCD models to realistically depict the decomposition of surface mass transfer dynamics. From the ACD model, there is a visible diurnal variation of thermal vapor flux qVT . It can be clearly identified that the isothermal vapor flux qVh contributed to most of the mass transfer during the freezing period. After 265 winter precipitation during the nighttime, there is a certain amount of isothermal vapor flux driven by the downward matric potential gradient. The reason is that the precipitation water immediately freezes on the soil surface. It should be noted that the sum of water/vapor fluxes at 0.1cm soil layer cannot balance the surface evaporation especially after the precipitation events (Fig. 6e). We assumed and attributed it to the surface ice sublimation process. Precipitation water was frozen on the soil surface, and only vapor exists in 270 the topsoil layers. Sublimation of surface ice may contribute to the gaps between liquid/vapor fluxes and LE (Yu et al., 2018). As more precipitation water was frozen on the soil surface from ACD-No-Snow model, the difference between the sum of water/vapor fluxes at 0.1cm soil layer and the surface evaporative water enlarged compared to ACD-Snow simulations (Fig. 6f). Thermal liquid water flux appears negligible to the https://doi.org/10.5194/gmd-2020-416 Preprint. Discussion started: 17 February 2021 c Author(s) 2021. CC BY 4.0 License. total mass flux during the whole simulation period. There is no significant difference recognized in the mass 275 transfer between the ACD-air and ACD during the freezing period.

Thawing period
During the thawing period, the diurnal variations of LE were well simulated by the models. There are some For ACD-Snow model, the precipitation induced evaporation enhancement was lagged compared to ACD-No-Snow model. This enhanced evaporation can be attributed to the snow sublimation and increased surface soil moisture content. The similar lag behavior of precipitation-enhanced evaporation was produced by the 295 ACD-air-Snow models.
During the thawing period, a certain amount of upward liquid water flux was produced by BCD model, supplying the water to evaporate into the atmosphere (Fig. 7b &c). Compared to the isothermal liquid flux, the thermal liquid flux was negligible to the total mass flux. For the ACD model, the diurnal variation of thermal vapor flux qVT was enhanced after precipitation, 300 producing a larger amount of upward/downward vapor flux during the night/day time (e.g., Fig. 7e). As the surface soil is relatively dry, the isothermal vapor flux contributes nearly all of the mass flux during the thawing period. Driven by the large downward matric potential gradient, a large amount of isothermal water vapor flux, accompanied by downward liquid water flux, can be found after the nighttime precipitation event.
This precipitation induced isothermal liquid/vapor flux was lagged and less intense from ACD-Snow model

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

Uncertainties of surface albedo simulations and Limitations
After winter precipitation, land surface albedo increases considerably (Fig. 2), indicating the formation of the snowpack. While such snowfall events were isolated with small magnitude, which are difficult to be well captured. Such difficulties can be partially attributed to the uncertainties in the representativity of 320 precipitation measurements. Due to the spatial variability of precipitation, the accurate observation of winter precipitation is proved to be a challenge, especially during windy winters (Barrere et al., 2017;Pan et al., 2017). On the other hand, the temporal resolution of precipitation measurements adopted in this study is relatively coarse (3-hour). 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 325 of precipitation measurements, the model may produce a time shift of snowfall events or even the malidentification of snowfall taken into account the effect of air temperature. 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 330 needs to be adapted to the specific site, especially regarding the shallow snow situations (Ueno et al., 2007;Ueno et al., 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 to capture the surface albedo variations with precipitation ( Fig. 2) and can be seen acceptable to conduct the analysis of snow cover effects in such a harsh condition. This amount of incoming water increased the evaporation after precipitation (Fig. 8 d). The other source for the enhanced evaporation flux after precipitation is snow sublimation, which is absent from the model without snow module. Sublimation, although not easy to observe, occurs readily under certain weather conditions (e.g., with freezing temperatures, enough energy). It can be further sped up at regions with low relative 345 humidity, air pressure and dry winds. Such amount of sublimation has been reported important from the perspective of climate and hydrology (e.g., Strasser et al., 2008;Jambon-Puillet et al., 2018), especially at high altitude regions with the low air pressure. During the freezing period, the evaporation enhancement can be also sourced from the sublimation of surface ice. The amount of ice sublimation appeared to be decreased during the freezing period as the presence of transient snowpack (Fig. 6). This is consistent with the results

350
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. Their results indicated that 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/temperature states, air temperature and the time and magnitude of 355 precipitation events. Under the conditions of the low pre-precipitation SWCL under the freezing soil temperature (e.g., Fig. 6e, 11 th vs. 12 th Days after 1 December), the precipitation falls down on the surface as snowfall and rainfall (mostly freezes as ice). The sublimation from surface ice can contribute to most of the total mass transfer (e.g., Fig. 6e, 11 th Days after 1 December). If the soil temperature rises above the freezing temperature, there will be no sublimation of surface ice contributing to the enhanced evaporation.   (Fig. 7).

380
Rendering from the aim to investigate the hydrothermal effect of the snowpack on the underlying soil system,

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.
where (m 3 m -3 ) is the volumetric water content; q (kg m -2 s -1 ) is the water flux; z (m) is the vertical direction coordinate (positive upwards); S (s -1 ) is the sink term for root water uptake; ρL (kg m −3 ) is the soil liquid 440 water density; K (m s -1 ) is the soil hydraulic conductivity; (m) is the soil water potential; t (s) is the time.
The heat conservation equation, considering the latent heat due to water phase change, can be expressed as: where Csoil (J kg −1 °C −1 ) is the specific heat capacity of bulk soil; T (°C) is the soil temperature; ρi (kg m −3 ) is the density of soil ice; Lf (J kg −1 ) is the latent heat of fusion; θi (m 3 m −3 ) is the soil ice volumetric water content. λeff (W m −1 °C −1 ) is the effective thermal conductivity of the soil; 445

A.1.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: where ρV and ρi (kg m −3 ) are the density of water vapor and ice, respectively; θL and θV (m 3 m −3 ) are the volumetric water content (liquid and vapor, respectively); qL and qV (kg m −2 s −1 ) are the soil water fluxes of liquid water and water vapor (positive upwards), respectively. KLh (m s −1 ) and KLT (m 2 s −1 °C −1 ) are the isothermal and thermal hydraulic conductivities, respectively. DVh (kg m -2 s -1 ) is the isothermal vapor conductivity; and DVT (kg m -1 s -1 °C -1 ) is the thermal vapor diffusion coefficient.

455
On the basis of De Vries (1958) and Hansson et al. (2004)'s work, the heat transport function in frozen soils, considering the fully coupled water and heat transport physics, can be expressed as: temperature Tr; W (J kg −1 ) is the differential heat of wetting (the amount of heat released when a small amount of free water is added to the soil matrix).

A.1.3 Coupled mass and heat physics with air flow
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 465 corresponding governing equations of liquid water flow, vapor flow and air flow.
where ℎ , , and (kg m -2 s -1 ) are the liquid water fluxes driven by the gradient of matric potential , temperature , and air pressure , respectively. ℎ , , and (kg m -2 s -1 ) are the water vapor fluxes driven by the gradient of matric potential , temperature , and air pressure , respectively. Pg (Pa) is the mixed pore-air pressure.

�(
where ε is the porosity; Sa (=1-SL) is the degree of air saturation in the soil; SL (=θL/ε) is the degree of saturation in the soil; Hc is Henry's constant; De (m 2 s -1 ) is the molecular diffusivity of water vapor in soil; Kg (m 2 ) is the intrinsic air permeability; µa ( kg m -2 s -1 ) is the air viscosity; θa (=θV) is the volumetric fraction of dry air in the soil; and DVg (m 2 s -1 ) is the gas phase longitudinal dispersion coefficient.

A.2.1 Mass balance equation
The increase or decrease of snow water equivalence with time equals the difference of income and outgoing water flux:

490
(m/s) is the meltwater outflow from the snowpack; and is the sublimation from the snowpack.

A.2.2 Energy balance equation
The energy balance of snowpack can be expressed as: where (W/m 2 ) is the net shortwave radiation; (W/m 2 ) is the incoming longwave radiation; (W/m 2 ) is the advected heat from precipitation; (W/m 2 ) is the ground heat flux; (W/m 2 ) is the outgoing 495 longwave radiation; ℎ (W/m 2 ) is the sensible heat flux; (W/m 2 ) is the latent heat flux due to sublimation/condensation; and (W/m 2 ) is the advected heat removed by meltwater.
Equations (8) and (9) form a coupled set of first order, nonlinear ordinary differential equations. Euler predictor-corrector approach was employed in UEB model to solve the initial value problems of these equations (Tarboton and Luce, 1996). where , and , are the bare soil/ground albedo for the visible and infrared band, respectively. is 505 the saturated soil albedo, depending on local soil color. is the surface volumetric soil moisture.

A.3.2 Vegetation albedo
The calculation of vegetation albedo is developed to capture the essential features of a two-stream approximation model using asymptotic equation. It approaches the underlying surface albedo , or the thick canopy albedo , when the is close to zero or infinity. where subscripts , , , , 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, 0.15 for visible and 0.85 for infrared band, respectively; is assigned as 0.5; is the sum of leaf area index LAI and stem area index SAI; , is the thick canopy albedo dependent on vegetation types. The bulk snow-free surface albedo, averaged between bare ground albedo and vegetation albedo, then is written as: where , is the averaged bulk snow-free surface albedo; is the fraction of vegetation cover.

A.3.3 Snow albedo
According to Dickinson et al. (1993), snow albedo can be expressed as a function of snow surface age and 520 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, which 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: where and represent diffuse reflectance in the visible and near infrared bands, respectively. (= 525 0.2) and (=0.5) are parameters that quantify the sensitivity of the visible and infrared band albedo to snow surface aging (grain size growth), (=0.85) and (=0.65) are fresh snow reflectance in visible and infrared bands, respectively.
is a function to account for aging of the snow surface, and is given by: where τ is the non-dimensional snow surface age that is incremented at each time step by the quantity 530 designed to emulate the effect of the growth of surface grain sizes.
where ∆ is the time step in seconds with = 10 6 s. r1 is the parameter to represent the effect of grain growth due to vapor diffusion, and is dependent on snow surface temperature: r2 describes the additional effect near and at the freezing point due to melt and refreeze: r3=0.03 (0.01 in Antarctica) represents the effect of dirt and soot.

535
The reflectance of radiation with illumination angle (measured relative to the surface normal) is computed as: where b is a parameter set at 2 as Dickinson et al. (1993).

810
---> means the relevant subroutines which are incoming to the current one, --> means the relevant subroutines for which the current subroutine is output to; agesn 2 and ALBEDO 2 , means the use of subroutines agesn and ALBEDO after solving the snowpack energy and mass conservation equations, to update the snow age and albedo.
, where , �, are the measured and model simulated values of the selected variable (snow albedo, LE, soil temperature/moisture); � is the mean values of the measurements of the selected variable (snow albedo, LE, soil temperature/moisture); n is the number of data points.
The correlation is all significant at the 0.01 level, except for "-", which indicates that the correlation is not significant.  and are soil liquid water and ice content; is soil hydraulic conductivity; is thermal conductivity. , , 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).

830
UEB, Utah Energy Balance module. Precip, Ta, HRa, Rn, and u are the meteorological inputs (precipitation, air temperature, relative humidity, radiation and wind speed). Model subroutines are in red fonts.