SnowClim v1.0: high-resolution snow model and data for the western United States

Seasonal snowpack dynamics shape the biophysical and societal characteristics of many global regions. However, snowpack accumulation and duration have generally declined in recent decades, largely due to anthropogenic climate change. Mechanistic understanding of snowpack spatiotemporal heterogeneity and climate change impacts will benefit from snow data products that are based on physical principles, simulated at high spatial resolution, and cover large geographic domains. Most existing datasets do not meet these requirements, hindering our ability to understand both contemporary and changing snow regimes and to develop adaptation strategies in regions where snowpack patterns and processes are important components of Earth systems.

We developed a computationally efficient process-based snow model, SnowClim, that can be run in the cloud. The model was evaluated and calibrated at Snowpack Telemetry (SNOTEL) sites across the western United States (US), achieving a site-median root-mean-squared error for daily snow water equivalent (SWE) of 64 mm, bias in peak SWE of −2.6 mm, and bias in snow duration of −4.5 d when run hourly. Positive biases were found at sites with mean winter temperature above freezing where the estimation of precipitation phase is prone to errors. The model was applied to the western US (a domain covering 3.1 million square kilometers) using newly developed forcing data created by statistically downscaling pre-industrial, historical, and pseudo-global warming climate data from the Weather Research and Forecasting (WRF) model. The resulting product is the SnowClim dataset, a suite of summary climate and snow metrics, including monthly SWE and snow depth, as well as annual maximum SWE and snow cover duration, for the western US at 210 m spatial resolution (Lute et al., 2021). The physical basis, large extent, and high spatial resolution of this dataset enable novel analyses of changing hydroclimate and its implications for natural and human systems.

1 Introduction

Seasonal snowpack shapes the climatic, hydrologic, ecological, economic, and cultural characteristics of many global regions. Snow is an important determinant of the surface energy balance through its effect on land surface albedo, partitioning of sensible and latent heat fluxes, near-surface atmospheric stability, and horizontal energy transport (Cohen, 1994; Rudisill et al., 2021; Stiegler et al., 2016). Hydrologic benefits of snow include natural water storage, delayed runoff, and cooler stream temperatures (Bales et al., 2006; Luce et al., 2014b). Ecologically, seasonal snow insulates flora and snow-dependent fauna, controls mobility and foraging opportunities, mediates nutrient cycling, and supplements plant water availability (Formozov, 1964; Grippa et al., 2005; Jones, 1999). Economically, seasonal snow helps to synchronize water and energy supply and demand, enables crop irrigation, fuels a multi-billion-dollar winter recreation industry in the United States (US) alone, and causes transportation delays and accidents (Burakowski and Magnusson, 2012; Qin et al., 2020; Seeherman and Liu, 2015; Sturm et al., 2017). Finally, seasonal snow is a defining aspect of many cultures globally, shaping language, traditions, and sense of self (Eira et al., 2013; Mergen, 1997).

In many mountain regions, recent decades have seen less precipitation falling as snow, lower peak snow water equivalent (SWE), shorter snow duration, and earlier snowmelt runoff (Choi et al., 2010; Fritze et al., 2011; Knowles et al., 2006; Mote et al., 2018). These developments are projected to continue in the coming decades, resulting in substantial declines (>50 %) in seasonal snowpack for areas such as the western US and significant impacts to human and natural systems (Fyfe et al., 2017; Huss et al., 2017; Marshall et al., 2019a; Siirila-Woodburn et al., 2021). In addition to these macroscale developments, there are important nuances to changing snow. Increased atmospheric water vapor due to warming is expected to enable larger snowfall events (Lute et al., 2015), which may buffer declines in snowpack (Kumar et al., 2012; Marshall et al., 2020). Changes in atmospheric circulation may affect snow accumulation, for example by diminishing orographic precipitation enhancement (Luce et al., 2013) or altering characteristics of atmospheric rivers (Dettinger, 2011). Decreasing snow cover will result in increased hydrologic importance of microclimates that serve as snow refugia, such as high elevations, deposition zones, and shaded areas (Marshall et al., 2019b; McLaughlin et al., 2017). A warmer and moister atmosphere will shift the relative importance of snowpack energy and mass budget terms, resulting, for example, in earlier but slower snowmelt (Musselman et al., 2017), changes to the partitioning of snow ablation between runoff and sublimation (Sexstone et al., 2018), and increasing rain-on-snow risk in regions that retain snow cover (Musselman et al., 2018).

Understanding these changes and their implications often requires snow models and modeled snow data products (hereafter snow data) that satisfy at least one of several criteria. These criteria include that the data (a) are simulated with physics-based representations of energy and mass transfer processes, (b) are spatially continuous, (c) have a high spatial resolution, (d) have a large extent, (e) are multivariate, and (f) are multitemporal. To address some questions about contemporary or future snow, the snow models themselves are needed and must be able to synthesize data that satisfies these criteria. Snow data developed from process-based equations for radiative and turbulent energy exchanges as opposed to temperature index approaches are argued to be necessary for both capturing the spatial variability of energy fluxes across the landscape and providing physically realistic simulations of the effects of climate change (Kumar et al., 2013; Raleigh and Clark, 2014; although see Lute and Luce, 2017). While it is increasingly clear that machine learning and artificial intelligence can emulate the net effect of physical processes to, for example, predict streamflow based on meteorological data (Fleming and Gupta, 2020), the ability of these approaches to predict snow under a changing climate has not been thoroughly evaluated. To assess changes in snowpack across a landscape, spatially continuous data are needed. In areas of complex terrain, high spatial resolution (<1 km²) data are necessary to resolve the effects of elevation and shading (Barsugli et al., 2020; Sohrabi et al., 2019; Winstral et al., 2014), which contribute to snow refugia that are important for species such as wolverine (Barsugli et al., 2020; Curtis et al., 2014). For some applications, such as water management and species distribution modeling, snow data may need to cover large geographic domains. Multiple snow metrics are needed for diverse applications. For example, SWE is commonly used for water management applications, whereas snow depth and density may be most relevant for wildlife applications at both large (e.g., for ungulate movement) and small scales (e.g., for wolverine denning sites). Finally, historical and future data are necessary to evaluate changes over time and to inform long-term planning and development of adaptation strategies for specific locales.

There are two major hurdles to the development of a snow dataset that meets all of these criteria: appropriate forcing data and computational cost. Presently, large-extent climate datasets only achieve horizontal resolutions of up to 1 km (e.g., Abatzoglou and Brown, 2012; Fick and Hijmans, 2017; Thornton et al., 2020) and the finer-resolution datasets cover limited domains or are restricted to historical periods (Dietrich et al., 2019; Holden et al., 2011, 2016). Second, even with appropriate forcing data, the computational expense of running snow models has generally forced the selection of some of these criteria at the expense of others (Winstral et al., 2014). For example, a temperature index model might be used for applications requiring rapid results over large domains (e.g., SNOW-17; Anderson, 2006), a process-based model might be run at high resolution over watershed-sized domains (Garen and Marks, 2005; Liston and Elder, 2006), or a process-based model might be run at coarser resolution over a larger extent (e.g., SNODAS, National Operational Hydrologic Remote Sensing Center, 2004; WRF, Rasmussen and Liu, 2017; Gergel et al., 2017; Wrzesien et al., 2018). There is potential for clever computational solutions and model formulations, such as variable resolution grids, to alleviate these trade-offs to some extent (Marsh et al., 2020).

