Evaluating the utah energy balance (UEB) snow model in the noah land-surface model

. Noah (version 2.7.1), the community land-surface model (LSM) of National Centers for Environmental Predictions-National Center for Atmospheric Research (NCEP-NCAR), which is widely used to describe the land-surface processes either in stand-alone or in coupled land-atmospheric model systems, is recognized to underestimate snow–water equivalent (SWE). Noah’s SWE bias can be attributed to its simple snow sub-model, which does not ef-fectively describe the physical processes during snow accumulation and melt period. To improve SWE simulation in the Noah LSM, the Utah Energy Balance (UEB) snow model is implemented in Noah to test alternate snow surface temperature and snowmelt outﬂow schemes. Snow surface temperature was estimated using the force-restore method and snowmelt event is regulated by accounting for the internal energy of the snowpack. The modiﬁed Noah’s SWE simulations are compared with the SWE observed at Cali-fornia’s NRCS SNOTEL stations for 7 water years: 2002– 2008, while the model’s snow surface temperature is veri-ﬁed with observed surface-temperature data at an observation site in Utah. The experiments show that


Introduction
The Noah LSM, a moderately complex community model, is widely used in weather and regional-climate models and is the operational land-surface scheme for National Centers for Environmental Predictions-National Center for Atmospheric Research (NCEP-NCAR) (Chen et al., 1996;Chen and Dudhia, 2001;Ek et al., 2003;Leung et al., 2005Leung et al., , 2006Jin and Miller, 2007;Jiang et al., 2008). It is also used in land-data assimilation systems such as the North America Land Data Assimilation System , the Land Information System (Peters-Lidard et al., 2007), and HRLDAS (Chen et al., 2007). The model has been advanced numerous times to accurately predict warm-and cold-season processes, which resulted in various versions of the Noah LSM. Noah version 2.7.1 (studied here) includes major modifications in snow processes by Koren et al. (1999) and Ek et al. (2003). However, the model has been noted for substantially underestimating snow water equivalent (SWE) (Jin et al., 1999a;Pan et al., 2003;Sheffield et al., 2003;Mitchell et al., 2004;Jin and Miller, 2007;Slater et al., 2001;Livneh et al., 2009;Wang et al., 2010;Barlage et al., 2010;Niu et al., 2011;Yang et al., 2011) by simulating less amounts of snow during peak winter season as well as melting the snow earlier in the spring. Physical processes that influence the model's prediction of SWE are primarily the (1) representation of snowpack and underlying half of top soil layer as a single bulk layer, (2) snow albedo parameterization, (3) lack of snow-water retention and refreeze, and (4) snowmelt based on residual energy from the surface energy balance, (Livneh et al., 2009;Barlage et al., 2010;Niu et al., 2011). First two processes control the availability of energy in the snowpack while the Published by Copernicus Publications on behalf of the European Geosciences Union.
last two processes regulate snowmelt. Single layer snowpack combined with vegetation and underlying soil layer underestimates ground heat flux followed by overestimation of snow surface energy . Further energy is added at the snow surface due to the model's snow albedo parameterization which does not consider high reflectivity of fresh snow and snow aging (Livneh et al., 2009;Barlage et al., 2010). The residual energy in the snowpack is directly used to melt snow instead of using some available energy to warm the snowpack. The warm snowpack does not retain any liquid water which can refreeze at night and thaw during day before finally draining from the snowpack. Livneh et al. (2009) and Barlage et al. (2010) suggested that inclusion of snow-aging processes in the snow-albedo decay scheme can reduce Noah's SWE estimation bias. Livneh et al. (2009) have also implemented snow-water retention algorithm which also improved the model's SWE prediction. The limitation of Noah's single layer snowpack has been considered by Niu et al. (2011). Since, Noah computes a single temperature for the entire snowpack disregarding the temperature variation within the snow depth, Niu et al. (2011) replaced the model's single-layer snowpack representation with multiple layers to explicitly capture the non-linear temperature gradient of the snowpack. Recognizing the difference in snow surface and bottom temperature improves prediction accuracy of snow surface temperature, surface fluxes and ground heat flux. Therefore, most complex snow models (e.g., SNTHERM -SNow THERmal Model, Jordan, 1991; CLM -Community Land Model, Dai et al., 2003; SAST -Snow-Atmosphere-soil Transfer, Jin et al., 1999b) also apply finite-difference models to simulate snowpack temperatures. In addition, snow accumulation and ablation processes are also affected by land covers, which are addressed by several research groups (e.g., Mahat and Tarboton, 2012). To enhance the model, Niu et al. (2011) has also tested the model by separately computing temperature and heat fluxes from the canopy layer and included frozen soil scheme to improve soil permeability. Wang et al. (2010) have shown that Noah SWE simulation can be improved by considering the vegetation shading effect, under-canopy resistance, and roughness length adjustment in boreal forests and other grasslands.
In this study, we address the problem of Noah's SWE bias and early snowmelt by implementing the snow surface temperatures and snowmelt processes of the Utah Energy Balance (UEB) model Tarboton, 1994;Tarboton and Luce, 1996;Luce and Tarboton, 2001) in the Noah LSM. Similar to the Noah model, UEB simulates snowpack as a single layer but applies the force-restore method, which, unlike the finite-difference methods, implicitly represents temperature profiles within the snowpack.
Single-layer snow models like Noah LSM are less complex to apply as it requires a small number of input variables compared to complex models. More sophisticated models (ex, SNTHERM, Jordan, 1991;CLM, Dai et al., 2003;SAST, Jin et al., 1999b) simulate snowpack in multilayer and changes in snow properties within each layers (Anderson, 1976;Colbeck and Anderson, 1982;Jordan, 1991;Arons and Colbeck, 1995) are estimated, which are useful in some applications. However, the only piece of information required for climate study and hydrologic prediction is the snowskin temperature, because the temperature gradient between the snow surface and atmosphere drives the turbulent fluxes (Dickinson et al., 1993;Luce and Tarboton, 2001).
In this study, we have chosen the force-restore-based UEB snow model as the target in the Noah LSM for benefiting from its effectiveness in snow surface temperature estimation and snowmelt process, while keeping the Noah landsurface model's snow sub-model as a single layer model. The primary objective of this research is to evaluate Noah's SWE prediction with the alternate snow surface temperature and melt model formulations. While studies have identified factors that improve snow surface fluxes followed by enhancement in SWE simulations, there has not been any study that has considered substitute modeling of snow surface temperature and melt processes by preserving Noah's single layer snow model. This research is an attempt to answer the question -can an alternate single layer snow model reduce Noah's SWE bias?
Noah LSM's (version 2.7.1) approach for simulating snow temperature and snowmelt conditions will be described in the next subsection. The UEB snow model is discussed in Sect. 3. In Sect. 4, the study areas and forcing information are given. Results and discussion on the modified Noah model are given in Sect. 5. Finally, the conclusions are drawn in Sect. 6.

