Numerical simulations of snow water equivalent (SWE) in mountain systems can be biased, and few SWE observations have existed over large domains. New approaches for measuring SWE, like NASA’s ultra-high-resolution Airborne Snow Observatory (ASO), offer an opportunity to improve model estimates by providing a high-quality validation target. In this study, a computationally efficient snow data assimilation (DA) approach over the western United States at 1.75-km spatial resolution for water years (WYs) 2001–17 is presented. A local ensemble transform Kalman filter implemented as a batch smoother is used with the VIC hydrology model to assimilate the remotely sensed daily MODIS fractional snow-covered area (SCA). Validation of the high-resolution SWE estimates is done against ASO SWE data in the Tuolumne basin (California), Uncompahgre basin (Colorado), and Olympic Peninsula (Washington). Results indicate good performance in dry years and during melt, with DA reducing Tuolumne basin-average SWE percent differences from −68%, −92%, and −84% in open loop to 0.6%, 25%, and 3% after DA for WYs 2013–15, respectively, for ASO dates and spatial extent. DA also improved SWE percent difference over the Uncompahgre basin (−84% open loop, −65% DA) and Olympic Peninsula (26% open loop, −0.2% DA). However, in anomalously wet years DA underestimates SWE, likely due to an inadequate snow depletion curve parameterization. Despite potential shortcomings due to VIC model setup (e.g., water balance mode) or parameterization (snow depletion curve), the DA framework implemented in this study shows promise in overcoming some of these limitations and improving estimated SWE, in particular during drier years or at higher elevations, when most in situ observations cannot capture high-elevation snowpack due to lack of stations there.
Mountain systems with snowmelt-driven hydrologic cycles are crucial to more than one-sixth of the world’s population that depends on seasonal runoff as its primary water source (Barnett et al. 2005). Such systems include the western United States, the Alps, and the Hindu Kush, among others, many of which have agricultural and urban areas that depend heavily on the water supply provided by seasonal snowmelt. In addition, the timing of snow ablation plays an important role in soil and vegetation conditions in these mountainous regions, affecting ecosystems and phenomena such as wildfire timing, extent, and length of season (Westerling et al. 2006).
Because of the complexity of mountainous terrain, snow mass and snow extent exhibit great spatial heterogeneity, with their variability being influenced by a range of factors such as elevation, vegetation cover, slope and aspect, topography, orographic effects, wind redistribution, and atmospheric circulation (e.g., Dettinger et al. 2004; Clark et al. 2011; Girotto et al. 2014a). Given the importance and highly variable nature of the mountain snowpack, our ability to accurately estimate how much snow water equivalent (SWE) is stored in these mountainous regions is not only greatly desirable, but also a major scientific challenge (e.g., Dozier 2011; Margulis et al. 2015; Dozier et al. 2016; Painter et al. 2016). In situ observations, remotely sensed data, modeling, or a combination of these data sources have been used to address fundamental questions such as how much water is stored in mountain snowpacks, how does the snowpack vary in space and time, and how various physiological and climatological factors impact snow accumulation and ablation (e.g., Dozier 2011; Margulis et al. 2015). Answers to such questions still have significant uncertainty associated with them. This is due in great part to the fact that mountain snow is one of the most seasonally changing components of the Earth system, and there is presently no one single method that can accurately and fully capture the extent and heterogeneity of mountain snow and the processes that control the system dynamics.
In the western United States, in situ measurements such as snow pillows and snow courses have long been used by water resource managers to estimate seasonal water storage, yet these left runoff forecasts with significant errors (Dozier 2011). While the American West is a densely sampled mountain region, point-scale measurements cover only a very small portion of the snow-dominated areas. For example, in the Sierra Nevada the snow sensor density is 1 for every 700 km2 (Guan et al. 2013), with most sensors covering the lower to midelevations, 6500–9500 ft (~1900–2900 m), leaving the higher elevations undersampled (Dozier et al. 2016). This can be a problem in particular during drier years or toward the end of the snow season, when most snow is concentrated at the higher elevations and the snow has melted from many in situ measurement locations, as well as in a warming climate.
Remote sensing can provide useful measurements for mapping snow, including the undersampled high-elevation locations. There are two main methods of collecting snow information using spaceborne instruments. Passive microwave (PM) has been explored by past and current studies (e.g., Durand and Margulis 2007; Tedesco et al. 2010; Li et al. 2012), primarily because it can retrieve snow depth, which is a more direct and often more precise way of obtaining SWE, and it experiences no interference from clouds. However, PM data are often too coarse in spatial resolution—often greater than 10 km—making it difficult to use to estimate SWE given the subpixel spatial variability of snow and factors controlling snow and snowmelt. In addition, PM signal can saturate over deep snow. Studies have found this threshold to range from 100 to 150 mm, and to be as high as 200 mm (e.g., Takala et al. 2011; Liu et al. 2014). The other main technique is using visible and near-infrared (Vis/NIR; e.g., Painter et al. 2009; Rittger et al. 2013), which can provide information at better spatial resolution (30–500 m), but which is hindered by cloud cover, and perhaps more problematically, only retrieves snow-covered area (SCA) and not depth-integrated metrics. Because instantaneous snow cover does not provide ample information on instantaneous snow depth or SWE, using remotely sensed snow extent for water estimation requires more indirect methods such as reconstruction to determine SWE (e.g., Molotch and Margulis 2008; Schneider and Molotch 2016).
Recent advancements in airborne retrievals of snow depth, snow albedo, and SWE across watershed-scale regions have been providing unparalleled spatial distributions at relatively high resolution. The Airborne Snow Observatory (ASO; Painter et al. 2016) provides SWE estimates at 50-m spatial resolution, taken at approximately monthly intervals during accumulation and approximately weekly during ablation. ASO data are currently available over a growing number of watersheds in California and Colorado, with limited retrievals over the Olympic Peninsula in Washington State. Despite the high accuracy and rich spatial distribution of snow information it collects, at this time ASO does not yet provide these measurements over entire mountain ranges, at regional or subcontinental scales, or during the full duration of the water year (WY). In this study we use it as a validation dataset but recognize its value in informing current basin-scale operational water management practices, as well as its potential in data assimilation (DA) techniques.
Given the above challenges and shortcomings in observations when it comes to estimating SWE across mountainous terrain, using models is necessary to bridge temporal and spatial gaps in measurements. Snow models, as components in hydrologic, land surface, or climate models, are based on physical processes, and in combination with additional physiological data such as land cover (soil and vegetation type) and elevation maps, and forced with meteorological data, can provide spatially and temporally continuous SWE estimates at varying horizontal scales. With the advancement in computing resources in recent years, we can simulate SWE (and other hydrologic variables) at increasingly higher spatial resolutions. However, uncertainties in the input data, model parameters, and/or model design can lead to uncertainties in model output.
To address some of these deficiencies, a combination of observed data and modeling can be employed. DA is one such technique, which, generally speaking, uses observed data (e.g., remotely sensed) to correct or improve model estimates, and provides the associated uncertainty. This results in a spatially and temporally continuous SWE estimation that leverages the benefits of both observed data and modeling capabilities, which, when used over a relatively long period of time (e.g., multiyear), leads to the development of reanalysis datasets. Such datasets include the North American Regional Reanalysis (NARR; 32 km; Mesinger et al. 2006), NCEP–NCAR reanalysis (2.5°; Kanamitsu et al. 2002), and Modern-Era Retrospective Analysis for Research and Applications (MERRA; 50 km; Rienecker et al. 2011). While these existing data can be a useful resource for larger-scale regional or global climate studies, their relatively coarse spatial resolution makes them inadequate for subgrid, basin-scale SWE estimation across mountainous terrain. The Snow Data Assimilation System (SNODAS; Barrett 2003) does produce snow data over the United States and parts of southern Canada at a relatively high horizontal resolution of 1 km × 1 km, by using a forward modeling approach that is updated via nudging from a variety of remote sensing and in situ data. However, it has been shown to either underestimate or overestimate SWE when compared with other SWE estimation methods such as reconstruction, interpolation, or recent airborne SWE measurements (Dozier 2011; Rittger 2012; Bair et al. 2016; Dozier et al. 2016). Independent validation against ground-based observations is difficult since SNODAS assimilates the vast majority of these data as part of the SWE estimation process. Higher-resolution reanalyses have been recently developed, but they only cover certain river basins or mountain ranges (e.g., Sierra Nevada reanalysis; Guan et al. 2013; Rittger et al. 2013; Margulis et al. 2016). As a result, there are no high-resolution, spatially distributed, temporally and spatially continuous SWE datasets over mountainous terrain across regional to continental scales and spanning multiple years (with the exception of SNODAS, whose shortcomings are presented above).
While significant strides have been made in combining models and observations of snow in complex topographic terrain, accurately and comprehensively estimating spatially distributed SWE remains one of the biggest challenges in mountain hydrology (e.g., Dozier et al. 2016). One reason for this can be traced to the fact that in assimilation, most approaches integrate snow cover products into models to estimate SWE. Because an instantaneous SCA measurement does not inform how much instantaneous SWE there is in the snowpack, previous studies that used sequential filtering schemes such as ensemble Kalman filters (EnKFs) to assimilate snow cover data have only attained limited improvements in SWE estimates (e.g., Andreadis and Lettenmaier 2006; De Lannoy et al. 2012). If the goal of snow DA is to produce an SWE reanalysis—a multiyear retrospective snow dataset as opposed to short-term seasonal streamflow forecast—one can use a batch estimation approach (or smoothing), which is ideal for historic data (Dunne and Entekhabi 2005). In this method, the model is first run forward over the entire water year to estimate prior SCA, and then observed SCA values from throughout the seasonal cycle are assimilated into the hydrologic model in batch (all at once), such that estimated SWE at a given time is determined by including both prior and subsequent SCA observations. This is effective because in many cases SCA evolution throughout the snow season is highly correlated with SWE and in particular in the ablation period (e.g., Durand et al. 2008; Girotto et al. 2014a; Margulis et al. 2015). This method is different from the more typically used sequential filtering, where the state variable (e.g., SWE) is updated sequentially as (SCA) observations become available (and therefore the updates do not include information on the future state of the observed or modeled quantity). In other words, an advantage to using a batch smoother approach to assimilate SCA in order to update SWE is that it alleviates shortcomings related to how much information instantaneous SCA can providing for estimating instantaneous SWE (i.e., snowpack depth), by incorporating longer-term, seasonal SCA evolution behavior, which is more closely tied to SWE evolution. The batch smoother approach has been used successfully with assimilation of Landsat-derived fractional SCA data over the Sierra Nevada mountain range (Margulis et al. 2016). Because mountain snow changes rapidly during the ablation season, assimilating daily MODIS Vis/NIR—and in particular the fractional SCA (fSCA) MODIS Snow Covered Area and Grain Size (MODSCAG; Painter et al. 2009) product—using a batch smoothing approach could prove to be beneficial in improving SWE estimates across continental-scale domains.
Another challenge is the computational power and resources required to simulate and produce data at higher spatial resolution over extended periods of time and space. Members of the hydrologic community have acknowledged that in order to effectively tackle both research and applied water cycle challenges, higher-resolution (on the order of 1 km) hydrologic data are necessary across continental to global domains (e.g., Wood et al. 2011). As computational capabilities increase with the growth in computer power, hydrological models are run at increasingly fine resolutions and across larger spatial scales. These types of simulations require a DA approach that is computationally efficient and easy to understand mechanistically.
In this study, we test and validate a snow DA method that aims at simulating a continental-scale domain at 1.75-km resolution, to provide a spatially and temporally complete estimate of mountain SWE across the western United States over the past 17 years. The local ensemble transform Kalman filter (LETKF; Hunt et al. 2007), employed here as a batch smoother, is chosen for its relative simplicity, ease of implementation, and computational efficiency. The Variable Infiltration Capacity (VIC) hydrologic model (Liang et al. 1994) is used to compute the water and energy budgets and produces a first estimate of SCA and SWE. We then assimilate the MODSCAG (Painter et al. 2009) daily satellite SCA product to constrain an updated model SWE estimate. MODSCAG has a daily temporal resolution, provides subpixel (fractional) SCA, and offers the most accurate snow products available from MODIS (Rittger et al. 2013). One of the major novelties of this study is that our high-resolution SWE estimates (1.75 km) are validated against the ASO SWE product (50 m), the only ultra-high-resolution, spatially distributed SWE observations, which facilitates a more accurate assessment of modeled SWE spatial variability. In section 2 we describe the method for this study, including the hydrologic model, the data assimilation technique, the assimilation and validation datasets, and experimental design. In section 3 we present results, followed by discussion and summary in section 4.
a. Hydrologic model
The model used here is the VIC model (Liang et al. 1994), a robust and widely used macroscale hydrologic model that provides estimates of hydrologic state variables. The meteorological forcings used here to run VIC version 4.2.d include the PRISM daily precipitation and minimum and maximum temperature (http://prism.oregonstate.edu), and MERRA-2 wind (Gelaro et al. 2017; Global Modeling and Assimilation Office 2015; data can be accessed through https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/data_access/). These datasets were chosen because of their long temporal span, adequate spatial coverage over the western United States, daily temporal frequency, and in the case of PRISM precipitation, a high enough spatial resolution to capture complex mountain precipitation regimes.
In this realization we employ VIC in the water balance mode only, not full energy. This means that while the model does not solve the surface energy balance directly, allowing it to be computationally faster, in the snow module the surface energy balance is still computed explicitly in order to estimate the fluxes that drive snow accumulation and ablation (http://vic.readthedocs.io). The snow model, described in detail by Andreadis et al. (2009), accounts for snow on the ground as well as snow in the vegetation canopy. Ground snow has two layers, a thin surface layer and a pack layer with accumulation being controlled by snowfall and ablation that are calculated by solving the energy balance at the snow–air interface. Snow albedo is assumed to decay with age using an established parameterization, while densification results from compaction and snow metamorphism. Snow cover fraction (i.e., SCA) is calculated as a function of snow depth based on a probability distribution, being equal to 100% for depths greater than a critical depth (Cherkauer and Lettenmaier 2003). This relationship (SCA–snow depth) is often referred to as the snow depletion curve (SDC). The SDC algorithm in VIC was developed by Cherkauer and Lettenmaier (2003) based on deriving uniform probability distribution from observed snow depth at the University of Minnesota’s Rosemount Agricultural Experiment Station during the winters of 1997/98 and 1998/99. This fitted uniform distribution ignores extreme minimum and maximum observed values (P < 0.1 and P > 0.9) and is able to represent most of the spatial variability in daily snow depth observations, largely because there is little spatial variability in observations in this Cherkauer and Lettenmaier (2003) study, which is to be expected give the relatively flat terrain in Minnesota. In the parameterization, when the mean snow depth is greater than or equal to a critical depth, the snow coverage is 100% for that model grid. If the mean snow depth is less than the critical depth, snow-free areas exist in the grid. The algorithm can also account for cases when there is melting of thin accumulating snowpacks, or when there is fresh snow on a partially snow-free grid cell. More details about this algorithm can be found in Cherkauer and Lettenmaier (2003). It is important to note that because these parameterizations were derived based on snow depth observations over relatively flat surfaces, using uniform distribution, the resulting SDC is shallow, meaning once the snow is melting and reaches a certain depth, fSCA will quickly go to zero. Our western U.S. domain is mountainous, often with heterogeneous terrain and therefore snow depth, implying the current VIC SDC parameterization is not entirely appropriate for this study. Ideally, an SDC parameterization should be derived for respective mountain ranges or basins, since the snowpack is heterogeneous in mountain environments; however, this is an undertaking worthy of a study in and of itself, and was thus beyond the scope of this study. We hope to address this shortcoming in future work. In this instance of the model, the frozen soil option is set to “FALSE.” While frozen soil is an important consideration for springtime runoff and streamflow (e.g., Cherkauer et al. 2003) (from a water balance perspective), the focus of this study is on improving SWE estimates via DA rather than downstream components. Furthermore, frozen soil state can contribute to ground heat–snow fluxes, but these are often small compared to the snow–atmosphere interface flux exchange. Therefore, having frozen soil set to FALSE is not expected to have a significant effect on SWE accumulation or ablation processes.
b. Data assimilation approach: LETKF
Data assimilation techniques aim to estimate the state of a dynamic system by using both models and observations. There are several methods for approaching DA, ranging from direct insertion to a simple Kalman filter, extended Kalman filter, or EnKF, to batch smoothers like the Kalman-based ensemble batch smoother (EnBS) or a particle batch smoother. Examples of such DA work include Liston et al. (1999), Rodell and Houser (2004), Slater and Clark (2006), Andreadis and Lettenmaier (2006), De Lannoy et al. (2012), and Margulis et al. (2015). Ensemble-based assimilation techniques (e.g., Evensen 1994; Evensen 2003; Hamill 2006) are used for nonlinear dynamic systems and are often preferred since they can represent complex error structures without being computationally cumbersome (i.e., calculates and propagates error estimates from the variance of the ensemble of model states). In this study we use a particular version of an EnKF, the LETKF (Hunt et al. 2007), and we do so within the context of a batch smoother. This technique is relatively easy to implement, is computationally efficient, and has the flexibility to be adapted to a variety of applications while maintaining formal simplicity (of other DA methods). The LETKF DA method used here is adapted from an implementation in the Regional Hydrologic Extremes Assessment System (RHEAS) software (Andreadis et al. 2017), which also employs the VIC hydrologic model to produce hydrologic estimates.
The LETKF decomposes the estimation problem into a set of local domains (spatially or temporally discretized) where the analysis equations are solved independently using only nearby observations. Defining the background (i.e., open loop) state vector as xb (here, SWE), the predicted measurements (here, modeled SCA) as yb, and the observations as yo (here, observed SCA), we can derive the analysis (post-DA SWE) mean as
where N is the ensemble size (here set to equal 10); b is the matrix of background perturbations (i.e., a priori modeled SWE); b is the matrix of the predicted measurement perturbations (i.e., a priori modeled SCA); is the measurement (i.e., observed SCA) error covariance matrix; ρ is a multiplicative covariance inflation factor, here set to 1.05; and is the identity matrix. The ρ value is chosen based on recommendation of Hunt et al. (2007) and can be thought of as adjusting the influence of observations on the current analysis, such that the influence of an observation on future analyses decays exponentially with time, and thus localizes the analysis in time. Given that the assimilation is applied locally for each pixel, our relatively small ensemble size can capture the variability (e.g., Ott et al. 2004). The analysis (state) ensemble perturbations a (i.e., post-assimilation modeled SWE) can be updated as
and where a can be thought of as the optimal weights between model and observed states given model and observed errors. The ensembles of background b and predicted measurement b perturbations are randomly sampled from a 30-yr modeled SWE and SCA climatology (1981–2010) that is only generated once. This avoids producing meteorological forcing perturbation every time, making the algorithm more efficient. The observational error is computed by randomly sampling a Gaussian distribution with mean zero and standard deviation set based on the observed data used. Here we set a standard deviation observation error of 0.10, corresponding to MODSCAG data.
As mentioned in section 1, if SCA observations are assimilated in order to obtain retrospective SWE estimates, using a batch smoother method is more appropriate because it includes both prior and subsequent observations in the estimation of SWE at a given time, thus capturing the stronger SCA–SWE seasonal temporal correlation not often seen for instantaneous SCA–SWE relationships (which is what a sequential filter would capture). In other words, SCA DA methods that employ sequential filtering only use the observed SCA value at a given point in time (instantaneous SCA) to update (instantaneous) SWE, which typically limits the efficiency of this method to shallower snowpack. On the other hand, assimilating SCA in batch (using a batch smother) includes an entire snow season’s SCA observations to update SWE, taking advantage of the stronger seasonal correlation between SCA and SWE evolution. For this study, we use the LETKF as a batch smoother (Fig. 1): the model is first run forward in time for each water year (from 1 October to 1 October) to estimate prior SCA, and then observed SCA over a specified interval (an assimilation window) is assimilated into the hydrologic model in batch (all at once), retroactively updating SWE estimates for that period. More specifically, VIC is first run in open loop, estimating prior SCA and SWE. When the modeling system is run in the assimilation scenario (see section 2d), we restrict the assimilation window to 1 January–1 August of each water year in order to avoid spurious correlations during the beginning of the accumulation season. Outside of that window the model is run in open-loop mode, and assimilation of observed SCA is excluded. All observations within this assimilation window are incorporated at once, as a batch, and not sequentially. This approach has been used by other batch smoother DA studies, such as Girotto et al. (2014a,b). Because SCA is quite variable during the early accumulation season, going from 0% to 100% quickly, it can lead to positive, but spurious, correlations that would misleadingly update the weights in the filter (i.e., a). Furthermore, based on our sensitivity studies (not shown), the choice of length and time of the year for the assimilation window also matters: if too much of the accumulation season is excluded from the assimilation (e.g., assimilation window is 1 March–1 July or 1 February–1 August), the updated SWE will be underestimated at the end of the snow season because we would miss capturing part of the period during which SCA evolution is highly correlated with SWE. This sensitivity highlights the time of the year remote sensing snow cover observations are most needed.
c. Observational datasets
The data assimilation method described above uses remotely sensed fractional SCA to update modeled SWE. The observational dataset employed here is MODSCAG, a MODIS-derived product that uses the MODIS Snow Covered Area and Grain size retrieval algorithm (Painter et al. 2009) to produce a per-pixel fractional SCA and snow grain size from MOD09GA surface reflectance. In this work, we use the daily fractional SCA data, which range from 0% to 100% and account for any vegetation and soil within each pixel. MODSCAG fractional SCA has a lower root-mean-square error (RMSE) compared to MOD10A1 fractional cover (0.10 vs 0.23, respectively), consistently maintaining its retrieval capability during crucial transition periods of accumulation and melt, across all land cover classes and land surface types (Rittger et al. 2013).
To retrieve and process the MODSCAG data, the RHEAS software downloads daily geotiffs from the MODSCAG data portal (https://snow-data.jpl.nasa.gov/) via ftp for the time period indicated by the user, stitches the native ~500-m tiles into one daily geotiff, regrids the data to a user-specified spatial resolution (here 1 km2, in this case using the nearest neighbor method), and trims the data for a region of interest (here, the western United States). During assimilation the resampled 1-km2 MODSCAG data get reaggregated to a user-specified model resolution, here 1.75 km, by selecting the nearest MODSCAG 1-km2 pixel to the center of the model grid. At the time of this study only the noncanopy adjusted daily MODSCAG SCA data were available for the full length of DA model simulations (WYs 2001–17), which means the observed SCA used is the viewable rather than the ground observation, likely underestimating true SCA. Due to this factor, and given this is a fractional SCA which adjusts for vegetation cover area (but does not see under vegetation canopy), the MODSCAG fractional SCA over a model pixel can often range between zero and a value less than 100%. This becomes important when comparing to modeled SCA, which is generally either zero or 100% (due to a shallow snow depletion curve parameterization). In addition, while MODSCAG has a high daily temporal resolution, it is affected by the instrument’s viewing angle, causing snow cover measurements to often jump from full coverage to little to no coverage within the span of a few days. The implications of this model-observed snow cover discrepancy will be further discussed in sections 3 and 4.
2) Airborne Snow Observatory
Because this study aims to break new ground in terms of running at high spatial resolution over a large domain to produce better SWE estimates, our framework needs to use a high-resolution SWE dataset for validation. The one dataset that currently exists, and which can be exploited effectively for that use, is the NASA Airborne Snow Observatory SWE data.
The ASO is a coupled scanning lidar and imaging spectrometer observational airborne platform (Painter et al. 2016). It retrieves snow depth directly by flying over the region of interest during the snow-free season to capture the baseline topography, and then again when snow is present to measure the snow surface elevation. Differencing retrievals when snow is melted away and when snow is accumulated determines the snowpack depth. Using spatially distributed, physically based energy and mass balance models that are constrained by measured densities at the time of flight to estimate snow density, SWE can be deduced. ASO SWE product is available at 50-m spatial resolution, with a 5%–13% uncertainty, which is based on measured snow depth and modeled snow density uncertainties. ASO acquisitions are roughly weekly (weather permitting) during the ablation period, and only one snow-free retrieval is necessary. ASO can measure SWE under vegetation canopy.
In this study, we use the ASO SWE data product from several basins, including the Tuolumne basin in the Sierra Nevada (California), Uncompahgre basin in the Rocky Mountains (Colorado), and Olympic Peninsula in Washington State. Each basin has a varying number of observations and has different observational dates. Table 1 summarizes the available dates that were used in this study. The Tuolumne basin was the first basin observed by ASO, and therefore has the longest record, starting in 2013 and available each year to the present. This longer record is useful for assessing modeled SWE during a variety of snow years, including the lowest snow on record for the area (WY 2015), a roughly average snow year (WY 2016) and one of the highest snow years on record (WY 2017).
To have a more direct comparison between ASO and modeled SWE, ASO data are reprojected from the native transverse Mercator (m) projection to VIC’s geographic coordinate system (latitude–longitude; °), and regridded from the original 50-m resolution to the model’s 1.75 km by aggregating ASO pixels within each VIC grid using area averaging. The modeled data are trimmed to the ASO basin shape (e.g., see Fig. 5).
3) In situ snow measurements (SNOTEL)
Previous modeling studies have used point measurements for validation (e.g., Guan et al. 2013; Rittger et al. 2016), so we include them here for completeness. However, we note the limitations that come with comparing them against a spatially distributed dataset. In situ measurements such as snow pillows do not explicitly or implicitly capture SWE spatial variability because of their location, or lack thereof—snow pillows do not exist at the highest elevations and are often placed on flat ground, which may not represent the surrounding complex terrain (Dozier et al. 2016). These types of observations also present a scale discrepancy, since a “point” measurement cannot fully capture SWE amount over an area (like the model grid space). In addition, interpolation of such point-based data can have large uncertainties (e.g., Fassnacht et al. 2003; Dozier 2011) or insufficiently explain SWE variability (e.g., Molotch and Bales 2005). Because of these reasons, we focus our validation on the more direct comparison of modeled data to ASO SWE observations, especially when assessing spatial variability, but do include snow pillow data for reference, and because its daily resolution allows the evaluation of the modeled SWE temporal evolution. In this study, we use snow course and snow pillow [i.e., Snowpack Telemetry (SNOTEL)] data from the National Resources Conservation Service and National Water and Climate Center website (https://www.wcc.nrcs.usda.gov/snow/snow_map.html) for the Tuolumne basin, averaging over the available stations within the validation domain.
d. Experimental design
The VIC model is run over the western U.S. domain (see Fig. 3) at a spatial resolution of 1.75 km (3 km2) and a daily time step. This domain is composed of five Hydrologic Unit Code 2 (HUC-2) areas (the regional-level USGS-defined hydrologic units): California, Upper Colorado, Lower Colorado, Great Basin, and Pacific Northwest, covering an area roughly bounded between 31° and 53°N and between 125° and 105°W. A significant portion of this region is characterized by complex mountainous terrain, with three main mountain ranges: the Sierra Nevada, the Cascades, and the Rocky Mountains, each defined by their own climate characteristics. Two datasets are produced over the western United States, a 37-yr open-loop simulation covering WYs from 1981 to present (1 October 1980–1 August 2017) where no observed data are assimilated, and a 17-yr assimilation simulation for WYs from 2001 to present (1 January 2000–1 August 2017), where MODSCAG SCA is assimilated during the 1 January–1 August window only, for each year (with the rest of the water year being run in open loop; see section 2b for more in-depth discussion of implementation of the assimilation window). In both cases we initialize the model at the start date (1 October 1980 for open loop and 1 January 2000 for DA) and use the restart functionality [at 1 October of each year in the open-loop case, and at the beginning (1 January) and end (1 August) of the assimilation window in the DA case], to provide continuous, consistent energy and water balance simulations.
To validate the snow DA framework employed here, we take advantage of the ASO observational data. Here we focus most of our model validation over the Tuolumne basin in California’s Sierra Nevada during the recent severe drought and recovery the region has experienced, WYs 2013–17 (a note that WY2017 is defined here as 1 January 2016–1 August 2017 due to the timing of forcing data availability, the writing of the study, and time needed to rerun WY2017). The Tuolumne basin is the main source of water for the city of San Francisco and other Bay Area communities and is also a source of hydroelectric power. The watershed reaches elevations as high as 3998 m (~13 100 ft) above sea level and is characterized by granite bedrock and steep slopes in the highest parts of the watershed, as well as alpine meadows and forested areas at mid- and lower elevations. ASO retrievals over the Uncompahgre basin in Colorado’s Rocky Mountains and the Olympic Peninsula in Washington are also used in the validation process, assessing model ability across a variety of regional climates and mountain ranges. The location and general coverage of these basins are marked in Fig. 3d with red boxes.
a. Overall results in western United States
Figure 2 shows western U.S. domain-averaged SWE for WYs 2001–17. Comparing the open-loop and assimilation SWE results, in general the assimilation of MODSCAG fractional snow cover lowers modeled SWE. This pattern is evident across all mountain ranges for springtime (April) modeled SWE (Fig. 3), with the exception of the Sierra Nevada, where some areas see an increase in SWE post-assimilation. This difference is largest over deeper snowpacks, such as those in the Cascade Mountains. The main reason for this is that MODSCAG SCA for a given pixel ranges from zero to a value often less than 100% (as detailed in section 2c). By contrast, modeled snow cover area is essentially binary, either 0% or 100% (due to a shallow snow depletion curve parameterization). (Here we refer to the value of SCA as being mostly binary because it changes from 0% to 100% relatively quickly; it does not represent a flag of snow cover existing or not existing.) In deeper snow, modeled SCA will be 100%, therefore adjusting VIC snow cover with MODSCAG fractional SCA observations will generally tend to lower post-assimilation SWE estimates. Exceptions to this occur during extreme snow years or late in the melt season. For example, during winter and spring of WY 2015 when drought dominated most of the American West (e.g., https://www.ncdc.noaa.gov/sotc/drought/201504), snow cover was lower, and so modeled SCA was mostly 0% (as opposed to 100%), whereas observations varied from 0% to 100%, therefore causing post-assimilation SWE to increase. This kind of impact from assimilating MODSCAG can be further understood by looking at Fig. 4, which shows maps of modeled and observed snow cover for 1 February and 1 May 2013. During the middle of the snow season (e.g., 1 February 2013) when snow coverage tends to be greater, the model shows a value of 100% everywhere there is snow, while MODSCAG SCA varies between 0% and 100% (Figs. 4a,c), therefore leading to a lower post-assimilation modeled SWE. One the other hand, toward the end of the melt season (e.g., 1 May 2013), or when snow coverage is low, the model will have much less snow (Figs. 4b,d), and inclusion of MODSCAG will tend to increase post-assimilation modeled SWE. (The impact to SWE can be seen in Fig. 7a below.) This behavior is further discussed in the next section.
b. California: Tuolumne basin, Sierra Nevada
To evaluate whether assimilating MODSCAG SCA using the LETKF batch smoother method is effective in updating SWE, both open-loop and post-assimilation SWE are compared to ASO-observed SWE. ASO retrieved spatially distributed snow observations during both the lowest snowpack (2015) and one of the largest snowpack (2017) years seen in California, as well as an average water year (2016), allowing for an unusually comprehensive evaluation of model results and snow DA technique across space and time and across a wide range of snow conditions. Figure 5 shows modeled and observed springtime SWE (mm) in the Tuolumne basin (California), averaged over ASO available dates for each given water year (see Table 1), which are mostly during the spring snow ablation period.
In general, the model (open loop) tends to underestimate SWE with respect to ASO observations (Figs. 5, 6), anywhere from a −36% to −92% difference, but assimilation of MODSCAG reduces that bias (Table 2). In drier years such as WYs 2013–15, assimilation of MODSCAG SCA adds snow (Fig. 5), reducing the basin bias from −68%, −92%, and −84% in the open-loop scenario to 0.6%, 25%, and 3%, respectively (Fig. 6). Under these drought conditions, when the snowpack is already low, the open-loop model scenario removes snow too quickly, underestimating peak SWE and reaching complete melt-out early compared to ASO retrievals (Figs. 7a–c). By adjusting the model snow to a higher value, DA causes the model to maintain a relatively deeper pack (compared to open loop) for longer (Figs. 7a–c), more closely matching ASO observations (albeit with a small positive bias). This DA-model behavior can be seen in Figs. 4b and 4d where modeled and MODSCAG observed spatial SCA are plotted for a melting season date, 1 May 2013. The open-loop modeled SCA is almost zero basinwide, with a 5% average snow cover, while MODSCAG basin-average snow cover is 18% with greater snow distribution. The inclusion of these observed MODSCAG values through assimilation nudges the model closer to observations, as seen in Fig. 7a. This improvement during drier years occurs at all elevations, with DA assimilation increasing SWE magnitude throughout all elevation zones (Fig. 8), and more importantly at higher ones (above 2400 m, i.e., ~8000 ft) where most snow tends to exist during drought years. The RMSE metric also reflects an improvement in modeled SWE via DA during these years (Table 2), with the DA scenario having a lower RMSE throughout most of the Tuolumne basin (Fig. 6) and as a basin average, and in particular at mid- to high elevations and along ridgelines. One exception is WY 2015, where DA had a higher basin-average RMSE value compared to open loop, despite the fact that basin-average SWE bias is greatly improved in the DA case. It appears that while DA increases SWE estimates in WY 2015 over open loop, generally improving SWE estimates, in some areas of the basin it increases a bit too much, leading to a slightly higher RMSE compared to open loop (Fig. 6).
In an average year (WY 2016), data assimilation reduces the negative bias (Fig. 6) relative to ASO by 28% (Table 2). During the late melt season (Fig. 7d) the improvement after DA is even greater, the percent difference being −27% compared to open loop’s −96%. The DA allows the model to capture the final melt-out during June 2016 compared to open loop, which removes the snowpack about a month earlier. During peak SWE, however, the model underestimates snow in both scenarios, with the DA case performing slightly better (−55%) than the open loop (−74%) (Fig. 7d). While DA adds some snow throughout all elevation zones (Fig. 8), it is not enough compared to ASO observations. This deficiency is likely due to the snow cover–SWE relationship parameterized in the model through the snow depletion curve. Overall, DA does improve SWE estimates throughout most of the Tuolumne basin compared to the open loop by increasing SWE and decreasing SWE bias and RMSE (Fig. 6, Table 2).
During an anomalously wet year (WY 2017), open-loop SWE has a −36% change bias (−313 mm), with assimilation of MODSCAG SCA leading to a greater SWE bias, a basinwide −49% difference (−423 mm; Fig. 6; Table 2). This lack of improvement in WY 2017 by the DA scheme is also reflected in the RMSE metrics (Table 2) and can be seen throughout all elevation zones (Fig. 8). While open-loop simulations tend to underestimate SWE with respect to observations (both in situ SNOTEL and ASO) in all years modeled here, particularly during peak snow, the difference in WY2017 is that DA does not lead to an increase in SWE as it does in other WYs (and therefore an improvement in the negative bias); rather, it further increases the negative bias with respect to observations. In this anomalously wet scenario in 2017, during the main part of the snow season (from January to mid-May), the snow depletion curve keeps VIC (open loop) modeled SCA at 100%; in other words, the snow depth is well above the rather shallow critical depth and SCA is then at 100%. On the other hand, MODSCAG SCA varies between 0% and 100% and tends to be less than 100% even when the snowpack is that large, due to vegetation, rocks, and instrument viewing angle lowering the overall pixel SCA (e.g., Figs. 4a,c exemplify this behavior). The difference in SCA between model (open loop) and assimilated data therefore leads to a lower post-DA SWE, both compared to open loop, and to ASO observations. Late in the season, during peak ablation (from mid-May through June) when snow depth is lower, the shallow snow depletion curve causes open-loop modeled SCA to jump to low SCA (near zero), and therefore assimilating MODSCAG leads to a post-DA SWE slightly greater than open loop (e.g., Figs. 4b,d exemplify this behavior), but still not reaching observed ASO values. In other words, during peak accumulation season when the snowpack is deep and snow coverage is relatively large, modeled SCA is already at 100% for a significant period of time (albeit biased low) and DA cannot further increase that, while during melt season DA, despite increasing SWE over open loop, is simply not able to compensate for the open-loop biases. The model SWE underestimation in this work is consistent with other current SWE model estimates (e.g., Raleigh et al. 2015, Broxton et al. 2016) and likely includes contributions from precipitation forcing biases; given that 2017 was an anomalously wet year, it is possible the precipitation gauges that inform the model forcing under capture the full extent of the snowfall, causing a negative precipitation forcing bias and an underestimated open-loop SWE.
Overall, when compared to spatially distributed ASO SWE observations over the Tuolumne basin in the Sierra Nevada, assimilating observed SCA improves modeled SWE when the snowpack is low, such as during dry years or toward the end of snow season, which is significant given this is a time when snow estimates are most uncertain given the lack of in situ point-based measurements at high elevations. Figure 8 supports this finding, showing DA increases snow at all elevation bands, but particularly important above 2400 m (~8000 ft). Having the ASO spatially distributed SWE data has enabled a more direct evaluation of model and DA performance in a spatial context. Temporally, DA improves correlation with respect to the daily SNOTEL measurements in the Tuolumne for all years, but most notably during WYs 2014–16, with open-loop Pearson correlation coefficients of 0.76, 0.78, and 0.85 and post-assimilation coefficients of 0.87, 0.89, and 0.91, respectively. If we consider the anomaly correlation coefficient, which calculates correlations with respect to climatology, DA improves SWE estimates for WYs 2015 (0.98 for DA and 0.95 in open loop) and 2017 (0.98 for DA and 0.96 in open loop). Temporal anomaly correlation is lower in the DA case for WYs 2013 (0.78 for DA, 0.94 for open loop), 2014 (0.77 for DA, 0.95 for open loop), and 2016 (0.44 for DA, 0.56 for open loop). [We calculated the anomaly correlation coefficients by computing the daily modeled SWE anomaly (for both open loop and DA) with respect to the corresponding monthly modeled climatology, the daily observed (SNOTEL) anomaly with respect to the corresponding monthly observed climatology, then computing the correlation coefficient of these anomalies (e.g., r(SWE′DA,WY2013, SWE′SNOTEL,WY2013). We used years 2001–16 to calculate monthly climatologies, since that was the period for which both open-loop and DA SWE were simulated.] It is important to note here that the improvements seen by assimilating MODSCAG are over the VIC model open loop, which perhaps produces suboptimal SWE estimates given that the model was run in water balance mode and not full energy mode. Employing the VIC model in its water balance mode means the skin temperature to air temperature and frozen soil parameterization is turned off, which can affect the accuracy of the snowpack simulation (implicitly via the ground heat flux). Nevertheless, despite these limitations of the open-loop simulations, assimilation of MODSGAC actually overcomes most of these model physics simplifications, thus adding value and the potential for improvement in estimating SWE.
c. Washington: Olympic Peninsula
Over the Pacific Northwest’s Olympic Peninsula, for the two dates in 2016 that ASO data are available, the model has a relatively low bias compared to observed SWE, with the open loop having a domain-average percent difference of 26%, while assimilation decreases that bias to a nearly 0% change (Table 2), albeit the biases have more spatial variability in their sign relative to the other basins studies. The decrease in the RMSE score in the DA case, both spatially throughout the basin (Fig. 9) and as a basin average (Table 2), supports the argument for DA improving SWE simulation in the Olympic Peninsula. Spatially, the model is able to capture the SWE distribution across this complex terrain (Figs. 9a–c). The model–ASO difference maps show an SWE bias pattern corresponding to the rain shadow typical of this region, with the model underestimating SWE in the wetter regions of the peninsula (southwest part), and overestimating in the drier areas (northeast part), to varying degree for open loop and assimilation (Figs. 9d,e). Comparing the model cases in Fig. 9f, DA tends to add snow in the midelevation regions, where snow might be more variable, and removes snow in areas where the snowpack is deeper (higher elevations), that is, where model snow cover tends to be 100% more consistently throughout the season. This behavior is consistent with the analysis over the Tuolumne basin in California, where incorporating observed SCA nudges the model toward a higher SWE value in drier, more variable snow conditions (or places), and toward a lower snow for deeper snowpacks. While it is difficult to draw any robust conclusions about model and DA performance given there are only two ASO measurements, Fig. 11a displays the temporal variability of modeled SWE in the Olympic Peninsula mountain range along with ASO observations, for completeness.
d. Colorado: Uncompahgre basin, Rocky Mountains
As with the other regions of the western United States, the model (open loop) underestimates SWE across the Uncompahgre basin in the Colorado Rockies. However, assimilating observed SCA does improve the spatial distribution of snow, as well as SWE magnitude (Fig. 10), reducing basin bias by 19% and decreasing the RMSE (Table 2). DA increases SWE across the higher model elevation bands in the southern part of the basin sampled here (Fig. 10f), bringing snow distribution closer to that observed by ASO. This is also reflected in the lower RMSE values seen in that part of the basin in the DA scenario compared to the open loop (Figs. 10g,h). With only two ASO observational dates available over the Uncompahgre, it is difficult to draw any significant conclusions, especially about the temporal evolution of the modeled snow, but as Fig. 11b suggests, DA improves late spring SWE values, keeping a deeper snowpack compared to the open loop. This too is consistent with the previous basins, with DA acting to increase SWE for conditions when model snowpack is lower. Just as we have seen in the Tuolumne basin, toward the end of the melt season, informing the model with observed MODSCAG SCA leads to a relatively deeper snowpack, bringing it closer to ASO SWE measurements (Fig. 11b).
It is also important to note that some of the spatial SWE differences between the model and observations for any of the regions analyzed could be due to differing topographic maps, that is, certain features such as high ridges, deep valleys, or peaks might not line up exactly in the model and ASO data, as well as differing elevation maps, that is, the model does not capture the highest of the elevation zones.
4. Summary and discussion
The capability to estimate time- and space-continuous SWE in mountainous terrain with snowmelt-dominated systems is greatly needed, but remains a significant challenge. In situ observations can provide direct SWE information throughout the water year, but fall short on capturing SWE spatial distribution. Observing the snow mountain systems remotely presents some opportunities, but clouds, temporal frequency or spatial resolution, and satellite instrument viewing angle can interfere with reliable, complete snow retrievals. Data assimilation, which combines modeling and observations, is one approach to achieve a spatially and temporally complete SWE state estimate. In this study, we have expanded upon previous efforts by implementing and testing a DA method for long-term snow estimation over a large domain at a high spatial resolution (1.75 km). We employ the VIC hydrologic model and the LETKF as a batch smoother over the western U.S. region and assimilate MODSCAG SCA to obtain SWE estimates, which we then validate against the unique ASO data. This is the first time a snow DA technique is assessed against spatially distributed SWE measurements, which provides a more reliable way to evaluate model performance when it comes to snow distribution.
General conclusions taken from our results indicate that our DA approach performs best under relatively low(er) snow conditions, whether that is a dry year or the end of a melt season. Post-assimilation Tuolumne basin-average model bias with respect to ASO is lowest in 2013 (0.6%) and 2015 (3%), in the midst of the California drought. In general, the DA SWE had a lower bias with respect to ASO than open loop, for all years and all basins analyzed, with the exception of Tuolumne basin in 2017, which is one of the wettest years on record. This different model behavior for low versus high snow conditions is also evident in 2016 over the Tuolumne. During peak SWE and early snowmelt (April–May) the model post-DA underestimates SWE, but at the end of the melt season (June) the data-assimilated SWE matches observations very closely both spatially and temporally, with a basin-average percent difference of −18%. These results are generally speaking consistent with previous Kalman-based ensemble batch smoother (EnBS) studies, where the DA improved SWE estimates significantly over open loop, especially for drier years, and sees a more modest improvement during wet years (Margulis et al. 2015), when as mentioned earlier, the DA leads to a smaller update when prior model bias is significant. For drier years, the LETKF batch smoother biases with respect to observations in this study are comparable (~13 mm) to those seen in the EnBS study [~200 mm estimated from Figs. 3–5 in Margulis et al. (2015)]. Of course, it is important to note that there are several differences between these studies aside from the batch smoother used, such as observational datasets used for validation (spatial SWE from ASO in current work and point-based snow pillows in EnBS study), the time period and domain of the simulations (albeit both looked at basins in the Sierra Nevada), model resolution, and assimilated data, which would contribute to some of the differences we see in the respective results. Temporally, we found DA improves SWE estimates and time evolution for all (five) water years when compared to daily snow pillow data, with Pearson correlation coefficients ranging 0.87 (WY 2014) to 0.97 (WY 2017) over Tuolumne. This is comparable to the temporal Pearson correlations of 0.89 achieved by other EnBS studies (e.g., Margulis et al. 2015) when compared to snow pillow and snow course data in a basin in the Sierra Nevada. In general, assimilation of MODSCAG reduces the RMSE in all the basins we looked at in this study, further supporting the added values of employing such a DA framework and technique, with one exception being the extreme snow years in the Tuolumne basin, WYs 2015 (dry) and 2017 (wet).
By comparison against the spatially distributed ASO observed SWE, we show the model is generally able to correctly capture snow spatial distribution in both open-loop and assimilation cases—at least with respect to elevation, one of the main factors influencing snow distribution. During dry and average snow years, while the open loop correctly reproduces snow distribution with elevation, it fails to capture the SWE magnitude throughout all elevation bands. Constraining the model with MODSCAG SCA improves SWE magnitude estimates during these years. This is an encouraging result given that most in situ snow observations do not capture high-elevation snowpack due to lack of stations at these locations (e.g., Dozier et al. 2016), yet this is where most of the snow exists at the end of the snow season, or during drought years. In abnormally wet years, when there is full snow coverage for an extended amount of time, and when the snowpack is deep, assimilation acts to increase the bias by lowering SWE. Our results suggest that this is because snow cover parameterization in the model allows for snow cover to have values of either 0% or 100%, while MODSCAG snow cover, being a fractional SCA product that accounts for vegetation and soil, has snow cover values ranging from 0% to 100%, most often not reaching or staying at full coverage. Thus, during dry periods when the snowpack is smaller, MODSCAG SCA is larger than the open-loop modeled snow cover, leading to a greater SWE post-assimilation, whereas during wet years like 2017, modeled snow is at 100% coverage, and MODSCAG SCA would cause a decrease in post-assimilation SWE. It is important to note that while DA in this case can lead to a decrease in SWE for pixels with deeper snow, there are locations where DA appropriately adds snow where snow may have not existed at all in open-loop simulation, or there was too little snow. An example of this is the 2016 Olympic Peninsula snowpack, decreasing the model bias with respect to ASO. The same is true in Colorado’s Uncompahgre basin, where assimilation of MODSCAG SCA resulted in a lowering of the negative model bias by increasing SWE for most grids in the domain (Fig. 10f).
Looking at the results for the western United States, DA generally tends to lower SWE compared to open loop, which at first might appear contradictory. However, these are consistent with the validation basins results. Areas that have deeper, more consistently snow-covered mountains (e.g., Cascade Mountains, northern U.S. Rocky Mountains) see a reduction in SWE after DA, given the model versus MODSCAG SCA behavior explained above (binary vs fractional). However, there are regions or times where DA adds snow, for example in the Sierra Nevada in April (e.g., Fig. 3, blue areas), or during dry years or late in the season (e.g., Fig. 2, years 2009, 2013, 2015, and 2016), consistent with ASO basin validation results. Because we validate with ASO, which is mostly available over the Sierra Nevada, we end up assessing the model and DA behavior over what might be an exception to the western United States as a whole rule, that is, the Sierra Nevada are different from the rest of the western U.S. mountain ranges. As for the western U.S. behavior, the fact that DA underestimates SWE over this area has also been found by other studies, where it was concluded that most reanalysis (DA) datasets have too much snow ablation at near-freezing temperatures (e.g., Broxton et al. 2016).
The data assimilation method implemented and tested in this study shows promise, in particular during critical times such as drought years or at the end of the snow season. However, the model, even prior to DA, shows a negative SWE bias with respect to observations, not capturing peak SWE during average to wet water years. This is likely due to an underestimated precipitation forcing, an inadequate snow cover–SWE parameterization in the model, and/or limitations from running VIC in water balance mode. For future work, improving the representation of the relationship between SWE and SCA (i.e., the snow depletion curve) in the hydrologic model should improve SWE estimates. Furthermore, the different behavior of DA results for low versus high snow conditions, rising in part from the discrepancy between modeled SCA estimation (0% or 100%) and observed fractional snow cover (ranges from 0% to 100%, but often less than 100%), points to the fact that if we want to take full advantage of the fractional SCA information MODSCAG provides (which captures subpixel snow variability and impact of vegetation and rocks), the model needs to be improved to also capture that subpixel snow distribution, even when the model is run at a relatively high resolution such as 1.75 km × 1.75 km. Another improvement in post-assimilation SWE estimates could come from adjusting the assimilation window based on length of snow season. In other words, instead of performing assimilation during 1 January–1 August for each water year, the end date for assimilation should be more dynamic, ending at about the time of the last snow.
Despite the potential shortcomings due to the VIC model setup (e.g., water balance mode instead of full energy; precipitation input) or parameterization (snow depletion curve), the DA method implemented and tested in this study shows promise in overcoming some of these limitations and improving estimated SWE using statistical methods, in particular during drier years or at higher elevations, that is, in situations when snowpack would be lower. Here we have described developing a modeling and DA framework that is suitable for high spatial and temporal resolution and over a large spatial extent while remaining computationally effective, which meant making some perhaps less than ideal choices when it comes to snow and surface energy and physics modeling (e.g., by running the VIC hydrologic model in water balance mode as opposed to full energy). To improve upon the promising framework and resulting snow estimates shown in this study, model settings such as SDC, precipitation input, and water versus energy balance mode need to be further tested through sensitivity studies in future efforts, and adjusted accordingly.
Connecting SCA and SWE, in particular when assimilating SCA to update SWE (i.e., the state estimate is not the same as the observed assimilated estimate), remains a challenge because the relationship between these two variables is not well represented in models (nor is it well defined in the observations, as the processes controlling SCA and SWE are different). This highlights the need to assimilate observed SWE where present in order to compute SWE estimates that are continuous in space and time over continental scale. The modeling and DA framework developed here for the western United States provides a path forward for assimilating such data, at a high spatial resolution, and continuous in time. A next step would be to assimilate SWE or snow depth measurements from the ASO program, which is expanding to the entire Sierra Nevada of California and into the Colorado River basin and Rio Grande basin. The ultimate goal is to assimilate remotely sensed SWE from spaceborne instruments, which would cover a greater area than airborne instruments can, and ASO can serve as pathfinding for spaceborne missions with the promise of its operational version providing robust calibration of the spaceborne retrievals.
This work was supported by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the U.S. National Aeronautics and Space Administration. We thank the three anonymous reviewers for their constructive feedback and suggestions, which improved this paper.
Current affiliation: Raytheon Company, Pasadena, California.
Current affiliation: Department of Civil and Environmental Engineering, University of Massachusetts Amherst, Amherst, Massachusetts.
Current affiliation: Global Institute for Water Security, University of Saskatchewan, Saskatoon, Saskatchewan, Canada.