We suggest that using a blended approach comprising process-based representations of the most crucial energy and mass balance fluxes (radiation and turbulent fluxes) and empirical simplifications for more computationally expensive (e.g., snow surface temperature) and/or typically minor (e.g., ground heat flux) components, implemented at high spatial resolution, can reduce computational cost relative to more process-based approaches and enhance accuracy of spatiotemporal snow simulations relative to coarser-scale implementations. The spatial resolution of the model implementation is a primary controlling factor on model accuracy in complex terrain, since the topographic smoothing inherent in coarser implementations can result in poor estimates of net shortwave radiation and precipitation-phase partitioning (Sohrabi et al., 2019). While subgrid parameterizations can achieve similarly low errors to fully distributed approaches (Luce and Tarboton, 2004), they do not provide spatially explicit simulations, which are useful for applications such as species distribution modeling and identification of topography-related snow refugia. Higher spatial resolution is especially important in the context of assessing future habitat, since snowpack response to climate change is strongly dependent on elevation and aspect (Barsugli et al., 2020). In contrast to existing process-based models that are implemented over large domains (e.g., VIC, Hamman et al., 2018; WRF, Ikeda et al., 2021), the present model does not calculate runoff and currently does not consider vegetation effects. Models that account for vegetation heterogeneity and dynamics may offer advantages over the present approach for specific applications.

In this study we developed a computationally efficient largely process-based snow model called SnowClim that has a flexible model structure and can be run in the cloud (Lute et al., 2021). The model retains the most important components of process-based models, including the complete energy balance and internal snowpack energetics, while omitting more computationally expensive components such as horizontal transport, multiple layers, and iterative solutions for snow surface temperature. Unlike existing models, this simplified process-based model is efficient enough to be run over sub-continental domains at high spatial resolution. We force the SnowClim model with pre-industrial (1850–1879), historical (2000–2013), and projected future (2071–2100) meteorological data from the Weather Research and Forecasting (WRF) model downscaled to correct for terrain effects. We then applied the model to the western US (a domain covering 3.1 million square kilometers) to create the SnowClim dataset, a multivariate, gridded, snow and climate dataset for three time periods at 210 m spatial resolution. Here we provide a description of the model and its application to the western US, including parameterization, calibration, climate forcing data preparation, and resultant datasets.

Figure 1 Snow model conceptual diagram. Solid black arrows indicate mass fluxes, and dashed grey arrows indicate energy fluxes. Fluxes are described in Sect. 2.

2 Model description

2.1 Model overview

The SnowClim model is a fully distributed energy and mass balance snow model. It simulates the snowpack as a single layer but accounts for different surface and pack temperatures (Fig. 1). The effects of vegetation, fractional snow cover, and snow redistribution via gravitational and wind-driven processes are currently not represented.

The model has a flexible structure to facilitate uncertainty analysis and application to new conditions. This flexible structure includes tunable parameters, a customizable spatiotemporal application, and process modularity. Key parameters (Table 2) are user defined (as opposed to hard coded in the model), allowing for calibration of the model to new conditions and regions as desired. The temporal and spatial resolution and extent are also user defined, which allows users to adjust to computational constraints and the requirements of the project. Finally, key processes such as albedo and turbulent fluxes are modularized to allow evaluation of alternative process representations.

The required forcings are described in Table 1. The model can be run in MATLAB (2020b) and requires the Statistics and Machine Learning Toolbox. The code is available at https://www.hydroshare.org/resource/dc3a40e067bf416d82d87c664d2edcc7/ (last access: 14 June 2022). The model can be run in the cloud using MATLAB Online through the HydroShare Platform hosted by the Consortium of Universities for the Advancement of Hydrologic Science, Inc. (CUAHSI).

Table 1 Required forcing data for the snow model.

Table 2 Parameters, their abbreviated names, the parameter values used in calibration, and their units. Additional parameter options, including the VIC model albedo option, were evaluated in preliminary work but were excluded from the full calibration due to consistently poor performance.

^a Essery et al. (2013). ^b Tarboton and Luce, (1996). ^c Parameter values chosen for the full model run by calibration at SNOTEL sites.

2.2 Energy balance

The SnowClim model evaluates the snow cover energy balance at each time step such that

$$\begin{array}{}\text{(1)}& Q_{\mathrm{net}}={\mathrm{SW}}_{\downarrow }-{\mathrm{SW}}_{\uparrow }+{\mathrm{LW}}_{\downarrow }-{\mathrm{LW}}_{\uparrow }+H+E_{\mathrm{i}}+E_{\mathrm{w}}+P+G,\end{array}$$

where Q_net is the net snow cover energy flux, SW_↓ is the downward shortwave radiation at the surface, SW_↑ is the upward shortwave radiation at the surface, LW_↓ is the downward longwave radiation at the surface, LW_↑ is the upward longwave radiation at the surface, H is the sensible heat flux, E_i and E_w are the latent heat fluxes of ice and water, P is the advected heat flux from precipitation, and G is the ground heat flux (Fig. 1).

2.2.1 Shortwave radiation

Upward shortwave radiation is equivalent to

$$\begin{array}{}\text{(2)}& {\mathrm{SW}}_{\uparrow }={\mathrm{SW}}_{\downarrow }\mathit{\alpha },\end{array}$$

where α is the spectrally integrated snow surface albedo.

Springtime snow model simulations are sensitive to the specific albedo algorithm (Etchevers et al., 2004; Günther et al., 2019). The SnowClim model provides three options for computing snow albedo (albedo_opt). In all options, albedo decays with time, and the albedo of shallow snowpacks (<100 mm depth) is diminished to account for the albedo of the ground surface, assumed to be 0.25 (Walter et al., 2005). A user-specified maximum albedo parameter (albedo_max) is used in each method.

The simplest albedo model (Essery et al., 2013; hereafter Essery) is empirical and sets albedo decay as a function of snowpack temperature. Snow albedo is augmented based on the occurrence and amount of new snow. Parameters other than the maximum albedo (minimum albedo, new snow threshold, linear and exponential albedo decline rates) are taken from Douville et al. (1995).

In the second albedo model (Hamman et al., 2018; Liang et al., 1994; hereafter VIC), snowpacks with new snow depth >10 mm and non-zero cold content receive the maximum snow albedo. Other albedo parameters are taken directly from VIC. Snow albedo decays more rapidly for melting snowpacks than cold snowpacks (cold content, cc<0).