Current method for snow-temperature estimation
The Noah model, originally developed by Mahrt and Pan (1984) and Pan and Mahrt (1987), applies energy and water balance to simulate land-surface conditions. The model's physical representation has been enhanced numerous times and updated versions of the model are periodically published at NCAR website (http://www.ral.ucar.edu/ research/land/technology/lsm.php). The model is driven by seven input variables -precipitation, air temperature, surface pressure, wind speed, relative humidity, downward and upward shortwave radiation. This stand-alone, 1-D column version (version 2.7.1) has a multi-layer soil model but a simple canopy and snow model. If air temperature is less than 0 • C, precipitation falls as snowfall. Snow cover area fraction within a model grid is determined as a function of SWE using a generalized snow depletion curve. When snow is on the ground, the model considers a bulk snow-soil-canopy layer and computes a single surface temperature for the bulk layer at every time step. The other state variables in the Noah's  (Ek et al., 2003). For each location, maximum snow albedo is derived from a database developed based on the work of Robinson and Kukla (1985). The data set covers the area of 25 • N at 1 • × 1 • resolution. In the Noah LSM, snow surface temperature (T s ) for the entire snowpack is estimated in a two-step process and the computation processes are graphically represented in Fig. 1 during a snow season at a California snow measuring station. In the first step, an intermediate temperature (T 12 ) is computed by applying energy balance between the snowpack, top soil layer and the overlying air (Koren et al., 1999;Ek et al., 2003). Detail of the energy balance equation is given in the Appendix A. The method allows this temperature to rise above freezing temperature (T freeze ) even when the model grid is 100 % covered with snow as can be seen in Fig. 1b.
In the second step, the effective temperature for a model grid (T 1 ) is adjusted by accounting for the fractional snow cover (f sca ) in the ground as (Koren et al., 1999): Equation (1) describes that, when the ground is completely snow-covered, T 1 is essentially snow surface temperature which can be below or at freezing temperature. On the other hand, T 1 can be above freezing temperature when the ground is partially covered with snow. Figure 1c shows the calculated snow surface temperature at the end of the two-step process. One of the primary deficiencies of the control model is that the snow sub model is conditioned to initiate snowmelt whenever the temperature T 1 is at or above the freezing point ( Fig. 1d). But, snow surface at freezing point is not the single factor to cause snowmelt. Snowmelt initiates only after the entire snowpack is isothermal at 0 • C which is known as warming phase (Dingman, 1994). This phase is followed by a ripening phase when melt water can retain in the snowpack until it exceeds liquid water holding capacity of the snowpack. Therefore, to initiate snowmelt, it is crucial to know at which state the snowpack is and that can be determined by accounting for the internal energy of the snowpack (Tarboton and Luce, 1996). Further, the warming and ripening processes are not considered in Noah's snow sub-model and, therefore, the model overestimates the net energy which is entirely used to control snowmelt outflow rate (Livneh et al., 2009;Lundquist and Flint, 2006) resulting in a faster melt rate.

