Assimilating satellite-based snow depth and snow cover products for improving snow predictions in Alaska

Advances in Water Resources 54 (2013) 208–227

Contents lists available at SciVerse ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Assimilating satellite-based snow depth and snow cover products for improving snow predictions in Alaska Yuqiong Liu a,b,⇑, Christa D. Peters-Lidard b, Sujay Kumar a,c, James L. Foster b, Michael Shaw b,c,d, Yudong Tian a,b, Gregory M. Fall e a

Earth System Science Interdisciplinary Center, The University of Maryland, College Park, MD, USA Hydrological Sciences Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD, USA Science Applications International Corporation, Beltsville, MD, USA d Air Force Weather Agency, Offutt, NE, USA e National Operational Hydrologic Remote Sensing Center, 1735 Lake Drive West, Chanhassem, MN, USA b c

a r t i c l e

i n f o

Article history: Received 30 August 2012 Received in revised form 1 February 2013 Accepted 12 February 2013 Available online 27 February 2013 Keywords: Data assimilation Satellite snow products Snow prediction Streamflow prediction

a b s t r a c t Several satellite-based snow products are assimilated, both separately and jointly, into the Noah land surface model for improving snow prediction in Alaska. These include the standard and interpreted versions of snow cover fraction (SCF) data from the Moderate-Resolution Imaging Spectroradiometer (MODIS) and the snow depth (SD) estimates from the Advanced Microwave Scanning Radiometer for the Earth Observing System (AMSR-E). The satellite-based SD estimates are adjusted against in situ observations via statistical interpolation to reduce the potentially large biases, prior to being assimilated using an ensemble Kalman filter. A customized, rule-based direct insertion approach is developed to assimilate the two SCF datasets. Our results indicate that considerable overall improvement on snow prediction can be achieved via assimilating the bias-adjusted satellite SD estimates; however, the improvement does not always translate into improvements in streamflow prediction. Assimilating the standard MODIS SCF is found to have little impact on snow and streamflow predictions, while assimilating the interpreted SCF estimates, which have reduced cloud coverage and improved snow mapping accuracy, has resulted in the most consistent improvements on snow and streamflow predictions across the study domain. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Owing to its unique characteristics of high albedo, low roughness, and large heat capacity, snowpack plays an important role in modulating the hydrologic cycle across multiple spatiotemporal scales and represents a critical driver for floods, droughts, and many other applications in snow-dominated areas [1,2]. It is hence essential to be able to accurately predict snow accumulation and ablation processes to support the various applications in these areas. Like the study of other environmental processes, the modeling and prediction of snow can be subject to substantial uncertainty arising from errors in the forcing data and the deficiencies in the physical representations of the snow model. Different strategies have been developed or adopted to improve snow simulation, including, for example, enhancing model physics to accommodate

⇑ Corresponding author at: Earth System Science Interdisciplinary Center, The University of Maryland, College Park, MD, USA. E-mail address: [email protected] (Y. Liu). 0309-1708/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.advwatres.2013.02.005

complex environments such as forest areas (e.g., [3–6]), evaluating multiple snow models to account for model structure uncertainty (e.g., [7,8]), optimizing model parameters to minimize parametric uncertainty (e.g., [9–11]), or using reconstruction techniques based on snow disappearance timing and back-calculated snowmelt to reduce the influence of precipitation uncertainty (e.g., [12]). With the increasing availability of in situ and remotely-sensed snow observations in recent years, a great deal of research has also examined the potential of assimilating observations of snow quantities, such as snow cover fraction (SCF), snow water equivalent (SWE), and snow depth (SD) to improve the prediction of snow accumulation and ablation, as well as snowmelt runoff (e.g., [13– 21]). For large-scale snow assimilation applications, the success of assimilation critically depends on the use of an appropriate snow prediction model, high quality snow observations, properly specified error characteristics (for all relevant uncertainty sources), and an effective assimilation algorithm. From a modeling perspective, apart from improving physical parameterizations, it is important to adequately represent the spatial variability of snow cover in the models to be able to accurately predict the magnitude, timing,

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