The final albedo model (Dickinson et al., 1993; Tarboton and Luce, 1996; hereafter Tarboton) accounts for the wavelength dependence of albedo by computing separate visible and near-infrared band albedos as a function of snow surface age and solar illumination angle using a parameterization of global radiation (i.e., separate visible and near-infrared radiation are not supplied). The maximum albedo parameter is set equal to the average of the maximum visible band and infrared band albedos. This is the only albedo model of the three that includes a correction for illumination angle.

2.2.2 Longwave radiation

Upward longwave radiation is a function of snow surface temperature (T_s) in degrees Celsius, snow emissivity (ε), and the Stefan–Boltzmann constant (σ) such that

$$\begin{array}{}\text{(3)}& {\mathrm{LW}}_{\uparrow }=\mathit{\epsilon }\mathit{\sigma }(T_{\mathrm{s}}+\mathrm{273.15}{)}^{\mathrm{4}}+(\mathrm{1}-\mathit{\epsilon }){\mathrm{LW}}_{\downarrow }.\end{array}$$

We assume ε=0.98 (Armstrong and Brun, 2008). We consider T_s to be a function of the dew point temperature (T_d; Raleigh et al., 2013) such that

$$\begin{array}{}\text{(4)}& T_{\mathrm{s}}=\min (\mathrm{0}{}^{\circ }\mathrm{C},T_{\mathrm{d}}+T_{\mathrm{add}}),\end{array}$$

where T_add is an augmentation parameter that increases T_s and improves simulations of sublimation. For further discussion of T_s, see Sect. 2.2.6.

2.2.3 Turbulent fluxes

The turbulent fluxes, H, E_i, and E_w, are estimated using a Richardson number parameterization of the exchange coefficient following Essery et al. (2013). The bulk formula are

$$\begin{array}{}\text{(5)}& {\displaystyle }H=-{\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{a}}{C}_H{U}_{\mathrm{a}}(T_{\mathrm{s}}-T_{\mathrm{a}}),\text{(6)}& {\displaystyle }{E}_{\mathrm{i}}=-{\mathit{\rho }}_{\mathrm{a}}{C}_H{U}_{\mathrm{a}}(Q_{\mathrm{s}}-Q_{\mathrm{a}}){\mathit{\lambda }}_{\mathrm{s}}\text{ for }{T}_{\mathrm{s}}<\mathrm{0},\text{(7)}& {\displaystyle }{E}_{\mathrm{w}}=-{\mathit{\rho }}_{\mathrm{a}}{C}_H{U}_{\mathrm{a}}(Q_{\mathrm{s}}-Q_{\mathrm{a}}){\mathit{\lambda }}_{\mathrm{v}}\text{ for }{T}_{\mathrm{s}}=\mathrm{0},\end{array}$$

where ρ_a is the air density, c_a is the specific heat capacity of air, C_H is the bulk exchange coefficient that accounts for near-surface atmospheric stability, U_a is the wind speed, Q_s is the specific humidity of the snow surface, and Q_a is the specific humidity of the air (which is a required forcing). The specific humidity of the snow surface is calculated from T_s. The exchange coefficient C_H is parameterized as a function of the near-surface atmospheric stability as captured by the bulk Richardson number (Ri_B) such that

$$\begin{array}{}\text{(8)}& {\displaystyle }{C}_H=F_H\left(R{i}_{\mathrm{B}}\right)C_{\mathrm{HN}},\text{(9)}& {\displaystyle }R{i}_{\mathrm{B}}=\left(g{z}_u\left(T_{\mathrm{a}}-T_{\mathrm{s}}\right)\right)/\left(T_{\mathrm{a}}{U}_{\mathrm{a}}^{\mathrm{2}}\right),\text{(10)}& {\displaystyle }{C}_{\mathrm{HN}}=k^{\mathrm{2}}\left[\ln \left(z_u/z_{\mathrm{0}}\right){]}^{-\mathrm{1}}\right[\ln \left(z_T/z_{\mathrm{h}}\right){]}^{-\mathrm{1}},\text{(11)}& {\displaystyle }{F}_H\left(R{i}_{\mathrm{B}}\right)=\mathrm{1}\text{ for }R{i}_{\mathrm{B}}=\mathrm{0},\text{(12)}& {\displaystyle }\begin{array}{rl}{F}_H& \left(R{i}_{\mathrm{B}}\right)=\mathrm{1}-\left(\mathrm{3}cR{i}_{\mathrm{B}}\right)/\\ & \left(\mathrm{1}+\mathrm{3}{c}^{\mathrm{2}}{C}_{\mathrm{HN}}{\left(-R{i}_{\mathrm{B}}{z}_u/z_{\mathrm{0}}\right)}^{\mathrm{1}/\mathrm{2}}\right)\text{ for }R{i}_{\mathrm{B}}<\mathrm{0},\end{array}\text{(13)}& {\displaystyle }\begin{array}{rl}{F}_H\left(R{i}_{\mathrm{B}}\right)& =[\mathrm{1}+\left(\mathrm{2}cR{i}_{\mathrm{B}}\right)/{\left(\mathrm{1}+R{i}_{\mathrm{B}}\right)}^{\mathrm{1}/\mathrm{2}}{]}^{-\mathrm{1}}\\ \text{ for }R{i}_{\mathrm{B}}& >\mathrm{0},\end{array}\end{array}$$

where g is gravitational acceleration, z_u is the height of simulated wind speeds, z_T is the height of simulated air temperatures, z₀ is the surface roughness length for momentum, z_h is the surface roughness length for heat and water vapor, and c is a constant assumed to equal 5 (Louis, 1979). Both z₀ and z_h are adjustable user-specified parameters (Table 2).

An optional windless exchange coefficient is available to counter large radiative losses, particularly during stable conditions (Helgason and Pomeroy, 2012; Jordan, 1991). Application of the windless exchange coefficient can be modified through three parameters: E0_value, E0_app, and E0_stability (Table 2). E0_value is the value of the windless exchange coefficient (in W m⁻²). E0_app controls the application of the windless heat exchange coefficient to the sensible and latent heat fluxes: an E0_app value of 1 applies the coefficient only to the sensible heat flux, whereas an E0_app value of 2 applies the coefficient to both the sensible and latent heat fluxes. E0_stability controls the type of conditions where the windless coefficient is applied: an E0_stability value of 1 applies the coefficient to all conditions, whereas an E0_stability value of 2 applies the condition only under stable atmospheric conditions.

2.2.4 Precipitation heat flux

The heat flux of liquid precipitation is

$$\begin{array}{}\text{(14)}& P=c_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{w}}{T}_{\mathrm{d}}{P}_{\mathrm{rain}},\end{array}$$

where c_w is the specific heat of water, ρ_w is the density of water, and P_rain is the rate of liquid precipitation. The heat flux of solid precipitation (S) is handled separately for diagnostic purposes and is added directly to the snowpack cold content:

$$\begin{array}{}\text{(15)}& S=c_{\mathrm{i}}{\mathit{\rho }}_{\mathrm{w}}{T}_{\mathrm{d}}{P}_{\mathrm{snow}},\end{array}$$

where c_i is the heat capacity of ice and P_snow is the rate of snowfall.