The UEB snow model
To overcome the deficiencies in Noah's snow model, snow surface temperature and snowmelt processes of the Utah Energy Balance (UEB) snow model are evaluated as an alternate method to the existing snow model. The UEB model, originally developed by Tarboton et al. (1994) and Tarboton and Luce (1996), is a physically based energy and mass balance model to simulate snow accumulation and snowmelt at a point location. Snowpack is defined in a single layer by three state variables -snow-water content, internal energy of the snowpack, and the dimensionless age of the snow surface. Input variables to the model are -air temperature, precipitation, wind speed, humidity and radiation (Tarboton and Luce, 1996). Tarboton and Luce (1996) assume neutral stability in the UEB model (Hellstrom, 2000). At every time step, the snow surface temperature is computed based on an energy balance between surface forcing and the snow surface capacity to conduct heat into or out of the snowpack (Mahat and Tarboton, 2013). Then snow surface temperature is applied to compute internal energy of the snowpack which in turn regulate snowmelt outflow. The model considers liquid water content to onset snowmelt. Since its development, the model has been tested and verified at different sites with additional efforts to enhance the model performance (Luce and Tarboton, 2001;You, 2004;Luce and Tarboton, 2010;Mahat and Tarboton, 2013). A detailed discussion of the UEB model can be found in Tarboton et al. (1994Tarboton et al. ( , 1995, Tarboton and Luce (1996), and Luce and Tarboton (2001), while a brief discussion of the model's physical processes pertinent to this paper is given below.

Snow surface temperature formulation
To compute snow surface temperature UEB model applies energy balance at the snow surface whereas Noah LSM computes net energy for the entire snowpack layer. In reality, the snow surface temperature is cooler (warmer) than the entire snowpack temperature at night (day) time and therefore results a non-linear temperature gradient within a snowpack layer. The UEB model approximates this temperature gradient by differentiating surface temperature (T s ) from the average snowpack temperature (T ). At every time step, snow-skin temperature (T s ) is numerically solved using the Newton-Raphson method by employing the following (Tarboton and Luce, 1996): where Q cs is the heat flux because of the temperature gradient within the snowpack and Q forcing is essentially the heat flux considering all of the energy components at the snow surface, Q snet is the net short-wave energy, Q lin is the incoming long-wave radiation, Q lout is the outgoing long-wave radiation from the snowpack, Q h is the sensible heat flux to/from the snow, Q le is the latent heat flux to/from the snow, and Q p is the energy advected by precipitation into the snow. The turbulent heat fluxes (Q h and Q le ) and the outgoing radiative flux (Q lout ) are functionally dependent on the surface temperature T s (Tarboton and Luce, 1996;You, 2004). Physically, snow temperature cannot be greater than 0 • C and, thus, the upper bound of T s is constrained to freezing temperature. Q cs in Eq.
(2) is derived from thermal-diffusion equation which describes how temperature changes with time along the depth of a layer and the equation is Tarboton, 2001, 2010;You, 2004): with the boundary condition for temperature as: where T is the temperature, t is time, z is the depth measured downward from the surface, c is the volumetric heat capacity, and λ is the heat conductivity, k is the thermal diffusivity of snow (= c λ ), T is the mean snow surface temperature, A is the amplitude of the diurnal snow-temperature wave at the surface, and ω is the angular velocity of the Earth's rotation (i.e., ω = 2π/24 radians h −1 ). An approximate solution of Eq. (5) for sinusoidal temperature fluctuation is (Carslaw and Jaeger, 1959): where d is the diurnal damping depth, and d = √ 2k/ω. The equation may be used to calculate the temperature gradient with depth which can be used with the heat conductivity (λ) to compute heat flux (Q c ) (Lin, 1980;Hu and Islam, 1995;Luce and Tarboton, 2001;Gao et al., 2008): Further, rearranging Eq. (6) and then differentiating with respect to time t can be represented as follows: Using Eqs. (8a) and (8b), Eq. (7) can be written as: At the surface boundary where z = 0, the heat flux at the surface (Q cs ) is: Equation (8c) is the basis for the force-restore method and with the finite-difference approximation for ∂T s /∂t in Eq. (9c) results in where t is the measurement time interval, and T slag1 is the surface temperature lagged by one time step, i.e., at t − t.
(2)-(10) has not been defined consistently in the literature (see examples of various definitions of T in Ren and Xue, 2004). In the UEB model, T is defined as the depth-average snowpack temperature, which is derived from two state variables: snow-water equivalent (W ) and internal energy (U ) of the snowpack. Internal energy U is defined as the energy to the melting point which means that, at 0 • C, ice possesses zero heat content Jin et al., 1999a). Internal energy has been also used as a prognostic variable by Lynch-Stieglitz (1994) and Jin et al. (1999b).
Snow temperature affects sublimation from the snow surface. UEB applies turbulent heat flux (Tarboton, 1994) while Noah LSM uses Penman equation (Wang et al., 2010) to compute sublimation.