and duration of snowmelt and the resulting runoff in snow-impacted regions (e.g., [22]). Aside from emphasizing the importance of sub-grid parameterization in snow modeling, this also suggests the potential of using high-resolution models (e.g., on the order of 1 km) with spatially distributed parameters to capture the effect of small-scale variability in the physiographic and topographic characteristics of the land surface on snowpack evolution (e.g., [11,23]). Nevertheless, despite the great deal of research into characterizing the spatial variability in snow distribution (e.g., [22]), large-scale high resolution modeling of snow accumulation and ablation still remain largely elusive. From the observing perspective, daily or hourly in situ snow measurements are available from many weather or snow monitoring stations across the continental United States, such as the National Oceanic and Atmospheric Administration (NOAA) Cooperative Observer Program (COOP) station network (http:// www.nws.noaa.gov/om/coop), and the Natural Resources Conservation Service (NRCS) SNOwpack Telemetry network (SNOTEL, http://www.wcc.nrcs.usda.gov/snow/). However, given the high spatial variability in snow distribution, these point-scale measurements alone may not be adequate for large-scale assimilation applications, especially when the in situ observing network is sparse. This highlights the need to integrate spatially-distributed snow information from, e.g., global snow estimates from remote sensing instruments onboard earth observing satellites (e.g., [24]). Examples of satellite-based snow observations include the SCF estimates from the visible or near-infrared sensors such as the Advanced Very High Resolution Radiometer (AVHRR), the Landsat Thematic Mapper (TM), and the Moderate Resolution Imaging Spectroradiometer (MODIS). In addition, SWE estimates are also available from microwave sensors like the Scanning Multichannel Microwave Radiometer (SSMR), the Special Sensor Microwave/Imager (SSM/I) and the Advanced Microwave Scanning Radiometer for the Earth Observing System (AMSR-E). The global terrestrial water storage (TWS) estimates derived from the observations of the Gravity Recovery and Climate Experiment (GRACE) satellite can also be used to infer snow accumulation for cold regions during winter time (e.g., [16,25]). It is important to understand that each of these sensors has its own strengths and limitations. For example, although their accuracy may vary by land cover and snow conditions (e.g., [26]), snow products based on visible and near-infrared sensors are usually considered to be reasonably reliable under cloud-free conditions and are also available at high resolutions (e.g., MODIS provides SCF observations at spatial resolutions of 500 m, 0.05° and 0.25°); they however suffer from data gaps caused by cloud obscuration that can become quite extensive in wintertime over snow covered areas. Passive microwave sensors, on the other hand, can detect snow through clouds, but have difficulties detecting thin snow covers and may overestimate snow in dry and cold regions (e.g., [27]). As a result, multi-sensor assimilation strategies can be designed to combine the strengths of different sensors and mitigate the impact of limitations, so as to achieve the optimal utilization of satellite data from different sources (e.g., [15,16]). In the case of SCF observations, since they provide no direct information on SWE, SD, or snow density, the assimilation of this type of data is expected to be most effective for regions with ephemeral snow, and during initial accumulation and springtime melting for seasonal snowpacks (e.g., [13,14]). Combining the assimilation of SCF with the assimilation of SWE or SD may provide a larger benefit compared to assimilating SCF only (e.g., [15]). It is also worth noting that satellite-derived estimates of SWE or SD usually contain large uncertainties, due to errors in the retrieval algorithm, the instrument, or the spatial representativeness of the measurement. It is hence desirable to adjust the satellite-based estimates (e.g., as is for soil moisture in [28]) to reduce the bias before assimilating them, or to scale the satellite-based

209

estimates to the model climatology and only assimilate the anomalies (e.g., [15]). To mitigate the impact of uncertainty in satellite snow retrieval algorithms, radiance-based data assimilation has also been explored by coupling radiative transfer models to hydrologic or land surface models to allow the prediction of radiance observations from the satellites (e.g., [29–31]). Various techniques have been developed to assimilate earth science observations into land surface and hydrologic models (e.g., [32–34]), including the widely used ensemble Kalman Filter (EnKF, [35]). The EnKF technique, along with similar techniques such as the ensemble Kalman Smoother (EnKS) and the extended Kalman Filter (EKF), has been successfully employed to assimilate SWE or SD observations to update the snow states (e.g., [14,18,25,36]). However, given the linear updating scheme used by the EnKF, it is not considered to be most effective for highly non-linear processes such as rainfall-runoff simulations (e.g., [37]). This is also the case for SCF assimilation that relies on the nonlinear snow depletion curves (SDCs) for updating the snow states in the model (e.g., SWE and SD). The SDCs, which relate SCF to SWE or SD, can be subject to substantial nonlinearity and high spatial heterogeneity (e.g., [22]). In addition, seasonal snowpacks usually involve extended periods of 100% snow cover or snow-free conditions for which the model ensemble spread is too small for EnKF to take effect ([15]). Hence, rule-based direct insertion methods have also been developed to assimilate SCF data (e.g., [17,19,38]) or to supplement the EnKF assimilation of SCF (e.g., [15,17]). Other assimilation techniques such as particle filtering (PF), ensemble square root filtering (EnsRF) have also been explored for snow data assimilation and have demonstrated promises in improving snow prediction (e.g., [29,39]). For an in-depth review on the various assimilation methods typically used in hydrology including the EnKF technique, the readers are referred to Liu and Gupta ([32]). As discussed above, previous studies have suggested that satellite-based snow products of SCF and SD/SWE (from, e.g., MODIS and AMSR-E) can be useful in improving snow and snowmelt-driven runoff predictions; however, proper preprocessing procedures need to be taken in order to realize the greatest potential of these products in data assimilation applications and to avoid the negative impacts of errors in the retrieval algorithms or data gaps caused by cloud cover (e.g., [13,15,40–43]). The goal of the present study is to investigate the potential for improving large-scale snow prediction in a region in Alaska, via integrating high-resolution land surface modeling with assimilation of preprocessed SCF and SD estimates from satellite observations. Given the large interannual variability in snow accumulation, complex physiographic environment and sparse in situ observation network in Alaska, the prediction of snow-driven hydrology in this region is particularly challenging. In this study, we assimilate into the Noah land surface model (version 3.2; [44]) the SCF estimates from the standard MODIS product MOD10A1, and the interpreted SCF estimates and the bias-corrected SD estimates from the global snow product developed using the Air Force Weather Agency (AFWA)–National Aeronautics and Space Administration (NASA) Snow Algorithm (ANSA, [27]). Impacts of assimilating these SCF and SD products on the prediction of SWE, SD, and streamflow are analyzed. The study domain, the Noah land surface model and the datasets used for evaluation are presented in Section 2, followed by a discussion in Section 3 on the experimental design, the satellite-based products used for assimilation, the ANSA interpretation algorithm for generating SCF, the bias adjustment of the ANSA SD estimates, and the assimilation methodologies. The results from the DA experiments are evaluated against SNOTEL and COOP observations, the daily snow depth analysis from the Canadian Meteorological Center, and the USGS streamflow observations in Section 4. Conclusions, implications, and limitations of the present study, along with future work, are discussed in Section 5.

210

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

2. The Noah model, study domain, and evaluation datasets

[44], Chen and Dudhia [45], Koren et al. [46] and references cited therein.

2.1. The Noah land surface model 2.2. Study domain The Noah land surface model (version 3.2, hereinafter referred to as Noah; [44,45]) is used for predicting snow accumulation and ablation in the study region. The Noah model includes advanced physically-based parameterizations to simulate energy, momentum, and water exchanges between the land surface and the overlying atmosphere. It uses a blended snow–vegetation– surface soil layer and simulates the snow accumulation, sublimation, melting, and heat exchange processes. Although based on a single-layer physical parameterization, the Noah model does include some advanced treatments of snow processes by explicitly accounting for, for example, snow aging (with decreasing albedo), the retention and refreezing of liquid water (from snowmelt or rainfall) within the snowpack, frozen ground (which is important for the present study domain), as well as dynamic snow depth and density adjustments based on snow compaction overtime ([46]). The current version also improves upon the earlier versions of the Noah model (e.g., version 2.7.1) by considering more realistic representations of snow surface albedo and roughness following Livneh et al. [4]. The Noah model simulates all snow states such as SWE, SD, and snow density, and calculates fractional snow covered area (or SCF) via an empirical areal depletion curve. For more details on the snow physics and parameterization of other processes, the readers are referred to Ek et al.

The focus area in this study includes a 12° by 7° region in southcentral Alaska (153°–141°W, 59°–66°N; Fig. 1(a) and (b)), with elevations ranging from 0 to 6100 m. The study of snow hydrology in this region is particularly challenging given the extensive glaciation and permafrost coverage and the non-uniform climatic conditions: North of the Alaska Range (i.e., the interior) is dominated by subarctic continental climate with relatively warm summers and cold winters, while South of the Alaska Range is maritime with relatively cool summers and warm winters. All snow classes except for ephemeral snow as defined in Sturm et al. [47] can be found in this region. The situation is further complicated by the fact that Alaska has experienced substantial climate change during the last couple of decades with warmer temperatures and increased precipitation especially during the winter months (e.g., [48]). Although the snowmelt onset can vary from year to year by almost a month, snowmelt in such high latitude regions (primarily driven by radiation) is a very quick process involving only a short period of time (typically one week to 10 days; [49]). Snowmelt and/or glacier outburst are the dominant drivers of surface runoff in the region (although in the coast areas summer rainfall can account for a substantial portion of the runoff), with streamflow in most rivers typically peaking in the summer. Land cover in the study region is

Fig. 1. (a) The location of the study domain in the broader area of Alaska; (b) the study domain (59°–66°N, 153°–141°W) with elevations ranging from 0 to 6100 m; (c) highlight (in white) of areas with an elevation above 1500 m in the study domain. Also shown in (b) are the locations of the SNOTEL, COOP, and USGS streamflow gage stations.

211

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

dominated by glaciers along the mountain ranges (18%), woodland and forest (44%), shrub land (33%), and grassland (5%). 2.3. Evaluation datasets The study domain includes 27 SNOTEL and 90 COOP stations (solid blue and red circles in Fig. 1(b), respectively) that have at least 100 daily non-zero observations of SWE or SD during the study period (10/1/2002–9/30/2009) for evaluating the corresponding simulations from the Noah model. The elevations of the SNOTEL sites range from 46 to 1082 m with a mean of 490 m, while the COOP sites are located at lower elevations ranging from 18 to 922 m, with a mean of 271 m. Most of these stations are located in either forests/woodlands (40% for SNOTEL and 58% for COOP), shrub lands (56% for SNOTEL and 15% for COOP), and grasslands (4% for SNOTEL and 15% for COOP). To be able to evaluate the model simulations across the entire study domain, the spatially distributed SD estimates from the Canadian Meteorological Centre (CMC) daily snow depth analysis dataset ([50]) are also used. The CMC daily SD analysis is only available at a coarse resolution of 25 km; when used for evaluating the 1-km model simulations, the CMC data is remapped to the model grid via a nearest-neighbor approach, i.e., for any given model grid, the SD value from the nearest grid of the CMC dataset is assigned. As an additional independent source of data for the evaluation, monthly streamflow observations from a selection of nine US Geological Survey (USGS) stream gages (black triangles in Fig. 1(b); Table 1) are also used, with drainage areas ranging from 141 to 25600 mi2. These gages are chosen to ensure that there is a continuous record of monthly discharge observations during the study period. Since the Noah model does not include a routing program to compute the stream discharge, a simple approach is adopted here to calculate the monthly streamflow following Schwanghart and Kuhn [51]. First, the daily gridded runoff simulation from the Noah model and the surface area of the grid cells (determined based on the latitude and longitude of the grid cell center) are used to calculate the gridded daily flows. Then, based on the 1-km global digital elevation model (DEM) GTOPO30 (http://www1.gsi.go.jp/ geowww/globalmap-gsi/gtopo30/gtopo30.html), a set of Matlab functions for topographic analysis from the TopoTooxbox ([51]) is used to calculate flow direction and accumulation patterns in order to identify the streamflow network and the contributing grid cells for each gage. The drainage areas of the gages, obtained by totaling the surface areas of all contributing grid cells for each gage, are reasonably close to the drainage areas provided by the USGS (i.e., the relative difference is less than 10% for all gages; Table 1). Finally, the gridded daily flow and the streamflow network are used together to determine the mean monthly flows for each gage. On a monthly scale, such approximation is considered to be reasonable for small- to medium-sized basins that are free of regulations from man-made or natural structures such as reservoirs or

in-land lakes. However, for large basins with long flow concentration times (e.g., Gage #4 in Table 1) or basins subject to regulations, this simple approach may introduce timing (or phase) and magnitude errors in the resulting streamflow time series. 3. Experimental setup, satellite datasets, and assimilation strategies 3.1. Experimental setup To illustrate the impact of assimilating satellite-derived estimates of SCF and SD and the impact of cloud cover on assimilating SCF, a series of experiments are carried out to assimilate different sources of SCF and SD estimates, separately or jointly, into the Noah model. All simulation runs are carried out within the NASA Land Information System ([52,53]) at 1-km, hourly resolutions for a period of seven water years, from October 1, 2002 to September 30, 2009. Here a water year (WY) is defined as the period from October 1 of the previous calendar year to September 30 of the current calendar year. Forcing data for precipitation, air temperature, and other hydrometeorologic fields required to run the Noah model are obtained from the National Centers for Environmental Prediction (NCEP) Global Data Assimilation System (GDAS), which produces meteorological forcing analyses at a spatial resolution of 0.703°, 0.469°, or 0.313° (depending on the time the analyses are issued) every 6 h at four synoptic hours (0000, 0600, 1200, and 1800 UTC) for lead times 0-, 3-, 6-, and 9-hours. The original GDAS dataset is processed and sub-setted for driving off-line land surface models, as part of the effort for the NASA Global Land Data Assimilation System (GLDAS; [54]). LIS interpolates the GDAS forcing spatially and temporally at runtime to the required resolutions for driving the Noah model, using lapse rate and hypsometric corrections described by Cosgrove et al. [55]. To minimize the impact of initial conditions (e.g., soil moisture and snowpack SWE contents), the Noah model is run from arbitrary initial conditions for 10 years by cycling the GDAS forcing data from October 1, 2000 to September 30, 2002. 3.2. Satellite-based products used for assimilation 3.2.1. MOD10A1 One of the SCF datasets used for assimilation in this study is from the MOD10A1 product derived from observations by the MODIS optical sensor onboard the Terra satellite ([26]). The retrieval algorithm uses the at-satellite reflectances from MODIS band 4 (visible, 0.545–0.565 lm) and band 6 (near-infrared, 1.628–1.652 lm) to compute the normalized difference snow index (NDSI) for the determination of SCF ([56]). The MOD10A1 product is provided at daily and 500 m resolutions and includes (via Collection 5 preprocessing) improved treatments for snow/cloud discrimination and

Table 1 USGS stream gages in the study domain used for evaluation. Gage

USGS_ID

Description

Latitude

Longitude

A (mi2)

A1 (mi2)

DA (%)

1 2 3 4 5 6 7 8 9

15484000 15493000 15511000 15515500 15292000 15266300 15272380 15271000 15292700

SALCHA R NR SALCHAKET AK CHENA R NR TWO RIVERS AK CHENA R NR FAIRBANKS AK TANANA R AT NENANA AK SUSITNA R AT GOLD CREEK AK KENAI R AT SOLDOTNA AK TWENTYMILE R BL GLACIER R NR SIXMILE C NR HOPE AK TALKEETNA R NR TALKEETNA AK

64.472 64.902 64.886 64.565 62.767 60.477 60.897 60.82 62.346

146.926 146.359 147.25 149.094 149.693 151.082 148.924 149.427 150.019

2170 937 372 25600 6160 1951 141 234 1996

2202 961 370 24184 5522 2018 144 219 2001

1 3 -1 -6 -10 3 2 -6 0

Note: A is the drainage area from USGS, A1 is the surface contributing area computed based on the 1-km global digital elevation model (DEM) GTOPO30, and DA (%) represents the percent difference between A and A1 relative to A.

212

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

thin snow detection issues ([26]). These improvements, along with the high spatial resolution, should enhance the product’s usefulness for SCF assimilation applications. Nevertheless, it is important to note that the MODIS SCF product has inherent limitations for snow mapping in high-latitude areas (such as the present study region) due to extensive data gaps resulting from polar darkness during winter months (November–January) and the frequent, extensive cloud cover over the entire snow season. Given that the SCF product is most informative during initial accumulation and melting periods (which are both very short for areas in Alaska), the polar darkness and cloud cover issues could dramatically reduce the information content of this product for snow assimilation. To be consistent with the model resolution, the MODIS SCF dataset is re-projected from 500 m to 0.01° using the NASA’s MODIS Reprojection Tool (MRT) with the nearest-neighbor resampling option. At high latitudes of the study domain, the 0.01° resolution in LIS’s geographic projection is actually finer than MODIS’s 500 m resolution in its equal-area projection. 3.2.2. ANSA SCF estimates Previous studies have reported that removing the cloud gaps in MODIS SCF data can increase cloud-free observations and improve the snow mapping accuracy by combining the strengths of MODIS SCF and AMSR-E SWE products (e.g., [57]). Hall et al. [42] also showed that assimilating the cloud-gap-filled MODIS SCF product (instead of the original product) can improve the snow predictions from a land surface model. In the present study, to investigate the impact of cloud cover and snow mapping accuracy on SCF assimilation, the SCF estimates from the ANSA product ([27]) are also explored for assimilation. The ANSA SCF dataset was obtained by

applying a snow interpretation algorithm to the MOD10C1 product, which provides maps of snow cover and cloud cover at 0.05° in the climate modeling grid (CMG) projection based on binning the 500 m MOD10A1 product (see more details on this in the MODIS Snow Products User Guide [58]). The ANSA snow interpretation algorithm employs a set of rules to adjust the SCF values from MOD10C1, based on the fractions of snow, cloud, and snow free land in a grid cell. For example, for any cell with a mix of snow and cloud, it is marked as cloud covered if the fraction of cloud cover is greater than 80%; otherwise, if the fraction of snow is greater than 20%, continuous snow cover is assumed for the cell and the SCF value is changed to 100. For any cell with a mix of snow, cloud, and snow free land, the SCF value is changed to zero if the snow fraction is less than 15%; otherwise, it is increased by some amount based on the cloud cover. The details on the snow interpretation algorithm are provided in Fig. 2. These adjustments are expected to reduce the cloud cover-induced data gaps in the SCF dataset, while at the same time improving snow mapping accuracy via alleviating snow and cloud errors of commission over problematic areas such as cloud shadowed vegetated surfaces, mixes of vegetation and water, coastlines or periphery of clouds. The ANSA product also includes a binary snow cover dataset with data gaps caused by cloud cover and polar darkness removed by blending information from the MODIS and AMSR-E instruments. In the current study, this binary dataset is not used; rather, the interpreted fractional dataset (ANSA SCF) is assimilated for comparison with the assimilation of the standard fractional dataset in MOD10A1. Fig. 3 shows the SCF distributions (along with cloud cover) for the MODIS product and the ANSA product for a typical day during the accumulation, peak snow, and melting seasons (12/15/2008, 3/

Fig. 2. Flowchart of the snow interpretation algorithm for generating the ANSA SCF dataset, where FSNOW, FCLOUD, and FSNOWFREE are the fractions of snow, cloud, and snow free land of a grid cell, respectively.

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

213

Fig. 3. Comparison of SCF from MODIS and ANSA for selected days during the accumulation, peak snow, and melting seasons (12/15/2008, 3/15/2009, 5/15/2009). Light grey = cloud cover, dark grey = missing data (e.g., due to polar darkness), and white = ocean.

15/2009, and 5/15/2009) during WY 2009. Note that for 12/15/ 2008, a large portion of the domain (dark grey) is covered by missing data due to polar darkness. It can be observed that the cloud cover (light grey) for these three days has been reduced to some extent in the ANSA product (as compared to the original MODIS product). The percent reduction in cloud coverage for ANSA with respect to MODIS is 14%, 32%, and 7% for 12/15/2008, 3/15/2009, and 5/15/2009, respectively. Over the entire analysis period (10/ 1/2002–9/30/2009), the average percent reduction in cloud cover for December, March, and May is 17%, 29%, and 27%, respectively. These reductions are remarkable given the extensiveness of cloud cover in this region. In addition, it can be observed that on 3/15/ 2009 (a typical day during the peak snow season), the snow distribution is more realistically represented in the interpreted ANSA SCF than in the standard MODIS SCF, with more spatially continuous (100%) coverage of snow as one would expect for the Alaskan study domain during the peak snow season. The average area of 100% snow cover in March over the 7-year analysis period is 40% for ANSA SCF, as compared 14% for the MODIS SCF. It is noted that the ANSA product is available at only a coarser resolution of 5-km, owing to the coarse resolution of the MOD10C1 product. It would be worthwhile to investigate how the gain from cloud cover reduction and improved snow mapping accuracy and the loss in spatial resolution interact with each other during the assimilation process. 3.2.3. ANSA SD estimates and bias correction Also used for assimilation in this study are the SD estimates from the ANSA product. Unlike the old AMSR-E algorithm ([59]) that derives SWE based on the differences in brightness temperature (Tb) from only the 18.7 and 36.5 GHz channels, the current algorithm in ANSA ([60]) uses Tb observations from several more channels (10, 18, 22, 89 GHz vertical and horizontal polarized signals) on the AMSR-E instrument to derive SD estimates first and then convert SD to SWE based on a snow density map for the seasonal snow classes delineated by Sturm et al. [47]. In addition, MODIS cloud-free observations are used for quality checking the microwave-derived SD/SWE values. In comparison to the previous SWE algorithm based on Chang et al. ([59]), the ANSA SD algorithm

is expected to provide increased dynamic range for deep snow, and improved capabilities for mapping shallow and forested snow covers. More details of the ANSA SD and SWE retrieval algorithms are provided in Foster et al. ([27]). Prior to the assimilation, the 5-km ANSA SD dataset is adjusted against the in situ SD observations from the COOP stations and SWE observations from the SNOTEL stations to reduce the potentially large bias in the SD dataset. In order to use the SNOTEL observations, the SWE values are converted to SD values based on the snow density values presented in [61] for each snow class. It is worth noting such assumptions of constant snow density for each snow class may introduce a certain level of uncertainty into the analysis; however, since snow depth typically varies over a much larger range than snow density ([61]), here we assume that the additional uncertainty introduced by using constant snow density has relatively limited impact on the bias adjustment. The statistical interpolation (or optimum interpolation) method presented in [21] for developing the CMC daily SD analysis is adopted here to blend the information from the ANSA SD dataset (the background) and the SNOTEL/COOP observations. The ANSA SD estimates are adjusted grid point by grid point with corrections computed from weighted sums of differences between the station observations and the background estimates at the station locations:

SbJ ¼ Sj þ

N X

wi;j ðOi  Si Þ

ð1Þ

i¼1

where Sj and Sbj denote the ANSA SD estimate and the adjusted SD estimate at grid cell j, respectively; Oi is the SD observation at station i; Si is the ANSA SD estimate for the grid cell at which station i is located; N is the number of stations used for the bias adjustment (here N = 117, including 27 SNOTEL stations and 90 COOP stations available in the study domain); and wi,j is the optimum weight associated with station i for calculating the adjustment to the SD value at grid cell j. The calculation of the weights takes into account both the horizontal and vertical separations between the background and the station locations, with an e-folding distance of 120 km in the horizontal direction and 800 m in the vertical direction (e.g.,

214

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

the elevation difference). The influence of errors in the station observations and the background field is also represented in the calculation of the weights. In the present study, the standard deviation of the observation error is chosen to be 24.5 mm (same as in [21]), while the standard deviation of the error in the ANSA SD estimates is arbitrarily specified to be 100 mm to reflect the large bias in the ANSA dataset. The specification of the error variances of the station and satellite data will have an influence on the weight factor (i.e., the bigger the error variance for satellite data relative to that of station data, the bigger the weight factor) but is noted to be not as important as in data assimilation applications ([21]). The readers are referred to Brasnet [21] for details on the calculation of the weights. In order to investigate the influence of using observations from the SNOTEL sites that are typically located along the snowlines at high elevations, two bias correction experiments are conducted: (1) including both SNOTEL and COOP observations in the bias adjustment and (2) including only COOP observations in the bias adjustment. The SD datasets produced from these two experiments are denoted here as ANSA.adj and ANSA.adj.coop, respectively. In Fig. 4, the mean annual peak SWE and SD estimates from ANSA (before and after bias correction) during peak snow seasons are compared with point observations from the SNOTEL and COOP stations, respectively. The least-squares trend lines are also plotted. Given that snow in most areas in the study region typically starts to melt in late April or early May, the peak snow season is defined here to be March 15–April 15. Fig. 4 shows that the original ANSA product considerably underestimates the mean annual peak SWE at all SNOTEL sites and the mean annual SD at most COOP sites. After bias adjustment, the underestimation is significantly reduced, regardless of using the SNOTEL observations in the

adjustment or not. Using both SNOTEL and COOP observations (i.e., ANSA.adj) appears to have achieved the best improvement (especially for the SNOTEL sites), with the correlation coefficient (R) between the station observations and the ANSA estimates improved from 0.42 to 0.94 for the SNOTEL sites and from 0.05 to 0.38 for the COOP sites (see also a comparison of the R values in Table 2). Also shown in Fig. 4 are the SWE and SD estimates from the open loop simulation of the Noah model (i.e., the simulation with no assimilation), which overestimates SWE for the deep snowpacks but underestimates it for the shallower snowpacks at the SNOTEL sites. At the COOP sites, the Noah model underestimates SD at most stations but overestimates at the rest, with an R value of 0.88 and 0.24 for the SNOTEL and COOP sites, respectively. It is important to note that, in addition to deficiencies in the satellite retrieval algorithm (or the model’s physical parameterization), a portion of the difference in the snow estimates between SNOTEL/COOP and ANSA (or model) may be attributed to uncertainty introduced by scale and representativeness issues. Given the overall superior performance resulting from bias adjustment using both SNOTEL and COOP observations, the discussion from now on will focus primarily on the ANSA.adj dataset. To examine the influence of bias adjustment on the ANSA SD estimates across the entire study domain, the maps of difference in mean SD between ANSA and ANSA.adj for December, March, and May over the 9-year analysis period are presented in Fig. 5. Although widely ranging from 1000 to +1000 mm, the adjustment is mostly positive (i.e., adding snow to the ANSA estimate) across the entire study domain, with larger corrections towards the southern part of the domain which is dominated by maritime climate and weather patterns and receives a much higher amount of winter precipitation than the northern part of the domain. The

Fig. 4. Comparisons of mean annual peak SWE (a) and SD (b) of in situ observations (from SNOTEL or COOP stations) vs. ANSA estimates before and after bias adjustment, and the open loop simulations from the Noah model. ANSA.adj indicates both SNOTEL and COOP observations are used in the bias adjustment, while ANSA.adj.coop indicates only COOP observations are used in the bias adjustment.

Table 2 Correlation coefficient of mean annual peak SWE and SD estimates at the SNOTEL and COOP stations.

SNOTEL COOP

ANSA

ANSA.adj

ANSA.adj.coop

OL

MODIS_DA

SCF_DA

SD_DA

SCF + SD_DA

0.42 0.05

0.94 0.38

0.74 0.33

0.88 0.24

0.88 0.25

0.9 0.29

0.86 0.63

0.86 0.63

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

215

Fig. 5. Maps of difference in mean SD between ANSA and ANSA.adj (ANAS.adj – ANSA) for December, March, and May (averaged across the analysis period WY 2003–2009). Note the data is categorized into 8 groups for better visualization.

Fig. 6. The mean SD maps for December, March, and May for the original 5-km ANSA product (ANSA), the new 5-km, bias corrected ANSA product (ANSA.adj), the 25-km CMC daily snow depth analysis, and the 1-km open loop simulations of the Noah model. Note the data is categorized for better visualization.

larger negative bias in the southern domain can be partially attributed to the increasing uncertainty in the passive microwave-based snow retrieval for deep snow packs due to various detection issues including saturation (see more discussion in the next paragraph). It can also be observed that the corrections are larger and more widespread for December and March than for May, when snow in low elevations has melted out. The bias adjustment process has also removed snow from the original ANSA SD estimates in some small areas, as indicated by the dark and light blue areas. The negative adjustment (i.e., removing snow) near the middle of the eastern edge of the domain is likely due to the sparse in situ station network there that includes several COOP stations and only one

SNOTEL station (Fig. 1). As a result, the bias-adjusted SD estimates may be less reliable in this region. Fig. 6 compares the mean SD maps for December, March, and May over the entire study period (10/2002–9/2009) for the original ANSA product (ANSA), the new, bias adjusted ANSA product (ANSA.adj), the CMC daily SD analysis, and the open loop simulation of the Noah model. When compared against the CMC analysis, the underestimation by the original ANSA product over the mountain ranges is quite substantial (as indicated by the isolated dark blue polygons with green–red edges). This underestimation is partly attributed to saturation – the AMSR-E retrieval algorithm is insensitive to snow depths greater than about 0.75 m

216

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

(equivalent to SWE values greater than about 300 mm). In addition, forests, especially coniferous forests, confound passive microwave algorithms since vegetation is emissive in passive microwave frequencies and can negate or minimize the scattering exhibited by an underlying snow pack. Moreover, mountainous terrain is known to affect microwave scattering. Snowpack scattering signals from various slopes, particularly steep slopes, are reflected differently and thus brightness temperature (Tb) values recorded at the satellite sensor are not always representative of the observed thickness of the snow pack. Glaciers in these areas may have also contributed to the complication. The bias adjusted ANSA product agrees better with CMC than the original ANSA product by providing more reasonable SD estimates for deep snowpacks along the major mountain ranges, especially during the accumulation season (December and March). The bias adjustment, however, tends to estimate more snow than CMC in the inland areas north of the Alaska Range during melting time (i.e., May). The Noah model is able to produce spatial snowpack patterns most consistent with that of the CMC analysis, partially attributable to the fact that the latter is derived from a modeling-based background field ([21]). Nevertheless, the Noah model tends to estimate more snow than the CMC analysis in the deep mountain snowpacks consistently through all seasons. This is consistent with earlier observations from Fig. 4 where the model tends to overestimate deep snowpacks at the SNOTEL sites. It is important to note here that the CMC analysis is considered to be a reliable source for evaluation but not necessarily the truth. Since the bias adjustment relies on the availability of in situ observations that may not be adequately representative in the horizontal and vertical dimensions, the adjusted ANSA SD estimates should not be expected to be free of bias everywhere across the study domain. Given that the elevation of the stations ranges from 18 to 1082 m, there is less confidence in the accuracy of the adjusted SD estimates at very high elevations (e.g., >1500 m). To avoid the potential negative impact of biased ANSA SD estimates at very high elevations, it is desirable to limit the assimilation of the adjusted ANSA SD estimates to areas with an elevation lower than a certain threshold. In the present study, this threshold of elevation is chosen to be 1500 m. The areas with an elevation higher than 1500 m account for about 10% of the entire domain and are highlighted in white in Fig. 1(c). For the purpose of a comprehensive investigation, assimilation experiments are carried out using SD estimates from both ANSA.adj and ANSA.adj.coop, and for both the entire study domain and only those areas with an elevation lower than 1500 m. The results (not shown) have confirmed that including the areas with a very high elevation in the assimilation would result in significant underestimation of SD in these areas and less improvement in the streamflow predictions for gages that receive snowmelt water from these high-elevation areas. For the convenience of presentation, the focus of discussion hereinafter for SD assimilation will be on assimilating the bias-adjusted ANSA SD estimates from using both SNOTEL and COOP observations (i.e., the ANSA.adj dataset) for areas with an elevation lower than 1500 m only.

3.3. Assimilation strategies Both SD and SCF estimates are assimilated at a daily frequency, regardless of whether they are assimilated separately or jointly. The SCF estimates are assimilated at 10:00 local solar time and the ANSA SD estimates are assimilated at 14:00 local solar time, to be close with the overpass times of MODIS and AMSR-E, respectively. Relevant details on the algorithms for assimilating SD and SCF estimates are given below. 3.3.1. SD assimilation In the present study, the standard EnKF technique ([35]) is used to assimilate the bias-adjusted ANSA SD estimates (the term ‘‘biasadjusted’’ is omitted hereinafter for convenience). Since high-resolution land DA runs are expensive and it is not practical to conduct a comprehensive sensitivity analysis to determine a preferable size for the ensembles, 20 members are used for the EnKF runs following De Lannoy et al. ([15]). During the assimilation, perturbations are introduced into the meteorological forcing fields, the snow states of the model (SWE and SD), and the ANSA SD estimates to represent the uncertainty in these fields (Table 3). Following recommendations in earlier studies (e.g., [15,62]), multiplicative perturbations with a mean value of one are imposed on shortwave radiation (SW) and precipitation (P), while additive perturbations with a mean value of zero are applied to near-surface air temperature (T) and long-wave radiation (LW). The Noah model states SWE and SD are also assumed to be subject to multiplicative noises; so are the SD estimates from the ANSA product used as the ‘‘observation’’ for assimilation. Cross correlation is imposed on the perturbations of forcing variables (P, T, SW, and LW) and state variables (SWE and SD). In addition, time series correlations are imposed via a first-order regressive model with a time scale of three hours for all variables. This may not be optimal since in reality some variables (e.g., SWE) may have a longer temporal correlation scale than 3 h. A one-dimensional version of the EnKF technique with no spatial correlation is used for the assimilation. It is noted that an ensemble open loop (OL) simulation run with similarly perturbed forcing and state variables is also carried out as the benchmark for evaluation. When evaluated against the observations from the SNOTEL/COOP stations, the CMC analysis and USGS streamflow measurements, the simulations based on the ensemble mean are used. 3.3.2. SCF assimilation For the assimilation of the SCF data (from MODIS or ANSA), a rule-based direct insertion strategy is developed. Such a strategy is considered to be more effective than an EnKF-based algorithm for continuous, deep seasonal snowpacks, since the latter does not work when the model simulates zero snow or 100% snow cover (e.g., [15]), which can last for an extended period of time (11 months) for any given year in the current study domain. Considering that the SCF data is most informative during the melting of seasonal snowpack that typically lasts for only 5–10 days in this area (e.g., [49]), and that the Noah model generally overestimates the deep snowpacks (see discussions in Section 3.1), the original

Table 3 Parameters used for the perturbation in the EnKF-based assimilation of the ANSA SD estimates. Variables

Perturbation method

Standard deviation

Cross correlation between variables

Forcing

SW LW P T

Multiplicative Additive Multiplicative Additive

0.1 (–) 15 (W/m2) 0.5 (–) 0.5 (K)

1.0 0.3 0.5 0.3

0.3 1.0 0.5 0.6

0.5 0.5 1.0 0.1

0.3 0.6 0.1 1.0

Model states

SWE SD

Multiplicative Multiplicative

0.01 (–) 0.01 (–)

1.0 0.9

0.9 1.0

– –

– –

Observation

ANSA SD

Multiplicative

0.2 (–)

–

–

–

–

217

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

technique presented in Rodell and Houser ([38]) is adjusted to enable the effective utilization of the information contained in the available SCF observations. Specifically, if the observation indicates very little snow on the ground at a specific date during the melt season while the model simulates a considerable amount of snow (e.g., larger than 20 mm of SWE), an arbitrary melt rate is assumed and a corresponding amount of SWE is removed from the snowpack for that specific day. Here the melt rate at any given day is assumed to be the constant rate it would take for the current snowpack to melt out completely in 10 days (with an upper boundary for 50 mm of SWE per day). Note that the melt rate is dynamic (usually decreasing) rather than constant, as the SWE amount is changing every day (usually decreasing during the melt season). Together with the maximum melt rate of 50 mm/day, the use of a dynamic melt rate will prevent unrealistically melting out all the snow/ice at locations where glaciers exist or where the snowpack is too deep to melt out during a single year. When there is less than 20 mm of SWE left, the snowpack is assumed to melt out completely within one day if the observation indicates no snow. On the other hand, if the model simulates less than 10 mm of SWE while the observation indicates a fractional snow cover larger than 40%, a nominal amount of SWE (20 mm) is added to the snowpack. A dynamic snow density (computed based on snowpack compaction algorithm in the Noah model) is used to convert SWE into SD for evaluation. A flow chart of the SCF assimilation strategy is provided in Fig. 7.

4. Evaluation of the results In this section, the results from four representative assimilation cases will be presented and discussed. These include separately assimilating MODIS SCF, ANSA SCF, and ANSA SD products, as well as the joint assimilation of ANSA SCF and SD products. For the convenience of presentation, these four assimilation experiments will be referred to as MODIS_DA, SCF_DA, SD_DA, and SCF + SD_DA, respectively. The results from the OL and DA simulations are compared against several sources of observations (i.e., SNOTEL/COOP, CMC analysis, and streamflow) to investigate the potential gain or loss in the skill of snow and streamflow predictions from assimilating satellite-based SCF and/or SD estimates. Several metrics are used for the evaluation, including Root Mean Square Error (RMSE), bias (or mean error), and Nash Sutcliffe Efficiency (NSE). These are defined below:

0PN RMSE ¼ @

i¼1 ðM i

i¼1 ðM i

A

¼ Oi Þ

2 i¼1 ðM i  Oi Þ PN 2 i¼1 ðO  Oi Þ

ð2Þ

ð3Þ

N PN

NSE ¼ 1 

10:5

N

PN Bias ¼

 Oi Þ2

! ð4Þ

Fig. 7. Flowchart of the rule-based direct insertion strategy for SCF assimilation for a period of N days. SDEN = dynamic snow density; MR = the constant melt rate (mm/day) at which the current snowpack would melt out completely in 10 days.

218

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

where M and O denote model simulations and observations, respectively. O is the mean observation over the entire simulation period. N is the total number of observation-model simulation pairs used for calculating the statistic. RMSE and Bias are used for evaluating SD or SWE predictions at the daily scale, while NSE and RMSE are used for evaluating streamflow prediction at the monthly scale. 4.1. Evaluation against in situ SWE/SD observations First, the SWE/SD simulations at individual SNOTEL/COOP sites from the different assimilation runs are examined to evaluate the skill improvement with assimilation (if any) relative to the open loop run. The results indicate that different snow assimilation scenarios may have remarkably different impacts on the SWE or SD simulations at a given location, depending on the physiographic situation of the location (e.g., elevation and land cover) and how the model and assimilation dataset perform at the location. Fig. 8 presents a comparison of SWE or SD time series from the station observations, the open loop simulation, and the simulations from the four different snow assimilation cases for three representative locations, including two SNOTEL stations (49M08S and 51K14S) and one COOP station (502107) that are located at different elevations and environments (see Fig. 1(b) and Table 4 for the locations and descriptions of the three stations). The time periods shown for

these three locations are WY 2004, 2009, and 2005, respectively. For both SNOTEL sties, the model OL simulation tends to overestimate the snow accumulation and also hold snow for too long during the May–June time frame, while the station observations indicate that snow has melt out completely a few weeks earlier. For 49M08S, the assimilation of ANSA SCF estimates (i.e., SCF_DA) via the rule-based strategy described in Section 3.2.2 effectively improves the snow accumulation and ablation processes so that the simulated SWE time series is in close agreement with the observation. The assimilation of MODIS SCF (MODIS_DA) using the same strategy, however, has little impact on the SWE simulation, due to the presence of cloud cover that has reduced the availability of observations for assimilation especially during melting (not shown). At 51K14S, both MODIS_DA and SCF_DA are able to effectively produce a snow depletion curve that is in better agreement with the observation, with the latter still having a larger impact (or better performance) owing to the reduced cloud coverage and snow mapping accuracy. Designed to work most effectively when the model is overestimating, no benefit from the SCF assimilations is observed at the COOP site 502107 where the OL simulation underestimates the observation. The EnKF-based assimilation of the ANSA SD estimates (SD_DA), on the other hand, considerably improves the SD simulation at 502107 with effectively reduced underestimation.

Fig. 8. Comparison of SWE or SD time series of the station observations and the simulations from OL and DA runs for two SNOTEL stations (49M08S and 51K14S) and one COOP station (502107) for the water year of 2004, 2009, and 2005, respectively.

Table 4 Site description and SWE/SD simulation statistics for three selected sites.

Latitude Longitude Elevation (m) Land cover RMSE (mm)

Bias (mm)

OL MODIS_DA SCF_DA SD_DA SCF + SD_DA OL MODIS_DA SCF_DA SD_DA SCF+SD_DA

49M08S (SWE)

51K14S (SWE)

502107 (SD)

61.067 149.467 716 Closed shrubland 161 159 (2%) 69 (57%) 168 (4%) 168 (4%) 117 114 (2%) 42 (64%) 116 (0%) 118 (1%)

59.733 151.250 402 Grassland 249 198 (21%) 177 (29%) 68 (73%) 69 (72%) 98 62 (37%) 46 (53%) 37 (62%) 39 (61%)

64.850 147.833 182 Evergreen Needleleaf Forest 120 115 (4%) 106 (12%) 42 (65%) 41 (66%) 71 68 (4%) 50 (30%) 16 (78%) 16 (78%)

Note: Bold values indicate the best performance; values in the parentheses refer to the percent reductions in RMSE or absolute bias relative to the open loop.

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

The SD assimilation tends to introduce underestimation at both SNOTEL sites (mostly during accumulation), likely due to AMSRE’s limitation in measuring deep, wet maritime snowpacks at these high elevation sites. Since only assimilating SD has already captured well the snow depletion during the melting time at these two locations, adding SCF to the assimilation (which helps most during melting) has not led to noticeable improvement on the SWE assimilation. These observations are also supported by the overall statistics (RMSE and bias) calculated over the snow season (October–June) of the entire study period (Table 4). Table 4 shows that for the three stations 49M08S, 51K14S, and 502107, the RMSE is reduced by as much as 57%, 73%, and 66%, respectively; and the bias is reduced by 64%, 62%, and 78%, respectively. The best performances are achieved by assimilating ANSA SCF for 49M08S and by assimilating ANSA SD or jointly assimilating ANSA SCF and SD estimates for the other two sites. When considering the mean annual peak SWE or SD simulated for all the individual SNOTEL or COOP sites (Fig. 9), the assimilation of SCF information has little impact on snow simulation with MODIS_DA, while some moderate improvements are achieved with SCF_DA based on the ANSA SCF estimates for both SNOTEL and COOP sites (see also the R values in Table 2). The SD-based assimilation experiments (SD_DA and SCF + SD_DA) have resulted in considerable improvements in the mean annual peak SD estimation at the COOP sites, where the correlation with the observations are improved from 0.24 for the OL run to 0.63 for both DA runs. At the SNOTEL sites, however, the SD-based assimilations shows little impact on the estimation of the mean annual peak SWE (in terms of correlation coefficient). To examine the average seasonal cycle of the DA performance, the mean statistics across all SNOTEL/COOP stations for each individual month are presented in Fig. 10. Since the SWE or SD values can vary dramatically from one station to another, the mean RMSE and bias of the SWE or SD anomalies (instead of those of the original SWE or SD quantities) are presented. This will make averaging across stations more meaningful, while at the same time reducing

219

issues of scale mismatch between model simulations and station observations. The anomaly time series are derived as deviations from the climatological seasonal cycle computed by smoothing the original model simulation or station observation time series with a 30-day window and then averaging over the seven water years of the study period. For the calculation of mean biases, the absolute values of the bias at individual stations are used to avoid cancellation of biases from different stations. Fig. 10 shows that there exists a strong seasonal pattern in the evolution of errors in snow estimation for both SNOTEL and COOP sties, as depicted here by RMSE and bias. The errors typically start small in October, gradually increase with each month to reach the maximum in the spring melt time, and then decrease and disappear by summertime as snow melts out. The error in the SWE estimates for the SNOTEL sites typically peaks in April; while at lower elevations of the COOP sites, the error in SD estimates tends to peak one month earlier in March. Overall, little or only marginal improvement over the OL run is achieved by MODIS_DA. With reduced cloud coverage and improved snow mapping accuracy, SCF_DA produces a larger improvement than MODIS_DA. The assimilation of SD estimates (SD_DA and SCF + SD_DA) has resulted in the largest reduction in both RMSE and bias for all months. Using both SD and SCF information in the assimilation does not seem to improve upon assimilating SD only, indicating again the dominant influence of SD in the joint assimilation experiment. 4.2. Evaluation against the CMC daily SD analysis Since the SNOTEL and COOP sites are mostly located at elevations below 1000 m (with only two SNOTEL sites located above 1000 m) and are not uniformly distributed over the study domain, the observations from these stations are not the best sources for validating the spatial distribution of the snowpacks. In this paper, the 25-km CMC daily SD analysis is used to assess how well the different snow DA experiments improves upon the OL simulation in predicting snow distribution over the entire study domain.

Fig. 9. Comparison of mean annual peak SWE (top row) and SD (bottom row) simulated by the OL (red) and DA (blue) runs against the observations from the individual SNOTEL or COOP stations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

220

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

Fig. 10. (a) The mean RMSE of SWE anomaly (mm) for the OL and DA simulations for each individual month (averaged across all SNOTEL sites); (b) same as (a) but for SD anomaly at the COOP sites; (c) same as (a) but for the mean absolute bias; (d) same as (b) but for the mean absolute bias.

In Fig. 11, the differences in RMSE distribution between the OL run and the different assimilation runs are shown. Again, when compared with the OL run, MODIS_DA seems to have very little impact on the SD prediction for all the three months (December, March, and May), except for a small area near the southwestern coast. For SCF_DA which uses the interpreted ANSA SCF, the impact is considerably increased across the entire domain and is mostly manifested as reductions in RMSE over the mountains, especially for the melt season (May). When SD estimates are used in the assimilation (SD_DA and SCF + SD_DA), the RMSE reduction is significant over most areas in the study domain (especially those along the mountain ranges). However, the assimilation of SD information also results in increased RMSE errors in some other regions, which appear to be more extensive but with smaller amounts during melting season (May) than the accumulation season (December and March). These observations are consistent with the spatial distribution of bias similarly calculated against CMC (not shown). Such inconsistent performance in SD prediction across the domain (i.e., improvement in some areas and degradation in other areas) from SD assimilation, as opposed to the more consistent performance of SCF_DA, may make it less straightforward and location dependent for the model to translate the skill in snow prediction into that of streamflow prediction, as will be discussed later in Section 4.3. Fig. 12 presents the seasonal cycle of domain-averaged RMSE and bias of the SD estimates from the OL and DA runs. For areas with an elevation below 500 m (accounting for 37% of the domain),

a similar seasonal pattern (as shown in Fig. 10 for evaluation against station observations) is observed for RMSE that increases as snow accumulates during winter time and decreases as snow melts out in spring and summer times, while large RMSE errors persist throughout the year for elevations above 500 m (accounting for 52% of the domain) where deep snowpacks or glaciers exist. Impact of the assimilation is largest for high elevation (>500 m) with significant reductions in the seasonal RMSE values especially when assimilating SD information. Similar conclusions can be reached when comparing the domain-averaged bias values for the different DA runs. The SD assimilation runs have drastically reduced the overestimation at high elevations (>500 m). For elevations lower than 500 m, all assimilation runs appear to exhibit limited effectiveness in reducing the bias, with SCF assimilations tending to introduce negative bias and SD assimilations positive bias during melting time (February–April). The positive bias introduced by assimilating SD at low elevations is likely due to using SNOTEL SWE observations in the bias adjustment process, which are typically located above the snowline where snowpack usually does not start melting until May. Overall, the mean RMSE and bias statistics presented in Fig. 12 indicate that the improvement in SD estimation by assimilating satellite-based SCF or SD data is more significant for high elevations than for low elevations, with SD assimilation having a larger impact than SCF assimilation. Joint assimilation of SCF and SD results in only marginal additional improvements upon assimilating just SD.

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

221

Fig. 11. Difference in RMSE in SD (against CMC) between the OL and DA runs for December, March, and May. Note the data is categorized for better visualization.

4.3. Evaluation of streamflow prediction As an integrated indicator of hydrologic fluxes within a basin, streamflow observations represent an additional independent source of reliable information for evaluating the impact of the snow assimilation experiments. Fig. 13 compares the mean monthly streamflow observations from the USGS with those computed based on the model OL and two DA runs (SCF_DA, SCF + SD_DA) for the nine selected USGS gage locations (Table 1). The streamflow simulations from MODIS_DA are very similar to those from the OL run and hence are not shown. The results from SD_DA are very similar to those from the SCF + SD_DA run and hence are not shown either. Surprisingly, despite that the SCF + SD_DA run outperforms SCF_DA in snow prediction as discussed in 4.2, the latter has resulted in larger improvements in streamflow prediction for all gages, suggesting that the improvement on snow prediction does not necessarily translate into improvement on streamflow prediction. At Gages 6, 8, and 9, the overestimation of monthly flow by the model OL is effectively reduced by assimilating ANSA SCF and/or SD estimates, suggesting the effectiveness of the SD bias adjustment and the ANSA SCF interpretation algorithm in these areas. However, at Gage 1 (Salcha River near Salchaket), Gage 2 (Chena River near Two Rivers), and Gage 3 (Chena River near Fairbanks, AK), SD_DA considerably overestimates the first flow peak in May. This may have caused, partially, by the use of observations from SNOTEL stations (which are typically located at high elevations with deep snowpacks) in the bias adjustment of the ANSA SD estimates, as indicated by the extensive positive adjustment in May in the northeast corner of the domain (where these three gages are located; Fig. 6). In addition, an assimilation run with the SD estimates from bias adjustment using only COOP observations (not shown) effectively removes the overestimation at these gages and exhibits improvement over the OL run; it however tends to introduce degradation in streamflow prediction at some other high-elevation gages (e.g., Gage 7 and Gage 8), where the use of observations from high-elevation SNOTEL stations for bias adjusting ANSA SD is critical. More discussion on this is provided later in this section. It is also shown in Fig. 13 that the model OL and assimilation runs all underestimate the flow at Gage 4 (Tenana River at Nenana,

AK) and especially Gage 7 (Twentymile River below Glacier River near Portage, AK) and peak too early in May, missing the larger observed flow peak in July. This is largely due to the summer-time glacier melt in the Alaska Range and the Chugach Mountains that feed into Gages 4 and 7 (respectively) and that the glacier driven runoff may not be handled properly in the Noah model. Although both SCF and SD assimilations have reduced the overestimation of snow in the glaciated areas, the melt water from glaciers (as a result of DA) are removed from the system without contributing to the runoff, which probably has also contributed to the underestimation at Gages 4 and 7. As mentioned earlier in Section 2, the monthly discharge is approximated by summing up runoff in all the upstream grid cells of a given gage. Such an approximation ignores the flow concentration time within a basin (which is a function of basin size, topography, geology and land use within the basin, among other factors); this may have caused the simulated monthly flow to peak earlier than the observed flow for a few gages including Gages 4 and 5, which have a large drainage area of 256,000 and 6160 square miles, respectively. Table 5 lists the RMSE, percent reduction in RMSE, NSE, and improvement in NSE for streamflow prediction at the nine gages computed on a monthly scale. It can be seen that SCF_DA outperforms all other assimilation cases and has resulted in an improvement in streamflow prediction at all gages, except for Gage 8 (Sixmile Creek near Hope) where the SD_DA and SCF + SD_DA have performed equally well as SCF_DA. For Gage 8, the RMSE of monthly streamflow is reduced by 33% and the NSE is improved by 0.87 (from 0.57 to 0.3). MODIS-based SCF assimilation has little impact on streamflow prediction, while the two SD assimilation experiments (SD_DA and SCF + SD_DA) exhibit every similar performance with improvements at only a few gages (6, 8, and 9) and degradations at the others. The encouraging performance of SCF_DA in streamflow prediction can be partially attributed to the its relatively consistent improvement in snow prediction across the entire study domain, as opposed to the inconsistent performance from assimilating ANSA SD (see Fig. 11). In addition, SCF_DA also benefits from the relatively conservative assimilation strategy that indirectly assimilates the SCF information to update the model states (SWE and SD) based on customized rules. Such a rule-based assimilation strategy is less sensitive to errors in the

222

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

Fig. 12. The domain-averaged RMSE of SD for (a) low-elevation areas (elevation 6 500 m) and (b) high-elevation areas (500 < elevation 6 1500 m); (c) and (d) are the same as a) and b), respectively, but for the domain-averaged bias of SD.

observations than the EnKF-based algorithm used for SD assimilation, where the model states are directly adjusted based on the actual values of ANSA SD estimates (after imposing random noise). To further decipher the insufficient performance of SD assimilation in predicting streamflow at most gages, the mean RMSE and bias (against the CMC analysis) averaged across the contributing area of each gage are examined for the months of March, April, and May in Figs. 14 and 15, respectively. Indeed, the SD assimilation runs have resulted in higher RMSE values and larger positive biases (than the OL and SCF_DA runs) on SD estimation across the contributing areas of Gages 1–3 for April and May. However, for March the RMSE and bias values from SD assimilation represent an improvement over the OL and SCF assimilation runs. An analysis of the elevations indicates that these three gages are located at relatively low elevations (especially for Gage 3, which has more than 60% of its contributing grid cells located at an elevation lower than 500 m and 96% lower than 800 m). At these low elevation locations, snow would have largely melt out by April, whereas at higher elevations (e.g., where the SNOTEL sites are located), the snowpacks would continue to accumulate throughout April and would not start melting until May. Consequently, using SNOTEL observations in the bias adjustment could introduce positive bias

in April and May in snow prediction in these lower elevation basins, which would subsequently lead to overestimation of the peak streamflow in May. Figs. 14 and 15 also show that the SD assimilation runs (especially SCF + SD_DA) have led to substantially larger improvements than the SCF_DA run on the mean RMSE and bias values in SD estimation across the contributing areas of Gages 6–8. The joint SCF and SD assimilation (SCD + SD_DA) has also resulted in large improvements over the OL run for Gages 4, 5, and 9. Interestingly, in these cases when SD_DA is not performing well (at Gages 4, 5, and 9), assimilating SCF in addition to assimilating SD does lead to substantially improved RMSE and bias values. This initially appears to be contrary to the previous observations that the SD assimilation tends to dominate the joint SCF + SD assimilation process. However, further investigation indicates that these three gages, unlike the other gages, have a considerable portion (i.e., 13%, 13%, and 23% for Gages 4, 5, and 9, respectively) of their contributing grid cells located in areas with an elevation above 1500 m, where SD assimilation is deactivated and consequently SCF assimilation becomes the only mechanism for updating the snow states in joint SCF + SD assimilation. The overall encouraging results from SD_DA and SCF + SD_DA shown in Figs. 14 and 15

223

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

Fig. 13. Mean monthly discharges from USGS observations and the OL, SCF_DA and SCf + SD_DA runs for the nine selected gages.

Table 5 Comparison of RMSE and NSE of the mean monthly discharge for each gage. Gage 1

Gage 2

Gage 3

Gage 4

Gage 5

Gage 6

Gage 7

Gage 8

Gage 9

RMSE (CFS)

OL MODIS_DA SCF_DA SD_DA SCF+SD_DA

1377 1364 1138 2317 2009

738 732 644 1193 1094

220 220 216 285 276

13486 13370 12356 16731 16226

5713 5695 5641 6956 6871

5502 5834 4528 5246 5142

1227 1230 1207 1416 1305

953 951 635 639 635

2727 2728 1830 2375 2202

[RMSE(OL)–RMSE(DA)]/RMSE(OL) (%)

MODIS_DA SCF_DA SD_DA SCF + SD_DA

1 17 68 46

1 13 62 48

0 2 30 26

1 8 24 20

0 1 22 20

6 18 5 7

0 2 15 6

0 33 33 33

0 33 13 19

NSE (–)

OL MODIS_DA SCF_DA SD_DA SCF + SD_DA

0.06 0.08 0.36 1.65 0.99

0.31 0.29 0.00 2.41 1.87

0.49 0.48 0.43 1.50 1.35

0.54 0.55 0.61 0.29 0.33

0.64 0.65 0.65 0.47 0.49

0.20 0.35 0.19 0.09 0.04

0.18 0.18 0.21 0.09 0.08

0.57 0.56 0.31 0.29 0.31

0.53 0.53 0.79 0.64 0.69

NSE(DA)–NSE(OL) (–)

MODIS_DA SCF_DA SD_DA SCF + SD_DA

0.02 0.30 1.71 1.05

0.02 0.31 2.10 1.56

0.00 0.06 1.01 0.86

0.01 0.07 0.25 0.21

0.00 0.01 0.17 0.15

0.15 0.39 0.12 0.16

0.00 0.03 0.27 0.10

0.01 0.87 0.86 0.87

0.00 0.26 0.11 0.16

Note bold values indicate the best performance.

suggest that snow prediction via assimilating bias-adjusted ANSA SD estimates is successful; however, improved snow prediction does not necessarily guarantee an improvement in snowmelt prediction or translate correspondingly into improvement in streamflow prediction. Other factors that can potentially contribute to

the unsatisfactory performance of SD assimilation in streamflow prediction include possible deficiencies in the model’s mechanism for translating snowmelt to runoff and the simplistic approach used in the present study for routing the runoff to the gages. Additionally, the current snow assimilation system (for either SCF or SD

224

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

Fig. 14. Mean RMSE in SD across the drainage area of each streamflow gage for April–June for open loop and assimilation runs.

assimilation) does not consider in the runoff generation process the amount of water that has been added to or removed from the system (as snow accumulation or melt) via assimilation, which interrupts the water balance of the hydrologic system (e.g., [17]) and may also have complications for the streamflow prediction process.

5. Summary and discussions In this study, influences of assimilating the MODIS SCF product and the SCF and SD estimates from the ANSA product are analyzed and discussed. An overall conclusion is that considerable improvements in snow and streamflow predictions in Alaska can be achieved via assimilating satellite-based SCF and/or SD products, provided that these products are properly preprocessed. As far as assimilating SCF is concerned, the interpreted ANSA SCF estimates tend to result in a substantially larger improvement in snow and streamflow predictions than the standard MODIS SCF estimates, owing to the considerably reduced cloud coverage and improved snow mapping accuracy in the ANSA SCF from applying the snow interpretation algorithm. Reducing the cloud cover in the SCF dataset is critically important for the Alaska region since snow melting typically lasts for only one week to 10 days and having SCF observations available in this short period is much advantageous. The rule-based assimilation strategy utilizes this information to approximate the melt rate and has effectively improved the modeling of the snow depletion process occurring in the region, resulting improvements on not only

predicting the snow amount, but also on predicting the snowmelt that effectively translates into improvement in streamflow prediction at all gage locations, with substantial improvements achieved at some locations. After bias adjustment against the snow observations from the SNOTEL and COOP stations, the ANSA SD estimates represent an overall improvement on the original ANSA SD estimates in the study region, especially over deep snowpack areas at high elevations (>500 m). This has resulted in an overall superior performance of SD assimilation in predicting SD in the study domain, especially in the coastal mountain areas. However, SD assimilation also introduces some overestimation in low-elevation areas (e.g., the northeast corner of the domain) in April and May, and subsequent overestimation in streamflow prediction in these areas. Unlike the rule-based strategy used for SCF assimilation, the EnKF-based assimilation of SD is easily impacted by the bias in the SD data and has resulted in the non-uniform performance of SD assimilation in snow prediction within the study domain (and likely unsatisfactory prediction on snowmelt), which has in turn led to improvements in streamflow prediction at some gage locations but degradations at the other locations. The results also suggest that when joint SD and SCF assimilation is implemented, the updating of snow states tends to be dominated by SD assimilation, leading to very similar results as those from assimilating SD only. However, in very high elevations where the quality of satellite-based SD estimates is in doubt, SCF assimilation can be used to improve the results. The current study is subject to a few limitations that need to be addressed in future work. For example, in the rule-based SCF

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

225

Fig. 15. Same as Fig. 14 but for the mean bias in SD.

assimilation, the model snowpack is prescribed to melt out in 10 days everywhere in the domain when the observation indicates no snow, while at some low elevations, the melting process may take less than 10 days. Using localized melting rates that are more relevant to the current ground and near surface conditions (e.g., air temperature and net radiation) should lead to improved assimilation results. In addition, since the current rule-based strategy is designed based on the observation that the Noah model mostly overestimates the snow amount in relatively high elevations in the study domain, it may not work well for other areas where the model underestimates the snowpack. Hence, it would be worthwhile to also explore other less location/model dependent approaches such as the forward-looking SCF assimilation strategy presented in Zaitchik and Rodell ([17]) that uses SCF observations up to 72 h ahead of the simulation time to adjust the major forcing variables (i.e., air temperature and precipitation) of the land surface model for improving predictions on snow cover extent and snow water equivalent. Moreover, it is noted that although the cloud cover is considerably reduced in the ANSA SCF dataset by the snow interpretation algorithm, extensive cloud cover still exists in the study domain; hence, it is desirable to explore other cloud cover removal approaches such as the spatial and temporal interpolation techniques discussed in Lo9 pez-Burgos et al. ([63]). Using both SNOTEL and COOP observations to adjust the bias in the ANSA SD estimates, although resulting in an overall improvement upon the original ANSA SD estimates (especially for high elevations), has introduced overestimation in low eleva-

tions in April and May and hence overestimation of streamflow in some low-lying basins. Future research should explore how the bias correction algorithm can be enhanced to achieve more consistent improvement in snow prediction over the entire domain. In other words, efforts should be made to avoid situations where improvement is achieved for parts of the domain while degradation is observed in some other parts. This is critical to ensure that the improvement in SD can be effectively translated into enhanced prediction of streamflow, which represents integration over the entire basin. In the current study, this is confirmed by the fairly consistent positive improvement across the domain from assimilating the ANSA SCF dataset that has led to improvement in streamflow prediction at every gage location. Future work could explore how excluding the SNOTEL observations from the bias adjustment of SD in low elevations (especially during the melt season) can help improve snow and streamflow predictions in these low elevation areas. Presumably, information on the elevation of the SNOTEL stations and the distribution of the snowline elevations (if available) can be useful in determining an appropriate elevation threshold for excluding SNOTEL observations for bias adjustment; this, however, is not always straightforward given the spatial variability in the snowline distribution and complications added by climate/weather patterns, topography, land cover etc. Alternatively, other sources of information such as the MODIS SCF product may also be used to constrain the bias correction to prevent over-adjusting the satellite-based SD or SWE estimates.

226

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227

The improvement in snow prediction and contrary degradation in streamflow prediction at some gage locations (e.g., Gage 5) by assimilating SD, as shown in Figs. 13–15, also suggest that snowmelt may have not been effectively modeled at these locations. Future work should consider assimilating satellite-based snow data to update the increments in the snow states (i.e., snow accumulation and melt) rather than the accumulated SWE/SD themselves, given that the uncertainty is really in the increments rather than the accumulated quantities. Such an approach might improve snow accumulation and melt in a more effective way than the current approach, leading to more effective improvements in streamflow prediction. In addition to bias-adjusting the satellite products, it is also worthwhile to explore the assimilation of the anomalies from these products to reduce the impact of bias (e.g., [15]). In addition, radiance-based snow assimilation has shown some promise in improving snow prediction (e.g., [29,31]) and can be explored for improving snow data assimilation and prediction in Alaska. In the present study, an approximation is used to compute the monthly streamflow values at the gages for evaluating the influence of snow assimilation on streamflow prediction; this may have introduced errors especially for gages with a large drainage area while also preventing the evaluation of streamflow prediction at finer than monthly (e.g., weekly or daily) scales. Hence, to use the streamflow observations more effectively, an appropriate routing scheme based on cell-to-cell routing or source-to-sink routing algorithms (e.g., [64]) should be used to properly transport the gridded runoff simulated by the model to the desired gage locations. Finally, in all the assimilation experiments explored in this study, runoff is not considered as a state variable of the model and hence is not updated along with SWE and SD during the assimilation process. In other words, snow amounts that are added to or removed from the system by the assimilation are not used to adjust the runoff amounts to achieve a better conservation of the mass balance. This may have implications for using the streamflow observations for evaluation of snow assimilation results. Li et al. [65] showed how the consideration of mass conservation in assimilating AMSR-E soil moisture retrievals can improve soil moisture prediction. It would be worthwhile to investigate how a similar mass conservation approach considering runoff in snow assimilation applications would influence the overall results, in particular in streamflow prediction. We expect future work along these lines will advance our understanding of the snow hydrology in Alaska and further improve the snow and streamflow predictions in this region.

Acknowledgement Support for this study was provided by NASA under Grant NNX08AU51G and as part of NASA’s contribution to the National Climate Assessment program. Additional funding comes from NOAA and the Air Force Weather Agency (AFWA). Comments from Martyn Clark and an anonymous reviewer helped to improve this manuscript.

References [1] Blanchet J, Marty C, Lehning M. Extreme value statistics of snowfall in the Swiss Alpine region. Water Resources Research 2009;45(5):1–12. http:// dx.doi.org/10.1029/2009WR007916. [2] Matsumura S, Yamazaki K, Tokioka T. Summertime land-atmosphere interactions in response to anomalous springtime snow cover in northern Eurasia. Journal of Geophysical Research 2010;115(D20):1–13. http:// dx.doi.org/10.1029/2009JD012342. [3] Wang Z, Zeng X, Decker M. Improving snow processes in the Noah land model. Journal of Geophysical Research 2010;115(D20):1–16. http://dx.doi.org/ 10.1029/2009JD013761.

[4] Livneh B, Xia Y, Mitchell KE, Ek MB, Lettenmaier DP. Noah LSM snow model diagnostics and enhancements. Journal of Hydrometeorology 2010;11(3):721–38. http://dx.doi.org/10.1175/2009JHM1174.1. [5] Barlage M, Chen F, Tewari M, et al. Noah land surface model modifications to improve snowpack prediction in the Colorado Rocky Mountains. Journal of Geophysical Research 2010;115(D22):1–15. http://dx.doi.org/10.1029/ 2009JD013470. [6] Andreadis KM, Storck P, Lettenmaier DP. Modeling snow accumulation and ablation processes in forested environments. Water Resources Research 2009;45(5):1–13. http://dx.doi.org/10.1029/2008WR007042. [7] Franz KJ, Butcher P, Ajami NK. Addressing snow model uncertainty for hydrologic prediction. Advances in Water Resources 2010;33(8):820–32. http://dx.doi.org/10.1016/j.advwatres.2010.05.004. [8] Rutter N, Essery R, Pomeroy J, et al. Evaluation of forest snow processes models (SnowMIP2). Journal of Geophysical Research 2009;114(D6). http://dx.doi.org/ 10.1029/2008JD011063. [9] He M, Hogue TS, Franz KJ, Margulis SA, Vrugt JA. Corruption of parameter behavior and regionalization by model and forcing data errors: a Bayesian example using the SNOW17 model. Water Resources Research 2011;47(7):1–17. http://dx.doi.org/10.1029/2010WR009753. [10] He M, Hogue TS, Franz KJ, Margulis SA, Vrugt JA. Characterizing parameter sensitivity and uncertainty for a snow model across hydroclimatic regimes. Advances in Water Resources 2011;34(1):114–27. http://dx.doi.org/10.1016/ j.advwatres.2010.10.002. [11] Kolberg S, Gottschalk L. Interannual stability of grid cell snow depletion curves as estimated from MODIS images. Water Resources 2010;46:1–15. http:// dx.doi.org/10.1029/2008WR007617. [12] Raleigh MS, Lundquist JD. Comparing and combining SWE estimates from the SNOW-17 model using PRISM and SWE reconstruction. Water Resources Research 2012;48(1):1–16. http://dx.doi.org/10.1029/2011WR010542. [13] Andreadis KM, Lettenmaier DP. Assimilating remotely sensed snow observations into a macroscale hydrology model. Advances in Water Resources 2006;29(6):872–86. http://dx.doi.org/10.1016/ j.advwatres.2005.08.004. [14] Clark MP, Slater AG, Barrett AP, et al. Assimilation of snow covered area information into hydrologic and land-surface models. Advances in Water Resources 2006;29(8):1209–21. http://dx.doi.org/10.1016/ j.advwatres.2005.10.001. [15] De Lannoy GJM, Reichle RH, Arsenault KR, et al. Multiscale assimilation of advanced microwave scanning radiometer–EOS snow water equivalent and moderate resolution imaging spectroradiometer snow cover fraction observations in northern Colorado. Water Resources Research 2012;48(1):1–17. http://dx.doi.org/10.1029/2011WR010588. [16] Su H, Yang Z-L, Dickinson RE, Wilson CR, Niu G-Y. Multisensor snow data assimilation at the continental scale: the value of gravity recovery and climate experiment terrestrial water storage information. Journal of Geophysical Research 2010;115(D10):1–14. http://dx.doi.org/10.1029/2009JD013035. [17] Zaitchik BF, Rodell M. Forward-looking assimilation of MODIS-derived snowcovered area into a land surface model. Journal of Hydrometeorology 2009;10(1):130–48. http://dx.doi.org/10.1175/2008JHM1042.1. [18] Slater AG, Clark MP. Snow data assimilation via an ensemble Kalman filter. Journal of Hydrometeorology 2006;7:478–93. [19] Yatheendradas S, Peters-Lidard CD, Koren V, et al. Distributed assimilation of satellite-based snow extent for improving simulated streamflow in mountainous, dense forests: an example over the DMIP2 western basins Soni Yatheendradas. Water Resources Research. in press. doi:10.1029/ 2011WR011347. [20] Barrett A, National operational hydrologic remote sensing center snow data assimilation system (SNODAS) products at NSIDC. NSIDC Special Report 11. Boulder, CO., USA; 2003:19. [21] Brasnett B. A global analysis of snow depth for numerical weather prediction. Journal of Applied Meteorology 1999;38:726–40. [22] Clark MP, Hendrikx J, Slater AG, et al. Representing spatial variability of snow water equivalent in hydrologic and land-surface models: a review. Water Resources Research 2011;47(7). http://dx.doi.org/10.1029/2011WR010745. [23] Wood EF, Roundy JK, Troy TJ, et al. Hyperresolution global land surface modeling: meeting a grand challenge for monitoring Earth’s terrestrial water. Water Resources Research 2011;47(5):1–10. http://dx.doi.org/10.1029/ 2010WR010090. [24] van Dijk AIJM, Renzullo LJ. Water resource monitoring systems and the role of satellite observations. Hydrology and Earth System Sciences 2011;15(1):39–55. http://dx.doi.org/10.5194/hess-15-39-2011. [25] Forman B, Reichle RH, Rodell M. Assimilation of terrestrial water storage from GRACE in a snow-dominated basin. Water Resources Research 2012;48(1):1–14. http://dx.doi.org/10.1029/2011WR011239. [26] Hall DK, Riggs GA. Accuracy assessment of the MODIS snow products. Earth Science 2007;21:1534–47. http://dx.doi.org/10.1002/hyp. 6715. [27] Foster JL, Hall DK, Eylander JB, et al. A blended global snow product using visible, passive microwave and scatterometer satellite data. International Journal of Remote Sensing 2011;32(5):1371–95. [28] Kumar SV, Reichle RH, Harrison KW, et al. A comparison of methods for a priori bias correction in soil moisture data assimilation. Water Resources Research 2012;48(3):1–16. http://dx.doi.org/10.1029/2010WR010261. [29] Dechant C, Moradkhani H. Radiance data assimilation for operational snow and streamflow forecasting. Advances in Water Resources 2011;34(3):351–64. http://dx.doi.org/10.1016/j.advwatres.2010.12.009.

Y. Liu et al. / Advances in Water Resources 54 (2013) 208–227 [30] Durand M, Kim EJ, Margulis SA. Radiance assimilation shows promise for snowpack characterization. Geophysical Research Letters 2009;36(2):1–5. http://dx.doi.org/10.1029/2008GL035214. [31] Toure AM, Goïta K, Royer A, et al. A case study of using a multilayered thermodynamical snow model for radiance assimilation. IEEE Transactions on Geoscience and Remote Sensing 2011;49(8):2828–37. [32] Liu Y, Gupta HV. Uncertainty in hydrologic modeling: toward an integrated data assimilation framework. Water Resources 2007;43:1–18. http:// dx.doi.org/10.1029/2006WR005756. [33] Reichle RH. Data assimilation methods in the Earth sciences. Advances in Water Resources 2008;31(11):1411–8. http://dx.doi.org/10.1016/ j.advwatres.2008.01.001. [34] Crow WT, Reichle RH. Comparison of adaptive filtering techniques for land surface data assimilation. Water Resources Research 2008;44(W08423). [35] Evensen G. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research 1994;99(C5):10143–62. http://dx.doi.org/10.1029/ 94JC00572. [36] Dong J, Walker JP, Houser PR, Sun C. Scanning multichannel microwave radiometer snow water equivalent assimilation. Journal of Geophysical Research 2007;112(D7):1–16. http://dx.doi.org/10.1029/2006JD007209. [37] Clark MP, Rupp DE, Woods RA, et al. Hydrological data assimilation with the ensemble Kalman filter: use of streamflow observations to update states in a distributed hydrological model. Advances in Water Resources 2008;31(10):1309–24. http://dx.doi.org/10.1016/j.advwatres.2008.06.005. [38] Rodell M, Houser PR. Updating a land surface model with MODIS-derived snow cover. Journal of Hydrometeorology 2004;5:1064–75. [39] Leisenring M, Moradkhani H. Snow water equivalent prediction using Bayesian data assimilation methods. Stochastic Environmental Research and Risk Assessment 2010;25(2):253–70. http://dx.doi.org/10.1007/s00477-010-04455. [40] Lee S, Klein AG, Over TM. A comparison of MODIS and NOHRSC snow-cover products for simulating streamflow using the snowmelt runoff model. Hydrological Processes 2005;19(15):2951–72. http://dx.doi.org/10.1002/hyp. 5810. [41] Déry SJ, Salomonson VV, Stieglitz M, Hall DK, Appel I. An approach to using snow areal depletion curves inferred from MODIS and its application to land surface modelling in Alaska. Hydrological Processes 2005;19(14):2755–74. http://dx.doi.org/10.1002/hyp. 5784. [42] Hall DK, Riggs GA, Foster JL, Kumar SV. Development and evaluation of a cloudgap-filled MODIS daily snow-cover product. Remote Sensing of Environment 2010;114(3):496–503. http://dx.doi.org/10.1016/j.rse.2009.10.007. [43] Corbari C, Ravazzani G, Martinelli J, Mancini M. Elevation based correction of snow coverage retrieved from satellite images to improve model calibration. Hydrology and Earth System Sciences 2009;13(5):639–49. http://dx.doi.org/ 10.5194/hess-13-639-2009. [44] Ek MB. Implementation of Noah land surface model advances in the national centers for environmental prediction operational mesoscale Eta model. Journal of Geophysical Research 2003;108(D22):1–16. http://dx.doi.org/10.1029/ 2002JD003296. [45] Chen F, Dudhia J. Coupling an advanced land surface – hydrology model with the penn state – NCAR MM5 modeling system. Part I: Model implementation and sensitivity. Monthly Weather Review 2001;129:569–85. [46] Koren V, Schaake J, Mitchell K, et al. A parameterization of snowpack and frozen ground intended for NCEP weather and climate models. Journal of Geophysical Research 1999;104(D16):19,569–85. [47] Sturm M, Holmren J, Liston GE. A seasonal snow cover classification system for local to global applications. Journal of Climate 1995;8:1261–83.

227

[48] Hartmann B, Wendler G. The significance of the 1976 Pacific climate shift in the climatology of Alaska. Weather 2005;18:4824–39. [49] Zhang T, Stamnes K, Bowling SA. Impact of the atmospheric thickness on the atmospheric downwelling longwave radiation and snowmelt under clear-sky conditions in the Arctic and subarctic. Journal of Climate 2001;14:920–39. [50] Brown RD, Brasnett B. Canadian meteorological centre (CMC) daily snow depth analysis data. Boulder, Colorado USA: Ó Environment Canada; 2010: Digital media. [51] Schwanghart W, Kuhn NJ. TopoToolbox: a set of Matlab functions for topographic analysis. Environmental Modelling and Software 2010;25(6):770–81. http://dx.doi.org/10.1016/j.envsoft.2009.12.002. [52] Kumar S, Peterslidard C, Tian Y, et al. Land information system: an interoperable framework for high resolution land surface modeling. Environmental Modelling and Software 2006;21(10):1402–15. http:// dx.doi.org/10.1016/j.envsoft.2005.07.004. [53] Kumar SV, Reichle RH, Peters-Lidard CD, et al. A land surface data assimilation framework using the land information system: description and applications. Advances in Water Resources 2008;31(11):1419–32. http://dx.doi.org/ 10.1016/j.advwatres.2008.01.013. [54] Rodell M, Houser PR, Jambor U, et al. The global land data assimilation system. Bulletin of the American Meteorological Society 2004;85(3):381–94. http:// dx.doi.org/10.1175/BAMS-85-3-381. [55] Cosgrove BA. Real-time and retrospective forcing in the North American land data assimilation system (NLDAS) project. Journal of Geophysical Research 2003;108(D22). http://dx.doi.org/10.1029/2002JD003118. [56] Salomonson V, Appel I. Estimating fractional snow cover from MODIS using the normalized difference snow index. Remote Sensing of Environment 2004;89(3):351–60. http://dx.doi.org/10.1016/j.rse.2003.10.016. [57] Gao Y, Xie H, Lu N, Yao T, Liang T. Toward advanced daily cloud-free snow cover and snow water equivalent products from Terra–Aqua MODIS and Aqua AMSR-E measurements. Journal of Hydrology 2010;385(1–4):23–35. http:// dx.doi.org/10.1016/j.jhydrol.2010.01.022. [58] Riggs G, Hall D, Salomonson V. MODIS snow products user guide collection. 2006;5:80. [59] Chang ATC, Foster JL, Hall DK. Nimbus-7 derived global snow cover parameters. Annals of Glaciology 1987;9:39–44. [60] Kelly REJ. The AMSR-E snow depth algorithm: description and initial results. Journal of the Remote Sensing Society of Japan 2009;29:307–17. [61] Sturm M, Taras B, Liston GE, et al. Estimating snow water equivalent using snow depth data and climate classes. Journal of Hydrometeorology 2010;11(6):1380–94. http://dx.doi.org/10.1175/2010JHM1202.1. [62] Peters-Lidard CD, Kumar SV, Mocko DM, Tian Y. Estimating evapotranspiration with land data assimilation systems. Hydrological Processes 2011;25(26):3979–92. http://dx.doi.org/10.1002/hyp. 8387. [63] López-Burgos V, Gupta HV, Clark M. A probability of snow approach to removing cloud cover from MODIS snow cover area products. Hydrology and Earth System Sciences Discussions 2012;9(12):13693–728. http://dx.doi.org/ 10.5194/hessd-9-13693-2012. [64] Zaitchik BF, Rodell M, Olivera F. Evaluation of the global land data assimilation system using global river discharge data and a source-to-sink routing scheme. Water Resources Research 2010;46(6):1–17. http://dx.doi.org/10.1029/ 2009WR007811. [65] Li B, Toll D, Zhan X, Cosgrove B. Improving estimated soil moisture fields through assimilation of AMSR-E soil moisture retrievals with an ensemble Kalman filter and a mass conservation constraint. Hydrology and Earth System Sciences 2012;16(1):105–19. http://dx.doi.org/10.5194/hess-16-105-2012.