2.2.5 Ground heat flux

The ground heat flux can be important in controlling the onset of seasonal snow accumulation, particularly in warmer environments (e.g., Mazurkiewicz et al., 2008). Many process-based snow models include or couple with a soil temperature model to simulate this flux. However, under most circumstances G is thought to provide a minor contribution to the energy budget (DeWalle and Rango, 2008). In the interest of model efficiency and to avoid uncertainties associated with estimating soil temperatures and thermal conductivities, we use a constant G of 2 W m⁻² (Walter et al., 2005), which is similar to other models (Etchevers et al., 2002).

2.2.6 Enhanced single-layer approach

Single-layer snow models typically provide less physically realistic snowpack simulations than multilayer models due to their simplified treatment of energy transfer within the snowpack (Blöschl and Kirnbauer, 1991; Waliser et al., 2011). Bulk single-layer conceptualizations treat the surface temperature and energy balance as synonymous with the pack temperature and energy balance, ignoring the contrast between the surface layer, which is highly sensitive to the near-surface atmosphere, and the deeper pack, which is characterized by thermal inertia, i.e., cold content. These distinctions are key to accurate modeling of snowpack heat fluxes (Blöschl and Kirnbauer, 1991) and snowpack ablation (Waliser et al., 2011).

To address these shortcomings, advanced single-layer snow models have differentiated between surface and pack temperatures, while attempting to maintain the parsimony of a single-layer model (Tarboton and Luce, 1996; You et al., 2014). To solve for snow surface temperature these approaches typically use iterative methods that can be computationally expensive (Wigmosta et al., 1994) or linearization approaches (Best et al., 2011).

The present model uses a two-step modification of the net snow cover energy flux to approximate the conduction of energy between the surface and the snowpack. This approach enables separate temperatures and energy balances for surface and pack components while retaining the computational efficiency necessary to accomplish the modeling objectives of both large spatial extent and relatively fine resolution. In this approach, the surface is conceptualized as a skin with zero depth.

First, we apply a temporal running mean to the net snow cover energy flux to approximate the attenuation with depth of the characteristic diurnal variations in energy at the surface, akin to the approach taken by You et al. (2014). The smoothed energy flux from the surface to the pack at each time step ($\stackrel{\mathrm{‾}}{Q_{\mathrm{net}}}$) is calculated as the average net snow cover energy flux over a period smooth_hrs, which is a tunable parameter (Table 2). This approach reduces unrealistic high-frequency modifications of the cold content and large-amplitude freeze–thaw cycles during the ablation season.

Second, we apply a progressive tax on the negative net energy flux to the snowpack to limit the excessive accumulation of cold content that results from all snow cover energy being directly translated to the pack. The net effect of the energy tax is to reduce snowpack cold content, resulting in more accurate cold content simulations relative to observations (Table A1) and similar to those from other, more complex process-based models (Jennings et al., 2018a). Other single-layer models have sought to limit cold content; however, they used approaches that required site-specific calibration (Blöschl and Kirnbauer, 1991; Braun, 1984). We apply a progressive tax such that negative energy fluxes to snowpacks with larger cold content receive larger taxes.

$$\begin{array}{}\text{(16)}& \begin{array}{rlr}& Q_{\mathrm{pack}}=\stackrel{\mathrm{‾}}{Q_{\mathrm{net}}}& \text{ for }\stackrel{\mathrm{‾}}{Q_{\mathrm{net}}}\ge \mathrm{0}\\ & Q_{\mathrm{pack}}=\stackrel{\mathrm{‾}}{Q_{\mathrm{net}}}\times (\mathrm{1}-\mathrm{tax})& \text{ for }\stackrel{\mathrm{‾}}{Q_{\mathrm{net}}}<\mathrm{0}\\ & \mathrm{tax}={\displaystyle \frac{\mathrm{cc}-{\mathrm{cc}}_{\mathrm{0}}}{{\mathrm{cc}}_{\mathrm{1}}}}\times \mathrm{maxtax}\text{ such that }\mathrm{0}\le & \mathrm{tax}\le \mathrm{maxtax}\end{array}\end{array}$$

In Eq. (16), $\stackrel{\mathrm{‾}}{Q_{\mathrm{net}}}$ is the smoothed net snow cover energy flux, Q_pack is the energy flux from the surface to the pack, and cc is the snowpack cold content, which uses a negative sign convention. cc₀, cc₁, and maxtax are tunable parameters that define the maximum (least negative) cold content to which the tax should be applied, the range of cold content over which the tax should be applied (cc₀ to cc₀+cc₁), and the maximum possible tax, respectively (Table 2). Negative energy fluxes to snowpacks with cold contents less negative than cc₀ receive zero tax, and negative energy fluxes to snowpacks with cold contents more negative than cc₀+cc₁ receive a tax equal to maxtax.

Q_pack is added to the snowpack cold content (cc) at each time step. Pack temperature (T_pack) can be obtained from cold content as follows:

$$\begin{array}{}\text{(17)}& T_{\mathrm{pack}}=\mathrm{cc}/\left({\mathit{\rho }}_{\mathrm{w}}\times c_{\mathrm{i}}\times \mathrm{SWE}\right),\end{array}$$

where SWE is the snow water equivalent.

2.2.7 Modification for shallow snowpacks

We developed a computationally efficient approach for controlling energy balance instabilities for shallow snowpacks. Marks et al. (1999) addressed the problem by shifting to progressively smaller time steps. In the interest of computational efficiency, we take an alternative approach. When modeled SWE is less than a threshold value, T_pack is set equal to the minimum of T_a and 0 ^∘C. Cold content is then updated according to this new temperature. The threshold for applying this correction is 15 mm of SWE for every hour in the time step (e.g., for a model run at a 4 h time step, the temperature correction would be applied to snowpacks with 60 mm SWE or less). Constraining T_pack and cold content in this way is reasonable given that surface and pack temperatures are likely to be similar for shallow snowpacks and the strong correspondence between T_s and T_a (Helgason and Pomeroy, 2012).

2.3 Mass balance

The mass balance of the solid and liquid portions of the snowpack are evaluated at each time step as follows:

$$\begin{array}{}\text{(18)}& {\displaystyle }{M}_{\mathrm{s}}=M_{\mathrm{snow}}+M_{\mathrm{ref}}-M_{\mathrm{melt}}+M_{\mathrm{dep}}-M_{\mathrm{sub}},\text{(19)}& {\displaystyle }{M}_{\mathrm{l}}=M_{\mathrm{rain}}-M_{\mathrm{ref}}+M_{\mathrm{melt}}-M_{\mathrm{runoff}}+M_{\mathrm{cond}}-M_{\mathrm{evap}},\end{array}$$

where M_s is the mass of the solid portion of the snowpack, M_snow is the mass of new snowfall, M_ref is the mass of liquid water in the snowpack that has been refrozen, M_melt is the mass of snow that has melted, M_dep is the mass of deposition, M_sub is the mass of sublimation, M_l is the mass of the liquid in the snowpack, M_rain is the mass of rain added to the snowpack, M_runoff is the mass of liquid water that has left the snowpack as runoff, M_cond is the mass of condensation, and M_evap is the mass of evaporation (Fig. 1).