Snowmelt formulation
In UEB model, whenever internal energy is positive, the snowpack attains sufficient energy to initiate snowmelt.The snowmelt outflow rate M r from ripened snow is simulated based on Male and  and is (Tarboton, 1994;Tarboton and Luce, 1996): where S is the relative saturation in excess of the liquid-water holding capacity, and K s is the snow-saturated hydraulic conductivity, which describes the water flux through the porous snowpack and is a function of snow density, porosity, and liquid-water holding capacity. The variation of K s with a saturated water content of natural snow is not clear (Iida et al., 2000) and, hence, K s is essentially a calibration parameter for each location (Tarboton and Luce, 1996). Different showed snowpack melts rapidly at the end of melting period with higher K s value. A constant K s value of 20 m h −1 for all sites was chosen which reasonably described the rate and timing of snowmelt during the accumulation and ablation period. The parameter S in Eq. (11) is derived from the following relationship:

S =
Liquid water volume-Capillary retention Pore volume-Capillary retention , where the value of variable S increases with increasing liquid water in the snowpack. Liquid water is the amount of water that can be retained in the snow pores against capillary forces, and consideration of capillary retention or liquidwater holding capacity can delay snowmelt during the ripening phase (Dingman, 1994). Amid the ripening phase, liquid water near the surface can refreeze with night-time cooling and thaw during day (Bargtsson, 1980). This refreeze and thaw cycle can continue for days if the liquid water does not exceed the water-holding capacity of the snowpack. During the day, this cycle might need several hours to warm up and resume melting again (Dingman, 1994). Snowmelt starts once the liquid water in the snowpack exceeds the waterholding capacity. Initially, snow melting is more uniform ("matrix flow" in porous media) but with increase in liquid water content and growth of snow grains, melt flow rate accelerates ("preferential flow"). Theoretical representation of preferential flow is difficult which can advance below freezing temperature, so snowmelt algorithm even in the sophisticated snow models (e.g., like SNTHERM; Jordan, 1991) is based on liquid water flow under isothermal conditions (Waldner et al., 2004). However, the parameter "liquid-water holding capacity" is difficult to measure from wet snow because, during the snowmelt metamorphism, snow can be supersaturated yet be below the liquid-water holding capacity due to the freeze-thaw cycle (Livneh et al., 2009). In various studies, the liquid water-holding capacity is quantified as 3-9 % of the volume of the snowpack (Denoth et al., 1984;Kattlemann, 1987;Kendra et al., 1994;and Albert and Krajeski, 1998). Jordan (1991) used 4 % of the pore volume in SNTHERM, Lynch-Stieglitz (1994) used 5.5 % height of the compacted snow layer, while Dingman (1994) suggested 6 % of the pore space as the liquid water-holding capacity. Following Jordan (1991) and Denoth (2003), Livneh et al. (2009) applied 4 % of the pore volume of the liquid-water holding capacity in the Noah model. Because the density of the snowpack is different for fresh snow compared to old snow, Livneh et al. (2009) showed that 4 % of the pore volume can range from approximately 2.5 % of SWE depth for old snow to approximately 10 % of SWE depth for fresh snow. Here, we have used 5 % of the total mass of the snowpack (liquid and ice) as the liquid water-holding capacity (Tarboton and Luce, 1996), and the fraction of liquid water L f is estimated as: where U is the internal energy of the snowpack, W is the snow-water equivalent, ρ w is the density of water, λ f is the heat of fusion of ice, and ρ w λ f W is the energy required to melt the entire snowpack at 0 • C. A sensitivity test of the liquid-water holding capacity is discussed in the results and discussion section.
3 Application of the method