2.3.1 Accumulation

Snowfall is calculated as an air-temperature- and relative-humidity-dependent fraction of precipitation using the bivariate logistic regression model of Jennings et al. (2018b). We use a non-binary formulation to allow for mixed-phase precipitation. New snowfall amounts less than 0.1 mm water equivalent per hour are set to 0. Rainfall is the difference between precipitation and snowfall. The temperature of new snowfall is set equal to the minimum of the dew point temperature and freezing point (0 ^∘C), whereas the temperature of rainfall is set equal to the maximum of the dew point temperature and the freezing point (Marks et al., 2013; Raleigh et al., 2013).

The density of new snowfall is calculated as a function of air temperature (Anderson, 1976) using constants identified by Oleson et al. (2004). Compaction of the snowpack is modeled as a function of SWE and snowpack temperature following Anderson (1976) and using constants from Boone (2002) for the ISBA-ES snow model. Snow depth is a function of SWE and density and is updated following changes in either variable.

2.3.2 Melt

Positive net energy flux must satisfy the snowpack cold content before melt can occur. Melt is equivalent to the minimum of the current SWE and the potential melt,

$$\begin{array}{}\text{(20)}& {\mathrm{melt}}_{\mathrm{pot}}=Q_{\mathrm{pack}}/({\mathit{\lambda }}_{\mathrm{f}}\times {\mathit{\rho }}_{\mathrm{w}})\text{ for }{Q}_{\mathrm{pack}}>\mathrm{0},\end{array}$$

where λ_f is the latent heat of freezing.

2.3.3 Liquid water content

Rainfall, melt, and condensation are added to, and evaporation is subtracted from, the snowpack liquid water content. Snowpack liquid water content in excess of the liquid water holding capacity of the snowpack contributes to runoff. The liquid water holding capacity of the snowpack is the product of snow depth and the maximum liquid water fraction (lw_max, Table 2). Liquid water content below this threshold but greater than the minimum liquid water content (equivalent to 1 % of snow depth; Marsh, 1991) is allowed to drain at a rate of 100 mm h⁻¹ (based on values in DeWalle and Rango, 2008).

2.3.4 Refreezing

Excess cold content can be used to refreeze liquid water in the snowpack. The amount of water refrozen is the minimum of the total liquid water content of the snowpack and the potential refreezing,

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathrm{refreeze}}_{\mathrm{pot}}=-\mathrm{cc}/\left({\mathit{\lambda }}_{\mathrm{f}}\times {\mathit{\rho }}_{\mathrm{w}}\right)\text{ for }\mathrm{cc}<\mathrm{0},\text{(22)}& {\displaystyle }{M}_{\mathrm{ref}}=\min \left({\mathrm{refreeze}}_{\mathrm{pot}},M_{\mathrm{l}}\right)\text{ for }\mathrm{cc}<\mathrm{0}.\end{array}$$

Energy released by refreezing is added to the snowpack cold content, and the refrozen mass is added to the SWE, increasing the snowpack density (we assume no change in snow depth).

2.3.5 Sublimation and condensation

Latent heat transfer results in sublimation or evaporation from (or deposition or condensation onto) the snowpack, such that

$$\begin{array}{}\text{(23)}& {\displaystyle }{M}_{\mathrm{sub}}=-E_{\mathrm{i}}/({\mathit{\lambda }}_{\mathrm{s}}\times {\mathit{\rho }}_{\mathrm{w}})\text{ for }{E}_{\mathrm{i}}<\mathrm{0}\text{ and }{T}_{\mathrm{s}}<\mathrm{0},\text{(24)}& {\displaystyle }{M}_{\mathrm{evap}}=-E_{\mathrm{w}}/({\mathit{\lambda }}_{\mathrm{v}}\times {\mathit{\rho }}_{\mathrm{w}})\text{ for }{E}_{\mathrm{w}}<\mathrm{0}\text{ and }{T}_{\mathrm{s}}=\mathrm{0},\text{(25)}& {\displaystyle }{M}_{\mathrm{dep}}=-E_{\mathrm{i}}/({\mathit{\lambda }}_{\mathrm{s}}\times {\mathit{\rho }}_{\mathrm{w}})\text{ for }{E}_{\mathrm{i}}>\mathrm{0}\text{ and }{T}_{\mathrm{s}}<\mathrm{0},\text{(26)}& {\displaystyle }{M}_{\mathrm{cond}}=-E_{\mathrm{w}}/({\mathit{\lambda }}_{\mathrm{v}}\times {\mathit{\rho }}_{\mathrm{w}})\text{ for }{E}_{\mathrm{w}}>\mathrm{0}\text{ and }{T}_{\mathrm{s}}=\mathrm{0},\end{array}$$

where λ_s is the latent heat of sublimation and λ_v is the latent heat of vaporization.

3 Model application to the western United States

The SnowClim model was evaluated and calibrated at a collection of automated snow stations across montane portions of the western US and further applied to the broader western US to create the SnowClim dataset. We describe the preparation and downscaling of the meteorological forcing data, the model calibration, and the model simulations for the western US. The model was calibrated at Snowpack Telemetry (SNOTEL) sites and model performance at these sites was used to select the parameters and temporal resolution at which to run the model over the full domain.

3.1 Spatial resolution

To balance the competing ambitions of high spatial resolution and computational feasibility over the western US domain, we used variable spatial resolutions. Regions of complex terrain were modeled at 210 m (hereafter “fine”). This high resolution enhances the model's ability to capture the effects of elevation, aspect, and slope on snowpack in complex terrain. Regions of less complex terrain were modeled at 1050 m (hereafter “coarse”). Terrain complexity was assessed for each coarse grid cell by examining the elevations and downscaled shortwave radiation values for the 25 colocated fine grid cells. If the elevation difference across the fine cells was less than 50 m and the maximum percent difference in shortwave radiation was less than 10 %, then snow simulations were completed at coarse resolution. Otherwise, simulations were completed at fine resolution. This resulted in approximately 30 % of the domain (920 605 grid cells) being modeled at coarse resolution (Fig. B1), with the remainder (64 310 454 grid cells) being modeled at fine resolution. Grid cells were defined using the 1 arcsec National Elevation Dataset Digital Elevation Model (DEM; Gesch et al., 2018) and aggregated to either 210 or 1050 m.

Figure 2 Performance metrics for an hourly model run with the selected parameterization.

Figure 3 Time series of observed and modeled SWE at the Summit Ranch, Colorado, SNOTEL site. Out of all 170 SNOTEL sites, errors at this site were closest to the all-station median errors reported in the text.

3.2 Forcing data preparation

Hourly meteorological data from the Weather Research and Forecasting model (WRF; Rasmussen and Liu, 2017) were downscaled to force the snow model (Table 3). Forcing data were developed for a historical period, future period, and pre-industrial period. The raw WRF data consisted of 4 km spatial resolution hourly simulations for 1 October 2000 to 30 September 2013 that used initial and boundary conditions from ERA-Interim (Dee et al., 2011), herein referred to as the historical period. A pseudo-global warming run was also performed by perturbing ERA-Interim by average differences from a suite of climate models participating in the Fifth Coupled Model Intercomparison Project (CMIP5; Taylor et al., 2012) between 1976–2005 and 2071–2100 under the high-warming RCP 8.5 scenario (Rasmussen and Liu, 2017). Pre-industrial forcing data were developed by perturbing the downscaled historical WRF data by monthly climatological differences in climate between pre-industrial (1850–1879) and the historical period using a pattern-scaling approach (Mitchell, 2003) based on spatially varying differences in variables from the CMIP5 models.

Table 3 WRF data used to derive forcing data for the snow model.

Spatial downscaling for all variables except shortwave radiation was accomplished using moving window lapse rates (i.e., the change in the variable with elevation). Lapse rate downscaling has been shown to perform well relative to other statistical downscaling approaches in mountainous terrain (Praskievicz, 2018; Wang et al., 2012). We estimated monthly lapse rates for each grid cell and each variable, except for temperature, for which we estimated hourly lapse rates for each grid cell. Windows of 7×7 WRF grid cells (or 28×28 km) were used to balance the competing objectives of sufficient data points and the ability to capture local phenomena (Lute and Abatzoglou, 2021). Lapse rate corrections were applied hourly using the elevation difference between the WRF grid cell and the target DEM grid cell. For air pressure, lapse rates were calculated from and applied to temporally averaged WRF data. Grid cells not classified as land by WRF were excluded from lapse rate calculations.

For precipitation, a modified version of the methods above was used. Prior to calculating lapse rates, WRF precipitation was bias corrected to monthly 4 km precipitation from PRISM (PRISM Climate Group, 2015) by calculating monthly correction ratios, i.e., the ratio of total monthly PRISM precipitation to total monthly WRF precipitation. Correction ratios were set to 1 (no correction) when monthly WRF precipitation was 0 or when the ratio was infinite. Monthly precipitation lapse rates were divided by the number of hours with precipitation each month in the underlying WRF data, and hours with zero precipitation were maintained in the downscaled data to avoid precipitation every hour due to non-zero monthly lapse rates.

Shortwave radiation was downscaled to the target DEM using the insol package in R (Corripio, 2015) following the approach of Lute and Abatzoglou (2021), which preserves the atmospheric effects (e.g., cloud cover) captured by WRF and also accounts for slope, aspect, self-shading, and shading by adjacent terrain. Parameters required by the algorithm, including visibility, RH, and temperature, were assumed to be constant. Terrain corrections were calculated for the midpoint of each hour of the middle day of each month, aggregated to the desired temporal resolution using a weighting scheme based on the amount of solar radiation each hour, and then interpolated to the full time period.

For model calibration at SNOTEL sites (see Sect. 3.3.2), the above downscaling procedures were applied, but values were adjusted based on the elevation difference between the SNOTEL site and the colocated WRF grid cell based on calculated lapse rates. Downscaled WRF precipitation was bias corrected to SNOTEL sites by applying a monthly correction factor consisting of the ratio of the total SNOTEL precipitation to the total WRF precipitation similar to Havens et al. (2019). We note that such bias correction approaches may not address issues of precipitation undercatch at SNOTEL sites.

Additional variables needed to force the snow model, including specific humidity, relative humidity, and dew point temperature, were derived from the downscaled WRF water vapor mixing ratio, air temperature, and air pressure data using standard methods. Dew point temperatures exceeding the air temperature were set equal to the air temperature.

Table 4 Spatial correlations (R²) between observations at SNOTEL sites and SnowClim simulations for various snow metrics over the model calibration period 1 October 2000–30 September 2013.

Table 5 Summary climate and snow variables included in the SnowClim dataset. Summary variables are available for pre-industrial, historical, and future time periods.

Figure 4 Snow model performance for different time steps using the parameter set selected in calibration of the hourly model. Points represent median values across 170 SNOTEL sites.

3.3 Model calibration

3.3.1 Calibration methods

The model was calibrated at SNOTEL sites across the mountains of the western US to select a single best parameter set across all sites. While the SNOTEL network may not be perfectly representative of the broader domain (Molotch and Bales, 2005), it represents the best available ground observation dataset for calibration purposes due to the number of sites, range of hydroclimate conditions monitored, length of record, and relatively consistent observational methods and equipment across sites. A total of 170 SNOTEL sites were selected meeting the following requirements:

1. an elevation difference of less than 75 m relative to the colocated WRF grid cell,

2. a dataset missing no more than 1 % of daily precipitation and SWE observations between October and May in every water year between 1 October 2000 and 30 September 2013,

3. a site located more than 25 km from any other SNOTEL site.

Missing SWE values were infilled using linear interpolation across time. Missing precipitation values were infilled using an inverse distance weighted average of the values at the three closest sites.

Calibration consisted of running the model across all SNOTEL sites for each possible combination of parameters listed in Table 2. Model performance was assessed using the mean absolute percent error (MAPE) of annual maximum SWE (maxswe), the MAPE of annual snow duration, and the root-mean-squared error (RMSE) of daily SWE at each site. Snow duration was defined as the duration (in days) of the longest period of consecutive days with SWE >0. RMSE was computed for days when observed SWE exceeded 10 mm. Additionally, we used the mean error (ME) and mean percent error (MPE) of maxswe and duration to visualize calibration errors. The optimal parameter set was selected using Pareto preference ordering (Khu and Madsen, 2005) based on the median of each statistic across stations.

The model was subsequently evaluated for different run time steps (1, 2, 3, 4, 6, 8, 12, and 24 h). Separate model calibration for each time step selected similar parameter sets to the hourly run, and thus the hourly parameter set was used for all time steps. Model performance was again assessed as described above.

3.3.2 Calibration results

Snow model calibration via Pareto optimization selected a single best parameter set (Table 2). The station median MAPE of maxswe, MAPE of snow duration, and daily RMSE for this parameter set were 15.9 %, 8.43 %, and 64.0 mm, respectively. The spatial distribution of ME and MPE in maxswe and duration lacked strong coherent spatial patterns (Fig. 2), and spatial correlations (R²) for a variety of snow metrics exceeded 0.75 (Table 4), suggesting that the model captured major climate-related effects and sources of large-scale spatial snow variability. The largest negative biases were found at drier sites with relatively shallow or intermittent snowpacks (Fig. B2). The largest positive biases were found at sites with mean winter temperatures at or above freezing, where snow accumulation is very sensitive to the partitioning of precipitation into rain vs. snow (Fig. B2). A time series of observed and modeled SWE at one site with error values close to the station median values illustrates the model performance on a daily scale (Fig. 3). The model also captured key components of interannual snowpack variability over the short historical period; the station median correlation coefficients for maxswe and for snow duration were 0.93 and 0.70, respectively. The station correlations did not demonstrate any clear geographic or climatic patterns. This lends confidence to the model's ability to accurately simulate snow dynamics across a range of climates. The parameter sensitivity of the model is shown in Fig. B3.