Study area and input data
The Noah LSM requires seven input variables, and the forcing data from the North American Land Data Assimilation System (NLDAS-2) are used to drive the model. NL-DAS forcing data are at 1/8 • -grid resolution (approximately 12.5 km) and available at hourly time scale and discussed in detail by Cosgrove et al. (2003) and extensively validated by Luo et al. (2003) and Pinker et al. (2003). The model outputs are evaluated against ground-observation data at various SNOTEL (SNOw TELemetry) stations located in the Sierra Nevada Mountains of California as well as at a SNOTEL site in Utah, which is close to a snow-surface temperature data collection site. Brief descriptions of these study areas are given below: (

1) SNOTEL stations in California
Snow-water equivalent (SWE) observations at SNOTEL sites in California are used here to verify the control and outputs of the modified models. These ground observations are at clear vegetated land and provide seasonal variation of snowfall with reasonable accuracy but is subjected to instrumentation accuracy (Serreze et al., 1999). SNOTEL stations have also been used by others for Noah-related studies Jin and Miller, 2007;Livneh et al., 2009;Barlage et al., 2010;Pederson et al., 2010). These stations measure daily SWE, 2-m air temperature, and precipitation, as well as soil moisture and soil-temperature data; a more comprehensive discussion about these sites is given on the Natural Resources Conservation Service (NRCS) website. Here, the data used are quality-controlled as described by Serreze et al. (1999), Over the study period, the maximum SWE is observed at the highest-elevated station (SNOTEL Station #574) and was 2628 mm in 2006. This year was an El Niño year, which usually generates higher than normal SWE values in the Southwest (Jin and Miller, 2007).
Land-surface characteristics were derived from the grid cell within which SNOTEL stations are located. SNOTEL stations are installed in open spaces but we wanted to simulate Noah with grid-scale attributes to facilitate model application at large scale. Table 1 lists elevation and average annual rainfall over four SNOTEL stations: #356, #508, #463,  and atmosphere are measured every 30 min at 12 monitoring weather stations and at an eddy covariance tower. Weather stations distributed across the site record air temperature, humidity, snow depth (when present), soil moisture, and sundry other quantities, while the eddy covariance system records atmospheric fluxes. Snow-temperature data collection started from 2008 and for this work, water year 2009 data were used. More information about this study site can be found at http://danielforest.usu.edu/Maps.aspx.
Near the experimental site, a SNOTEL station (#1098) was also installed in July 2007 where daily precipitation, temperature, snow-water equivalent, snow surface temperature are recorded at hourly rate. The model was simulated using the station recorded precipitation and temperature. Other meteorological inputs for the model were used from the NL-DAS grid data where this SNOTEL station resides. Modelsimulated snow surface temperature data were compared with the site recorded of snow surface temperature.

Nash-Sutcliffe model efficiency coefficient
To assess the goodness-of-fit of a model, the Nash-Sutcliffe model efficiency (NSE) coefficient is widely used and is defined as follows: where S o is observed SWE, S m is modeled SWE, and S o is the mean of observed SWE during the total time period T . NSE can range from −∞ to 1. An efficiency of less than zero (i.e NSE < 0) denotes that the model is not a good predictor of the variable of interest (Krause et al., 2005) whereas an efficiency of 0 (i.e., NSE = 0) indicates that the model predictions are as accurate as the mean of the observed data. While, Gupta et al. (1999) suggested the NSE values larger than 0 for minimally "acceptable performance", in literature, various threshold values are used for model's "satisfactory performance" (see Table 2 of Moriasi et al., 2007). Because NSE values are not easily interpretable if sampling distribution is not given and users can only provide subjective interpretation (McCuen et al., 2006). In essence, the closer the model efficiency is to 1, the more accurately the model matches with observation. Here, NSE values greater than 0.7 (Gupta et al., 1999) and less than 0.5 (Moriasi et al., 2007) are considered as "good" and "unsatisfactory" model performance.

Root Mean Square Error
Root Mean Square Error (RMSE) is commonly used model evaluation index and is defined as: where n is the total number of observations. Values of RMSE can be 0 to ∞ and for a model that perfectly fits the observation, value of RMSE is 0. There is no generally accepted threshold RMSE value to evaluate model performance but large RMSE values indicate large model error. In this study, followed by Singh et al. (2004), half the standard deviation of observed data will be used as a threshold RMSE value for model assessment. In addition, the model enhancements will be evaluated using RMSE-observations standard deviation ratio (RSR) (Moriasi et al., 2007) which is the ratio of the RMSE and standard deviation of observed data. RSR is calculated as follows: Model run before precipitation bias correction is called "control", and the model run after precipitation bias correction is called "control-bias corr".
RSR can range from 0 to very large value. A value of 0 indicates 0 RMSE and perfectly fit model simulation. The benefit of using RSR is in its normalization factor that allows developing a model performance rating. Followed by Moriasi et al. (2007), model performance is rated as "very good" if RSR value is between 0 and 0.5. Model will be regarded as "good", "satisfactory", and "unsatisfactory" for RSR values between 0.5 to 0.6, 0.6 to 0.7, and greater than 0.7, respectively. Figure 3 shows a comparison between the Noah LSM simulated SWE (blue line; called "control" in Fig. 3) and the observed SWE at the four SNOTEL stations. In general, the control model underestimates the SWE, although the underestimation of the modeled SWE differs from year to year and from station to station. At some stations and during some years, simulated SWE compares relatively well to that of observed SWE than other locations and years. For example, at Station #356, the SWE during water years 2002, 2003, 2004, and 2008, were certainly less than observed but modestly captured (∼ 40-60 % of the maximum SWE at the ground). Elsewhere, simulated SWE is almost negligible, particularly at Station #463 (less than 20 % of the ground observed maximum SWE). The primary reason for the Noah land-surface model's negative bias in SWE estimation is because of imperfection in its current snow-physical processes, as discussed earlier. In addition, uncertainty in input can be quite substantial as well, especially in the mountainous environments. Precipitation, a primary input for quantifying snowfall, can be extremely variable in space and time in high-elevated areas.  discussed this issue and pointed out that precipitation data in the NLDAS system are based on the National Weather Service (NWS) precipitation gauges. These gauges, located mostly in valleys, are known to underestimate higherelevation precipitation . Pan et al. (2003) compared the NLDAS precipitation with SNOTEL precipitation from September 1996-September 1999 and found that SNOTEL precipitation is, on average, more than twice the amount of the NLDAS precipitation data. On the other hand, differences in other forcing data between NLDAS and those of stations were found to be insignificant .

Precipitation bias correction
Therefore, a simple precipitation bias-correction method was applied to NLDAS precipitation data, while no corrections were made to other NLDAS forcing data. Precipitation data were adjusted by first determining the ratio of total yearly winter precipitation (October-May) from SNO-TEL stations to that of NLDAS grid. Then, NLDAS precipitation was scaled by the corresponding ratio. This simple bias correction shows that, in general, NLDAS precipitation is less than that recorded at the studied SNOTEL sites ; thus, a substantial increase in SWE can be seen in years and stations where SWE was very poorly modeled (e.g., Station #463). There are a few years in which NLDAS precipitation data were more than the total precipitation recorded at the SNOTEL stations and, therefore, precipitation bias correction resulted in reduced simulated SWE when compared to the control model (e.g., water year 2004 at Station #356, water year 2008 at Station #539). However, the number of snow-covered days has not been affected significantly (red line in Fig. 3 and termed as "control-bias-corr") with the bias correction.
Model bias can also increase at sites where additional snowdrifts can result from wind or at sites with precipitation under-catch which is common problem at mountainous climate stations. A study by Gaudet and Cotton (1998) in Colorado mountain region found more than 20 % precipitation under-catch in climate stations. With additional bias correction model forecast skill can be further enhanced.
UEB model considers effect of wind drift by wind drift parameter but the Noah land-surface model does not incorporate any physical processes for the effect of wind drifts; consequently, the UEB's drift parameter was not included and no additional processing was not done for cases when accumulated SWE exceeded accumulated rainfall. From this point on, the simulated SWE after bias correction will be referred to as the "control" run of the Noah LSM and will be used for evaluation of the modified approach, which is termed as "Noah-T s " run.

SWE simulation
Simulation of the Noah LSM modified with the UEB snow model is compared at SNOTEL stations and is shown in Fig. 4. The modified Noah shows substantially improved SWE estimation in terms of increasing the amount of maximum SWE as well as delaying snowmelt. However, water year 2007 was a dry year (a moderate La Niño year) and snowfall was less than 7-year average (study period: water year 2002 to 2008). For this year, SWE simulation in the modified model is not significant compared to the control run because for shallow depth of snow (less than 0.1 m), the modified model applies all the available energy for snowmelt. While the modified model enhanced SWE simulation by using the Noah LSM, it also shows delayed SWE melting in few years, for example, in water year 2004 at Stations #463 and #508 and in water year 2008 at Station #508. This late melting can be partly explained by comparing the simulated SWE by the control-bias-corr model and the control model (Fig. 3) at Station #508. At this station, NDLAS precipitation was more than that observed in 2004 and 2008 and, therefore, after precipitation bias correction, the controlbias-corr model predicted less snow than the control model. In general, forcing data from NLDAS, other than precipitation, are well validated but, at this location and in these years, forcing uncertainty may still prevail. Figure 5  maximum and minimum daily temperatures observed at Stations #508 and #463 with NLDAS daily temperature from 1 April to 20 May 2004. During this period, NLDAS temperature data were comparatively cooler than observation, and the difference in maximum air temperature can affect the snow-melting process and time (Hamlet et al., 2005). Nonetheless, the modified approach has improved SWE estimation. Figure 6 shows the components of water balanceprecipitation, sublimation, and snowmelt for the winter of 2001-2002 at California SNOTEL stations #356, #508, #463, and #539. During the accumulation period (December through March), loss of snow in the control model is because of both sublimation and snow melting processes but in the modified model, the loss is mainly from sublimation. As discussed earlier, the control model simulates snowmelt whenever the temperature of the snowpack reaches freezing point, and then melt water immediately becomes runoff. On the contrary, in the modified approach, snowmelt commences only when the net energy relative to the melting point is positive and snow melting do not start until later in the spring season (Fig. 6d, h, l, and p). During the accumulation period, loss due to sublimation is less in the modified model compared to the control model (Fig. 6c, g, k, and o) which can be attributed to the application of different formulation in the models.
In the modified model, liquid-water holding capacity of the snowpack was considered before the melt water becomes runoff and the effect of the liquid water content is shown in Fig. 7. The modified model was simulated with 0 and 5 % liquid water content, and the difference in response is seen only during the ablation period. There is no significant change in melt outflow rate until the beginning of the melt period. The simulation run with 5 % liquid-water holding capacity delays the onset of snowmelt, compared to that of 0 % liquid-water holding capacity for less than a day to approximately a few days at some study sites.
An additional review of the effectiveness of the modified model in predicting maximum SWE is presented in Fig. 8, where maximum-modeled SWE as a percentage of ground-observed maximum SWE is shown at the four California SNOTEL stations for 7 years of the study period. Although precipitation bias correction improved model SWE estimation, the overall enhancement in maximum SWE prediction by the modified model is evident (Fig. 8). A similar comparison is shown in Fig. 9, which includes the 21 California SNOTEL stations over the Sierra Nevada Mountains. The control model can reasonably predict SWE at locations where maximum SWE is relatively less (SWE max < 500 mm). However, bias is more pronounced at locations with higher snowfall which are typically located at higher elevation where uncertainty of input variables is greater. The modified model has enhanced SWE estimation at all locations, but improvement is more prominent for observation stations where the maximum snowfall is between 500-1000 mm.
Additionally, the modified model's overall predictive power to simulate SWE is described by RMSE, RSR in Table 2, and the Nash-Sutcliffe efficiency in Fig. 10. For all stations, RMSE from control model simulation is much higher than the threshold RMSE value (i.e., half the standard deviation of observed data). Although, the modified model simulation lowers the values of RMSE than that of the control model, majority of the station RMSE is still considered high. RSR performance ratings describe most of the higher RMSE simulations as "satisfactory" or "good". Figure 10 shows the NSE values of the modified model simulation which is positive for the four stations with most of the values above 0.5. However, only few years show NSE values were >= 0.7 which supports, like RMSE and RSR values, the modified model performance is satisfactory in SWE estimation.

Simulation -other variables
At the California SNOTEL sites, other properties of snowpack, such as snow surface temperature, the surface energy flux measurements are not available. So a comparison between the control model and the modified model simulated snow surface temperature and turbulent fluxes at SNOTEL Station #356 are shown in Fig. 11. The figure compares the model simulation from 1 April to 1 May 2002. This time period was chosen to present the difference between model simulations at the beginning of the melt period (1 April is considered as the beginning of snow ablation period). The figure shows that both the models have similar temperature distribution, although the control model is simulating a colder surface compared to the modified model (Fig. 11a). However, a significant difference between the models can be seen in  the snowmelt outflow rate (Fig. 11b). During this time period, frequent melt events in the control model reduced the snowpack from the ground while modified model did not simulate any snowmelt. The differences in sensible heat flux (Fig. 11c) and outgoing long-wave radiation (Fig. 11d) between the control model and the modified model are found to be small. Latent heat flux computed by the modified model is larger compared to that of the control model (Fig. 11e) because of the modified model simulating a warmer snow surface.
Snow-covered ground can affects soil temperature, as well as moisture content of the underlying soil. Although the scope of the paper is only limited to evaluating alternative processes for snow temperature and snowmelt, improvement in soil temperature and moisture in the modified model is Hydrol. Earth Syst. Sci., 18, 3553-3570, 2014 www.hydrol-earth-syst-sci.net/18/3553/2014/ Table 2. Model error statistics at four SNOTEL stations. Model performance ratings are given in the parenthesis where h, l, vg, g, s, and us refers to "high", "low", "very good", "good", "satisfactory", and "unsatisfactory" performance, respectively.
Station 2001-2002 2002-2003 2003-2004 2004-2005 2005-2006  compared with the respective observations at the SNOTEL stations. Figure 12 shows the comparison of model outputs with observed soil temperature and soil moisture at 5 mm below the ground surface at Stations #508 and #463. The control model predicts less snow and therefore, the ground is more exposed to the cooler atmosphere and soil is below freezing point throughout the snow season ( Fig. 12c and  d). Modified model simulation shows comparatively warmer soil temperature and observation shows soil temperature was above the freezing point during most of the snow season. Similarly, because of simulating less snow on the ground and more frequent snowmelt, control model predicts a higher soil moisture fraction compared to observation. On the contrary, with less snowmelt, the modified model's soil moisture simulation agrees better with observation ( Fig. 12e and f). Note that the observation sites initialize the measuring instrument at every water year and so for the first few months of the water year, the soil moisture is recorded as zero. Additional analysis of improving the soil temperature and moisture content is suggested but is beyond the scope of this study.

Utah site
The modified model was also assessed at Utah SNOTEL Station #1098 near the TWDEF forest and simulation result is shown in Fig. 13. At this site, the SWE predicted by the control model shows negative bias similar to that of California sites. The modified model shows improvement over the control model by reducing the SWE bias and delaying ablation period. Figure 14 snows the comparison of simulated and recorded snow surface temperature for early March (accumulation period) and May (ablation period), 2009. The snow surface temperature simulated by both the control and Noah-T s reasonably agree with observed snow surface temperature although the later predicted slightly warmer temperature than that of the control model. The comparison supports that the control model does not benefit from the snow surface temperature formulations and enhancement in SWE simulation is mainly from the controlling snow melting based on internal energy of the snowpack.

Discussion
Noah is a 1-D model which can be simulated using grid or point observational data. As Noah LSM is coupled with regional-climate models for weather simulation, the model was simulated using NLDAS 1/8 • grid data for this study. Simulation results are compared with point SNOTEL recorded SWE because of lack of reliable large scale observations of snow data . Therefore, to account for the scale issues, precipitation forcing data has been bias corrected using California SNOTEL's precipitation which improves control model's performance. The applied model modifications enhance SWE simulation by increasing amount of snowfall and delaying ablation period. Similar improvement in Noah's SWE modeling is seen at a Utah SNO-TEL site where the modified model was simulated using observed precipitation and temperature forcing data. Overall, the modified model removes the control model's deficiency of early melting but not adequately efficient in removing SWE bias. Modeling physical processes of snow accumulation/sublimation is much more complex in the mountainous terrain and therefore consideration of other factors, such as time-varying snow albedo parameterization and effect of canopy shading are important. Snow albedo parameterization along with other modifications suggested by various studies has been incorporated in the recent versions (current version is 3.4.1). Therefore, applying the modifications of Noah-Ts in the current Noah version is encouraging to study their combined effect on snow process modeling.

Conclusions
The Noah LSM has been identified to under predict snowwater equivalent throughout the snow season and melt away all the snow early in spring. The control model simulates snowpack as a single layer and estimates the snow surface temperature based on the energy flux on the snow surface and the air temperature above the snow surface. When the snowpack is at freezing temperature (0 • C), snow starts to melt. In reality, snow temperature is not the only determining factor for snowmelt. To begin snow melting, the entire snowpack is required to be isothermal and the modified Noah model, like UEB model, determines this by accounting for internal energy of the snowpack. This melt water does not instantly become runoff but remains within the warm snowpack up to the liquid-water holding capacity. Once this water-holding capacity is exceeded, snowmelt becomes runoff.
In addition to considering internal energy of the snowpack, the force-restore method was applied in modified Noah model to compute snow-surface temperature. The results show that snow surface temperature is similar to the control model and sometimes warmer than the temperature computed by the control model and thus the model did not benefit from the applied temperature scheme. The primary factor for improvement in modified model's SWE estimation is the regulating early season snow melting.
The new scheme adds only two prognostic variables and is compatible with other physical processes within the model. One effort of this study was to preserve the simplicity of the Noah model but remove model deficiencies. In general, the alternate single layer snow sub-model in the modified Noah LSM outperformed that of the control model by delaying snow melting and moderately removing SWE bias.
where T s is the surface temperature (K), T o is the air temperature (K), k is the thermal conductivity of soil (J m −1 s −1 K −1 ), ρ is the density of the air (kg m −3 ), C H is the turbulent exchange coefficient (m s −1 ), z is the depth of the snowpack plus half of the top soil layer (m), T soil is the soil temperature at half of the top soil layer (K), R dn is the net short-wave radiation (J m −2 s −1 ), ε is the emissivity of the surface (unitless), σ is the Stephen Boltzman's constant (J m −2 s −1 K −4 ), L f is the latent heat of fusion (J kg −1 ), C p is the specific heat of air (J kg −1 • K −1 ), P r is the precipitation rate (Kg m −2 s −1 ), θ o is the potential air temperature above the ground (K) (usually at 2 m above the ground), L v is the latent heat of condensation (J kg −1 ), β is the ratio of actual to potential evaporation (unitless), and E p is the potential evaporation rate (Kg m −2 s −1 ). Applying Taylor's expansion on εσ T 4 s : Dividing both sides by ρC p C H :