Model performance deteriorated as temporal resolution coarsened from 1 to 12 h but improved slightly for the 24 h time step (Fig. 4). Model performance was somewhat sensitive to the hours used for aggregation; other aggregation windows showed continuous performance deterioration with coarsening temporal resolution (not shown). The sensitivity of model performance to the aggregation window is likely related to how diurnal energy fluxes, particularly shortwave radiation, are aggregated. A time step of 4 h was selected for the full western US model run to balance the objectives of computational efficiency and model performance. The station median MAPE of maxswe, MAPE of snow duration, and RMSE for the 4 h time step were 17.6 %, 8.31 %, and 69.5 mm, respectively. We note that simulations without the modification for shallow snowpacks (Sect. 2.2.7) degraded more consistently and significantly with coarsening temporal resolution (Fig. B4).

Figure 5 (a) Historical and (b) future maxswe (mm), (c) historical and (d) future snow duration (days). Historical values are averages over the period 2000–2013. Future values represent averages during the period 2071–2100 under RCP 8.5. In (a) and (b), white land areas denote areas that had less than 50 mm maxswe. In (c) and (d), white land areas denote areas where snow duration was 0. Note the nonlinear color scale in panels (a) and (b).

Figure 6 (a) Absolute and (b) percent change in maxswe between historical and future periods. (c) Absolute and (d) percent change in snow duration between historical and future periods. In (d), small boxes in Utah and Colorado indicate the regions highlighted in Figs. 7 and 8, respectively.

Figure 7 Example of simulations of changing maxswe for a portion of the Uinta Mountains, Utah (location is marked in Fig. 6). The elevation (m) of the domain is shown in (a). The percent change (%) in maxswe between historical and late 21st century periods as simulated by a 4 km WRF product (Rasmussen and Liu, 2017) is shown in (b), and the same metric but from the SnowClim dataset is shown in (c).

Figure 8 Snow depth on 20 May 2010 in the Boulder Creek Watershed area indicated in Fig. 6 from (a) lidar (described in Harpold et al., 2014) and (b) SnowClim.

3.4 Model results for the western United States

The SnowClim model was applied to the western US (contiguous US west of 104^∘ W) using the parameters identified above, a temporal resolution of 4 h, and a variable spatial resolution as described previously (210–1050 m horizontal resolution). The model was run in parallel on a high-performance standalone server (Dell Poweredge R730) with 34 cores and 128 GB RAM. The compute time for downscaling the climate forcings and executing the snow model was 10 and 2.5 d, respectively, for the historical period. For reference, the model took less than 0.03 s per site year when run using a single core across the 170 SNOTEL sites.

Large-scale spatial patterns of climatologies and changes in maxswe and snow duration (Figs. 5 and 6) were broadly similar to those from existing products developed from a wide range of modeling approaches (Luce et al., 2014; Lute et al., 2015; Ikeda et al., 2021). Historical maxswe was 112 mm, spatially averaged across the full western US domain, and locations with historical maxswe <50 mm were found in the warmer southern and southwestern regions and in the northeastern portion of the domain where winters are relatively dry (Fig. 5a). Under the future scenario, spatially averaged maxswe declined to 54 mm, and the areas with maxswe <50 mm expanded to encompass many lower-elevation areas (Fig. 5b). Historical snow duration averaged 83 d, spatially averaged across the full western US domain (Fig. 5c), but declined to 43 d in the future scenario (Fig. 5d). There were only a handful of locations with increases in maxswe or duration in the future period compared with the historical period, and these increases were small (Fig. 6). The largest relative declines in maxswe and duration were found at low elevations. On average, maxswe decreased by 49 % across locations with at least 50 mm maxswe historically, and snow duration decreased by 61 % across locations with historical snow duration greater than zero. Summaries of historical and future maxswe and duration by four-digit hydrologic unit code (HUC) are provided in Table B1.

Compared to existing large extent, multitemporal, process-based snow datasets such as that from the 4 km WRF runs (Rasmussen and Liu, 2017), SnowClim provided a much more nuanced picture of changing snow, particularly in areas of complex terrain. For example, Fig. 7 shows relative changes in maxswe for the Uinta Mountains in northeastern Utah as simulated directly by the 4 km WRF product and by SnowClim. SnowClim captured effects of elevation and aspect, including greater percent reductions in maxswe at lower elevations and on south-facing aspects, similar to Barsugli et al. (2020). Nuanced results such as these are only possible with high-resolution snow modeling that explicitly simulates spatiotemporal variations in the dominant snow cover energy fluxes.

Comparison of SnowClim snow depth with a snapshot of finer-resolution lidar-based observations further highlights some of the strengths and limitations of SnowClim. For example, in the Boulder Creek Watershed, Colorado, SnowClim captured the broad-scale spatial patterns of snow depth that are present in lidar-derived depth observations (described in Harpold et al., 2014) aggregated from the original 1 m resolution to the resolution of SnowClim on 20 May 2010 (Fig. 8). However, SnowClim snow depth was less spatially variable, particularly in the higher-elevation western portion of the domain. Evaluations such as this are particularly challenging as they simultaneously evaluate the spatial and temporal fidelity of the forcing data, the snow model, and in particular the snow density algorithm. The muted spatial variability in SnowClim can be attributed to a combination of factors, chief of which may be the lack of wind–snow interaction in the current SnowClim model formulation. Snow redistribution by wind and blowing snow sublimation have a significant effect on snowpack heterogeneity in this region (Knowles et al., 2012; Sexstone et al., 2018; Winstral et al., 2002); one study showed that a model incorporating wind redistribution captured 8 %–23 % more of the spatial variability in snow depth than a model without these processes (Winstral et al., 2002). In order to better simulate snowpack in windy environments such as the alpine area shown here, incorporation of blowing snow transport and sublimation into future SnowClim model formulations should be considered.

4 Discussion and conclusion

Through the development of a new computationally efficient snow model, SnowClim, and novel forcing data, we have overcome the two major hurdles to achieving snow data that meets the criteria outlined in Sect. 1. SnowClim's distinctive balance of mostly process-based and some empirical components allows it to capture contrasts in radiative loading in complex terrain, timing and rate of ablation, and responses to future climate, while maintaining computational efficiency. The SnowClim dataset is spatially continuous across the western US at sub-kilometer resolution in complex terrain, enabling both high-resolution and large-extent analyses. In particular, the SnowClim model and dataset highlight the effects of elevation and aspect on snowpack in a changing climate. The inclusion of multiple snow variables and compatible climate variables across multiple time periods will empower analyses of hydroclimatic responses to changing climate and complement existing coarser-resolution products.

The SnowClim model excludes some processes that might be included in more complex, computationally expensive large-scale models (e.g., VIC, WRF) such as vegetation-related processes, blowing snow transport and sublimation, and gravitational redistribution. In some contexts, these processes may be necessary for accurate modeling of the snowpack (Freudiger et al., 2017; Musselman et al., 2008; Pomeroy et al., 1993). The western US SnowClim dataset was calibrated at SNOTEL sites, which are often in forest clearings, and it is therefore expected to be relatively accurate in similar environments. However, user judgment should be applied when using the model or dataset in vegetated areas. Given the complexity of vegetation–snow processes, incorporation of vegetation effects may add significant computational expense and is hindered by the need for vegetation-related data and parameters that are expected to change between the time periods considered here. Future vegetation change is subject to large uncertainty stemming from disturbance regimes and species composition pathways. However, incorporation of an optional vegetation routine to be used when data and computational resources are available is a logical next step. Approaches for incorporating blowing snow transport either require (i) high-resolution wind fields input to semi-empirical or 3D turbulent-diffusion models (summarized in Mott et al., 2018), requiring more sophisticated downscaling of wind fields than what was done here and substantially increasing computational cost, or (ii) require calibration of terrain parameters (e.g., Winstral et al., 2013), which would be possible but both challenging and computationally intensive for a large-scale model such as SnowClim. Simple algorithms do exist for modeling gravitational redistribution of snow (Bernhardt and Schulz, 2010); however, incorporation of either blowing snow transport or avalanching would necessitate restructuring of the model as a semi-distributed or fully distributed model with spatial interaction, which by itself would likely reduce the computational efficiency of the model. The model also includes simplified representations of the ground heat flux and snow surface temperature, which may be better captured by more process-based approaches. In particular, a more nuanced treatment of the ground heat flux may be desired in warmer snow climates (Mazurkiewicz et al., 2008).

Contrasts between modeled and observed snow metrics stem from several factors, including (but not limited to) uncertainties in climate forcings, SNOTEL-site-specific factors that the model neglects such as fine-scale topographic and vegetation patterns, and errors in model specification such as process representation and calibration. Despite these factors, errors at SNOTEL sites from the hourly SnowClim model run were relatively small and compared well with errors reported for other gridded snow products. Ikeda et al. (2021) evaluated the snow simulations from the same 4 km WRF model runs that we used to source our raw climate forcings (Rasmussen and Liu, 2017). Relative to SNOTEL sites, they found a −26.2 % bias in maxswe. In contrast, the SnowClim model achieves a maxswe bias of only 0.15 %. Wrzesien et al. (2018) compared maxswe at SNOTEL sites to maxswe from 9 km WRF simulations. Across sites, they found a correlation coefficient of 0.55 and a bias of −89 mm. SnowClim achieves a correlation coefficient of 0.94 and bias of −11 mm. In the Californian Sierra Nevada, Guan et al. (2013) blended modeled, remotely sensed, and observed data to capture SWE at six sites. Their method achieved a SWE RMSE of 205 mm compared to snow surveys. The SnowClim mean RMSE of daily SWE was 77 mm across all sites and 166 mm at Sierra Nevada sites. While errors at SNOTEL sites were generally low, the model did tend to overestimate maxswe and duration at some warm and wet sites and underestimate these metrics at dry sites (Fig. B2). Further evaluation of the parameters used here in more marginal snow environments would lend additional confidence to the application of SnowClim data in these areas. While the model's excellent performance relative to SNOTEL observations is in part due to the fact that the model was calibrated to SNOTEL data, the model could be calibrated to other observations for application in other contexts.

The flexible, modularized structure of the SnowClim model lends itself to calibration, parameter sensitivity assessment, and experimentation. In the western US, model performance was particularly sensitive to the choice of albedo algorithm and snow surface temperature parameterization, in line with previous findings (Etchevers et al., 2004; Günther et al., 2019; Slater et al., 2001; Fig. B3). Given the importance of impurities (e.g., tree litter, dust, and black carbon) on snow albedo and consequently snow melt (Waliser et al., 2011), a future step will be to add albedo algorithms that account for these effects. The modular structure of SnowClim would make this relatively straightforward.

Given the multifaceted importance of snow and ongoing snowpack changes due to climate change, there is a need for models that can accurately and efficiently simulate snow to generate spatially extensive, high-resolution datasets to meet the diverse requirements of different applications. We anticipate that the SnowClim model and data will be powerful tools for researchers and managers across a range of disciplines, including ecology and wildlife biology, recreation, transportation, hazard planning, and glacier and hydrologic modeling. The SnowClim model and data will be particularly useful for applications requiring high spatial resolution data, such as species distribution and refugia modeling, and will complement the existing array of snow models and datasets by providing a novel balance of process-based elements and computational efficiency.

Appendix A

We compared observed cold content and cold content simulated using SnowClim (both with and without the tax approach) at two sites in the Niwot Ridge Long-Term Ecological Research site (LTER), Colorado: an alpine site (Saddle) and a subalpine site (C1). The cold content observations are presented in Jennings et al. (2018).

For the comparison, we ran SnowClim at each of the sites using hourly observed climate forcings from the sites (Jennings et al., 2021) for water years 2001–2013. The model was first run with the cold content tax described in Sect. 2.2.6 (i.e., as reported in this paper). Following this, an additional run was performed in which the tax approach was turned off. The goal of this experiment was to evaluate how well the model simulates cold content relative to observations and to determine whether the tax approach improved cold content simulations compared to the same model without the tax approach. The same parameters from the western US run were used for modeling these sites, and the model was not recalibrated. Simulated cold content was smoothed using a 2-week moving mean to provide a more robust comparison to observations, which were collected at weekly to monthly intervals. Peak (minimum) cold content values were averaged across years to get a peak cold content value for each site and model run, and these were compared with the peak cold content values from observations reported in Jennings et al. (2018), Table 1.

Despite scale discrepancies between point observations and 210 m grid cells, SnowClim peak cold content compares well with observations when the tax approach is used. In contrast, when the tax approach is not used and all negative energy fluxes are directly added to the pack cold content, cold content becomes an order of magnitude greater than the observations.

Table A1 Comparison of peak seasonal cold content values from observations, SnowClim simulations including the tax approach described in Sect. 2.2.6, and SnowClim simulations without the tax approach at two sites from the Niwot Ridge Long-Term Ecological Research site in Colorado.

Appendix B

Figure B1 Map of the 3.1 million square kilometer modeling domain with locations modeled at 210 m spatial resolution shown in blue.

Figure B2 Performance of the best hourly model at SNOTEL sites in temperature–precipitation space. Each point represents a SNOTEL site.

Figure B3 Parameter sensitivity of hourly model performance.

Figure B4 Performance of the snow model without shallow snow correction for different time steps using the parameter set selected in calibration of the hourly model with shallow snow correction. Points represent median values across 170 SNOTEL sites.

Table B1 Summary of snow metrics by four-digit HUC. Total and modeled areas are in square kilometers, maxswe has units of millimeters, and snow duration has units of days. Values are averaged over all modeled grid cells within each HUC. Total snow metric values (in the last row) are averaged across all grid cells across the western US modeling domain, and total changes are computed from these total snow metric values.