Impact of Surface Albedo Assimilation on Snow Estimation

Abstract

Surface albedo has a significant impact in determining the amount of available net radiation at the surface and the evolution of surface water and energy budget components. The snow accumulation and timing of melt, in particular, are directly impacted by the changes in land surface albedo. This study presents an evaluation of the impact of assimilating Moderate Resolution Imaging Spectroradiometer (MODIS)-based surface albedo estimates in the Noah multi-parameterization (Noah-MP) land surface model, over the continental US during the time period from 2000 to 2017. The evaluation of simulated snow depth and snow cover fields show that significant improvements from data assimilation (DA) are obtained over the High Plains and parts of the Rocky Mountains. Earlier snowmelt and reduced agreements with reference snow depth measurements, primarily over the Northeast US, are also observed due to albedo DA. Most improvements from assimilation are observed over locations with moderate vegetation and lower elevation. The aggregate impact on evapotranspiration and runoff from assimilation is found to be marginal. This study also evaluates the relative and joint utility of assimilating fractional snow cover and surface albedo measurements. Relative to surface albedo assimilation, fractional snow cover assimilation is found to provide smaller improvements in the simulated snow depth fields. The configuration that jointly assimilates surface albedo and fractional snow cover measurements is found to provide the most beneficial improvements compared to the univariate DA configurations for surface albedo or fractional snow cover. Overall, the study also points to the need for improving the albedo formulations in land surface models and the incorporation of observational uncertainties within albedo DA configurations.

1. Introduction

Surface albedo, representing the reflectivity of a surface, is a key parameter in determining the land surface energy balance [1]. The net radiation on the surface is influenced by surface albedo as it controls the amount of reflected solar (shortwave) fluxes. Natural variability, as well as anthropogenic impacts, can influence albedo changes on the land surface [2,3]. The green-up from vegetation growth or precipitation on dry soils often leads to reductions in surface albedo whereas the presence of snow leads to significantly higher albedo in both vegetated and non-vegetated areas. Similarly, land use changes due to deforestation, agriculture, and natural hazards can lead to significant changes in vegetation characteristics and surface albedo. An increase in albedo leads to reduced shortwave absorption and reduced net radiation at the land surface. Over snow-covered areas, this can lead to reduced snowmelt, reduced runoff and increased snow-mass buildup. For example, the reduction in albedo from the deposition of airborne dust on glaciers has been shown to cause enhanced melting [4]. The reduced net radiation can also lead to reduced evaporation. Thus, the changes in surface albedo can have a propagating influence on the terrestrial water and energy budget components.

In recognition of its importance as a key modulating factor of the surface energy balance, albedo formulations in land surface models (LSMs) have evolved in complexity over the years [5]. Generally, these schemes have evolved from being vegetation based look-up tables to deriving albedo from canopy radiation transport formulations [1,6,7,8,9]. For example, the initial formulations of surface and snow-free albedos in the Noah LSM, used operationally at the US National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Prediction (NCEP), employ simple parameterizations. In this model version, the snow-free albedo values are either prescribed as a function of vegetation-based look-up tables or climatological inputs that limit the representation of interannual and seasonal variability in albedo [10]. The surface albedo is then parameterized as a function of vegetation seasonality and snow cover fraction [11]. Over deep snow, surface albedo is either set to a maximum value (a function of vegetation type) or an annual maximum snow albedo climatology. More recent changes to the Noah LSM have included time-varying decay formulations to model snow albedo [12], changes in fresh snow albedo and explicit computations of exposed roughness length in the presence of snow [13]. The Noah multi-parameterization (Noah-MP; [14,15] model, which represents the next generation of Noah LSM development, includes the calculation of albedo as a function of radiation transport within the canopy. A two-stream formulation for solar radiation exchange within the plant canopy [16,17] is used in Noah-MP, which computes the surface albedo as a function of the solar elevation, diffuse and direct solar radiation, vegetation properties, and bare soil/ground albedo. Despite these enhancements, the need to further improve the time-varying formulations of snow albedo and below-canopy radiative transfer for the accurate simulation of snow states have been emphasized in studies such as Chen et al. [18].

Remote sensing-based measurements of albedo have been derived from various geostationary and polar-orbiting satellites, relying on the angular distribution of reflectance measurements. Such retrievals have been generated from sensors such as Advanced Very High Resolution Radiometer (AVHRR; [19], Landsat Thematic mapper [20], Meteosat [21], and Moderate Resolution Imaging Spectroradiometer (MODIS; [22]. A review of available remote sensing-based broadband albedo products is provided in Qu et al. [23]. Though these retrievals provide valuable estimates of an essential climate system variable [24], they suffer from observational gaps related to cloud coverage, orbital configurations, and other sensing limitations. In addition, the spatial and temporal resolutions of these products also vary significantly. The MODIS-based albedo products, for example, are available at a temporal resolution of approximately eight days, which may not be adequate for characterizing the physical changes on the land surface.

Data assimilation (DA), which incorporates observational information to constrain model estimates, is an attractive option to blend the information from spatially and temporally infrequent observations with modeled estimates. The use of remote-sensing-based albedo observations for DA within LSMs is reported in several studies. For example, Malik et al. [25] examined the assimilation of MODIS-based snow albedo products over three selected locations in Colorado, within the Noah LSM. The assimilation was performed using a direct insertion approach, which led to improvements in simulated snow depth and duration of the snow season at these locations. A direct insertion approach to assimilate satellite-based albedo data into the ORCHIDEE LSM was also used by Wang et al. [26] to investigate the impacts on offline and coupled model simulations. They reported larger impacts in the coupled simulation from albedo assimilation due to the snow-albedo feedbacks. Yin et al. [27] examined the impact of using near real-time MODIS bidirectional distribution reflectance model albedo and green vegetation fraction within the Noah LSM. Their study focused on the evaluation of the net longwave and shortwave radiation at several locations in the continental US and reported positive impacts of assimilation. More recently, an optimal interpolation analysis of satellite-based surface albedo and leaf area index in a numerical weather prediction system was employed by Boussetta et al. [28]. This study reported small impacts on the surface fluxes and near-surface air temperature and humidity from albedo assimilation. It is notable that despite the direct linkages between surface albedo and snow evolution, only a few of these studies quantitatively examined the utility of albedo updates on snow simulation. Overall, there is limited work in the literature on the evaluation of the impacts of albedo assimilation at regional or continental scales.

Compared to the DA studies on incorporating albedo, there have been more efforts on the assimilation of snow cover fraction estimates from optical sensors, presumably because snow cover is a more direct measure of snow presence [29,30,31,32,33]. These studies employ heuristic approaches by correcting the mismatches in snow presence between the model and the observations either by removing the snow in the model or through the nominal increase of the model snow. The assimilation of snow cover estimates is most useful in the transitional snow periods, during snow onset and ablation. Over locations and times with extended periods of complete snow cover, the utility from snow cover fraction assimilation is small [34]. Though both surface albedo and snow cover estimates are available from optical sensors, very few efforts have evaluated their relative utility in assimilation. Xu and Shu [35] examined the impact of assimilating MODIS-based snow cover fraction and albedo within the Common Land Model over a small region in Northern China. In their study, the influence of albedo assimilation on the snow states was less impactful compared to that of the snow cover fraction assimilation. They attributed this finding to the indirect update of the snow variables from albedo assimilation. Similar assessments over large spatial extents, however, are currently lacking.

In this study, we present the assimilation of remotely sensed surface albedo measurements in the Noah-MP LSM over the continental US In addition, we also compare the relative utility of the surface albedo assimilation configuration to that of incorporating fractional snow cover ($fsc$) measurements. Spatially distributed estimates of albedo and $fsc$ from optical sensors are used for data assimilation. Specifically, we used the land surface albedo estimates [36] from the Global Land Surface Satellite (GLASS) project and the snow cover data from MODIS aboard the Terra satellite [37]. The GLASS product provides high-quality, gap-filled estimates of directional hemispherical reflectance (referred to as the black-sky albedo) and bi-hemispherical reflectance (known as the white-sky albedo). For the assimilation of $fsc$ estimates, the collection 5 (C5) version of the MODIS snow cover product (known as MOD10A1) was used. The assimilation integrations were conducted within the North American Land Data Assimilation System (NLDAS) configuration, which includes high-quality boundary conditions and land surface parameters.

The study explores the following three specific research objectives:

• Quantify the spatio–temporal impact of assimilating surface albedo estimates on the simulation of snow states.
• Evaluate the utility of albedo assimilation on improving runoff/streamflow estimation.
• Compare and contrast the utility of assimilating albedo inputs relative to that from the assimilation of $fsc$ estimates.

The first two objectives address a key gap in the current land data assimilation efforts as continental-scale studies of assimilating albedo retrievals have not been done in prior research. Since most studies on snow estimation using optical sensor-based data are focused on exploiting the information from $fsc$ retrievals, objective 3 listed above provides an assessment of the relative utility of $fsc$ and albedo measurements. The impact of assimilation on the simulated snow depth, snow cover, and streamflow is quantified through comparisons against a number of independent reference data products.

The article is organized as follows. Section 2 provides the details of the model domain, albedo formulations in Noah-MP, data assimilation setup, remote sensing and evaluation datasets. The analysis and discussion of results are presented in Section 3 and the summary and major conclusions are described in Section 4.

2. Study Settings

2.1. Model Configuration

All model simulations are conducted over a domain similar to the one used in the North American Land Data Assimilation System phase-2 (NLDAS-2) project, using the Noah-MP LSM version 3.6. The simulation domain spans the Continental United States (CONUS; 25–53${}^{\circ }$N and 125–67${}^{\circ }$W) at 1/8${}^{\circ }$ spatial resolution (Figure 1). The two locations A and B shown in Figure 1 are regions with different topographical characteristics and snow evolution patterns and are used for time series comparisons in Section 3. Site A is Cavalier 7NW, North Dakota, a site over the upper High Plains, which is at a low elevation (272 m) and where the snow evolution is typically ephemeral. Location B, Wild Horse Reservoir, Nevada, is over an area at a higher elevation (1892 m) with significant seasonal snowpacks.

The NLDAS-2 forcing dataset [38] that includes daily gauge-based precipitation (temporally disaggregated using radar estimates), bias-corrected shortwave radiation, and a surface meteorology analysis is used as boundary conditions for the LSM simulations. The model integrations are conducted using the NASA Land Information System (LIS; [39,40], a widely used, open-source land surface modeling and data assimilation system. LIS encompasses several land surface models, a comprehensive data assimilation subsystem, and support for a large suite of remote sensing data [41]. Through an interoperable design, LIS facilitates the configuration of sequential and non-sequential data assimilation configurations, as evidenced by numerous prior studies [42,43,44,45,46,47].

The model and data assimilation integrations are conducted during the time period from 2000 to 2017. The initial conditions for the LSM runs are generated through a long spinup of the Noah-MP model starting with uniform conditions based on the recommendations of [48]. Specifically, the model is run from 1979 to 2017 and restarted in 1979 using the initial conditions at the end of 2017 (essentially running the model for 78 years). A third run is conducted by reinitializing the model in 1979 using the climatological average conditions from this 78-year spinup. The third run is simulated without data assimilation from 1979 to 2000, after which data assimilation of the MODIS-based surface albedo and $fsc$ products began and continued to the end of 2017. In addition to Noah-MP, the Hydrological Modeling and Analysis Platform (HyMAP; [49] streamflow routing model is used to develop estimates of streamflow. The influence of lakes and reservoirs and the impacts of reservoir operations are not modeled in this study.

2.2. Noah-MP Land Surface Model

Noah-MP, developed to augment the Noah LSM, includes multiple physics options for surface water infiltration, runoff, groundwater, vegetation phenology, and multilayer snowpack, among others. Noah-MP also includes a number of modeling enhancements to address the early snowmelt and the shallow snowpacks in the traditional Noah LSM [12,50]. A modified two-stream radiation scheme [17] with representation of shading effects of the vegetation canopy is used to model the surface albedo, fluxes, and vegetation on the land surface. Among the two options available for modeling snow surface albedo in Noah-MP, a formulation from the Canadian Land Surface Scheme (CLASS) model [51] is used in this study. The snow albedo (${\alpha }_{snow}$) is modeled as a function of the snow surface albedo from the previous timestep (${\alpha }_{snow}^{-}$), the density of the fresh snow ($f_{sn}$) and the model timestep ($dt$). Note that the CLASS scheme does not account for the effects of snow aging and impurities on snow and albedo.

$$\begin{array}{cc}\hfill {\alpha }_1& =0.55+({\alpha }_{snow}^{-}-0.55)e^{-0.01dt}\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\alpha }_{snow}& ={\alpha }_1+f_{sn}dt(0.84-{\alpha }_1).\hfill \end{array}$$

(2)

The direct and diffuse components of ground surface albedo for visible (VIS) and near-infrared (NIR) bands (${\alpha }_g^{d,\mathrm{VIS}/\mathrm{NIR}}$ and ${\alpha }_g^{i,\mathrm{VIS}/\mathrm{NIR}}$) are then parameterized as a weighted average of the snow albedo and bare soil albedo using the $fsc$ estimate from the model. The bare soil albedo is computed separately for the direct (${\alpha }_b^{d,\mathrm{VIS}/\mathrm{NIR}}$) and diffuse (${\alpha }_b^{i,\mathrm{VIS}/\mathrm{NIR}}$) components as a function of the saturated soil albedo and the dry soil albedo, prescribed as a look-up table of soil color.

$$\begin{array}{cc}\hfill {\alpha }_g^{d,\mathrm{VIS}/\mathrm{NIR}}& ={\alpha }_{snow}^{d}fsc+{\alpha }_b^{d,\mathrm{VIS}/\mathrm{NIR}}(1-fsc)\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\alpha }_g^{i,\mathrm{VIS}/\mathrm{NIR}}& ={\alpha }_{snow}^{i}fsc+{\alpha }_b^{i,\mathrm{VIS}/\mathrm{NIR}}(1-fsc).\hfill \end{array}$$

(4)

A modified version of the two-stream radiative transfer model is used to calculate the radiation transfer through the canopy as well as the total albedo for the downwelling direct (${\alpha }^{d,\mathrm{VIS}/\mathrm{NIR}}$) and diffuse (${\alpha }^{i,\mathrm{VIS}/\mathrm{NIR}}$) components, for visible and near-infrared bands [17].

The total albedo components are used to calculate the reflected radiative fluxes from the visible ($R_{\mathrm{VIS}}$) and near infrared ($R_{\mathrm{NIR}}$) bands as:

$$\begin{array}{cc}\hfill R_{\mathrm{VIS}}& ={\alpha }^{d,\mathrm{VIS}}\ast SWdow{n}_{VIS}^d+{\alpha }^{i,\mathrm{VIS}}\ast SWdow{n}_{VIS}^i\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill R_{\mathrm{NIR}}& ={\alpha }^{d,\mathrm{NIR}}\ast SWdow{n}_{NIR}^d+{\alpha }^{i,\mathrm{NIR}}\ast SWdow{n}_{NIR}^i,\hfill \end{array}$$

(6)

where $SWdow{n}_{VIS}^d$, $SWdow{n}_{VIS}^i$, $SWdow{n}_{NIR}^d$ and $SWdow{n}_{NIR}^i$ represents the direct visible, diffuse visible, direct near-infrared and diffuse near-infrared downwelling radiation, respectively. The total albedo is then calculated as the ratio of the reflected radiative flux ($R_{VIS}+R_{NIR}$) to the total downward solar radiation ($SWdown$).

2.3. Remote Sensing Data

As noted earlier, the gap-filled black-sky and the white-sky albedo estimates from the GLASS project and the MOD10A1 $fsc$ measurements are employed for data assimilation. The black-sky component is the albedo in the absence of a diffuse component, whereas the white-sky component represents the albedo in the absence of the direct component. The GLASS products generated from the MODIS Terra and Aqua surface reflectance measurements are used in this study. Prior studies that compared the GLASS albedo estimates against ground-based and other remote sensing-based estimates quantify the high-quality of the product for different aerosol conditions, landcover regimes and meteorological conditions [52,53]. Generally, the GLASS albedo estimates have low errors under clear sky, snow-free conditions whereas the uncertainty increases under cloudy, snow/ice conditions. The GLASS albedo product has been documented to be of comparable quality to that of the standard MODIS product MCD43 [53]. The gap-filled nature of the GLASS product offers an advantage over MCD43 for modeling and assimilation purposes. The GLASS datasets are available at eight-day intervals on a 0.05${}^{\circ }$ regular latitude–longitude global grid.

The $fsc$ retrievals based on the measurements from the MODIS instrument on the Terra spacecraft are used for snow cover assimilation. Specifically, we use the MOD10A1 collection 5 [54] data product. The $fsc$ retrievals in MOD10A1 is generated based on the algorithm of Salomonson and Appel [55]. The original data available at 500 m is resampled to 1 km spatial resolution for the purposes of this study. The MOD10A1 data files are available daily from February 2000 to present.

2.4. Data Assimilation Configuration

As described in Section 2.2, the total albedo formulations in Noah-MP are diagnostic in that the ${\alpha }^d$ and ${\alpha }^i$ terms are calculated every time step as a function of ground surface albedo, presence of snow cover, and vegetation, among other factors. As a result, at every timestep, the albedo values generated by Noah-MP default back to the values generated as a function of the model look-up tables and default model parameters in the absence of observational constraints. In order for the assimilation updates to have a persistent impact in the model simulations, the GLASS albedo estimates temporally interpolated to the model timestep are used for data assimilation. At each timestep, the albedo inputs are generated by linearly interpolating between the GLASS albedo data available at every eight days. This is similar to the practice of using temporally interpolated values from monthly albedo climatology with the standard Noah model, employed routinely in research and application configurations [11]. The assimilation of the interpolated albedo values is performed using a direct insertion approach. Note that direct insertion is used here because the albedo formulations in Noah-MP are diagnostic. The availability of a prognostic albedo representation will allow the use of advanced DA algorithms that incorporate the relative errors in models and observations.

The sequence of updates when surface albedo retrievals is assimilated in Noah-MP is shown in Figure 2. The model physics calculations described above are performed first, for each grid point and at each timestep. This is followed by the calculation of the scaling factors representing the ratio of the ground surface albedo to the total albedo, for both direct and diffuse components.

$$\begin{array}{cc}\hfill {\gamma }^{d,\mathrm{VIS}/\mathrm{NIR}}& =\frac{{\alpha }_g^{d,\mathrm{VIS}/\mathrm{NIR}}}{{\alpha }^{d,\mathrm{VIS}/\mathrm{NIR}}}\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\gamma }^{i,\mathrm{VIS}/\mathrm{NIR}}& =\frac{{\alpha }_g^{i,\mathrm{VIS}/\mathrm{NIR}}}{{\alpha }^{i,\mathrm{VIS}/\mathrm{NIR}}}.\hfill \end{array}$$

(8)

The total and the ground surface albedo values are then updated using the input observations of black-sky (${\alpha }_o^{d,\mathrm{VIS}/\mathrm{NIR}}$) and diffuse (${\alpha }_o^{d,\mathrm{VIS}/\mathrm{NIR}}$) albedo values for visible and near-infrared bands, as follows.

$$\begin{array}{cc}\hfill {\alpha }_g^{d,\mathrm{VIS}/\mathrm{NIR}}& ={\gamma }^{d,\mathrm{VIS}/\mathrm{NIR}}\ast {\alpha }_o^{d,\mathrm{VIS}/\mathrm{NIR}}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\alpha }_g^{i,\mathrm{VIS}/\mathrm{NIR}}& ={\gamma }^{i,\mathrm{VIS}/\mathrm{NIR}}\ast {\alpha }_o^{i,\mathrm{VIS}/\mathrm{NIR}}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\alpha }^{d,\mathrm{VIS}/\mathrm{NIR}}& ={\alpha }_o^{d,\mathrm{VIS}/\mathrm{NIR}}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\alpha }^{i,\mathrm{VIS}/\mathrm{NIR}}& ={\alpha }_o^{i,\mathrm{VIS}/\mathrm{NIR}}.\hfill \end{array}$$

(12)

Using these updated values, the two-stream computations are repeated to reflect the impact of changes to the total and ground surface albedos on the flux calculations.

Snow cover assimilation is performed using the approach of Rodell and Houser [29], using MOD10A1 based $fsc$ estimates. The model and MOD10A1-based $fsc$ values are compared first. If the observation indicates no snow and the model indicates the presence of snow, the snow states within the model are reduced to zero. On the other hand, if the model has no snow and observation includes snow, then a nominal amount of snow (10 mm of snow water equivalent) is added to the model. Data assimilation is performed at 10.30 AM locally, on a gridcell-by-gridcell basis, to match the overpass time of the Terra satellite.

2.5. Observation Thinning for DA

Thinning of the observational datasets to ensure the use of high-quality observations is an essential step when satellite datasets are employed in data assimilation. Since direct insertion approaches that preclude the specification of the relative errors and uncertainties in the observations are used here for DA, it becomes more imperative to develop reliable screening procedures for assimilated datasets. As prior studies have established that the albedo retrieval uncertainty increases in the presence of cloud cover and thick vegetation [56,57,58,59], observational thinning procedures reflecting these limitations are applied. First, GLASS albedo retrievals that are labeled with an uncertainty larger than 0.08 (which reflects the area proportion of 1 km sub-pixels with the good or acceptable quality flag within the 5 km pixel) are excluded from assimilation. In addition, observational screens based on the information from the LSM are used by excluding retrievals when the LSM indicates dense vegetation (when green vegetation fraction is greater than 0.8) or when the landcover type is either Evergreen needleleaf forest, Evergreen broadleaf forest, Deciduous needleleaf forest or Deciduous broadleaf forest. These observational screening thresholds are developed through sensitivity experiments to provide a reasonable compromise between eliminating too many observations from assimilation vs. allowing low-quality observations within assimilation.

2.6. Methods and Datasets Evaluated

Quantitative assessments of the impact of data assimilation of albedo and snow cover measurements were developed by comparing the modeled snow depth, snow cover, and streamflow fields against independent reference datasets. All evaluations were performed using the NASA Land surface Verification Toolkit (LVT; [60], a formal environment for land model skill assessment.

For the evaluation of modeled snow depth fields, we used three reference datasets: (1) the daily in-situ measurements from the Global Historical Climate Network (GHCN; [61], (2) the daily, gridded snow depth analysis from NOAA National Weather Service’s National Operational Hydrologic Remote Sensing Center (NOHRSC) Snow Data Assimilation System (SNODAS; [62], and (3) the daily gridded snow depth estimates developed by University of Arizona (UA; [63]). A selected, quality controlled subset of the GHCN stations, consistent with prior studies [43,46] is used here. SNODAS data products at 1 km spatial resolution were developed by assimilating satellite-derived, airborne, and ground-based observations of snow covered area and snow water equivalent (SWE) into a physically-based, multilayer energy and mass balance snow model. MODIS snow cover data was also used within SNODAS to make corrections to the presence or absence of snow (Carrie Olheiser, NOHRSC, personal communication). Limited evaluations of the SNODAS product [64,65] using regional in-situ or airborne datasets document the reasonable performance of SNODAS, with larger errors over alpine areas. In a comparison of several model products, Dawson et al. [66] note large errors in the snow density estimates from SNODAS, particularly over ephemeral and maritime climates. The UA analysis produced at 4 km spatial resolution was developed partly in response to the limitations of modeled snow analysis products. This product was developed using an empirical temperature index snow model with ground measurements from networks such as the National Resources Conservation Services’s SNOTEL and the National Weather Service’s COOP. Based on the availability of the datasets, the snow depth evaluations presented in Section 3.2.1 encompass time periods of 2000–2017, 2003–2017, and 2000–2016, for GHCN, SNODAS, and UA, respectively.

The simulated snow cover estimates are compared against the MOD10A1 product using categorical metrics of probability detection (POD) and false alarm ratio (FAR). POD is the fraction of snow cover presence that is correctly simulated whereas FAR represents the fraction of no-snow events that is incorrectly simulated. For the calculation of POD and FAR, only snow cover fraction values that are greater than 20% are considered. In addition, grid points in MOD10A1 with cloud obscuration were excluded in the computations. Finally, locations where less than 40% of the maximum available data was present during the 2000–2017 time period are also excluded.

The model simulated streamflow estimates are compared against the daily data obtained from the US Geological Survey over 572, small basins where the impact of reservoir operations is considered minimal [46]. Similar to the strategy used in prior studies, the impact from DA on streamflow is quantified through the Normalized Information Contribution (NIC) metrics [67]. The NIC of Nash Sutcliffe Efficiency (NSE) and RMSE represents normalized measures of improvements from DA with positive and negative NIC values indicating improvements and degradations from DA, respectively. The NIC metrics are defined as follows:

$$NI{C}_{RMSE}=\frac{(RMS{E}_o-RMS{E}_a)}{RMS{E}_o}$$

(13)

$$NI{C}_{NSE}=\frac{(NS{E}_a-NS{E}_o)}{(1-NS{E}_o)},$$

(14)

where the subscripts o and a denote open loop and assimilation, respectively.

3. Results and Discussion

This section focuses on the evaluation of the impact of updating the modeled surface albedo and snow cover with remote sensing information. The analysis focuses primarily on four separate model integrations: (1) an open-loop (OL) simulation that employs look-up table values of albedo and does not involve any assimilation, (2) a simulation that assimilates the GLASS total surface albedo datasets, (3) a simulation that assimilates the MODIS $fsc$ data, and (4) a simulation that jointly assimilates GLASS surface albedo and MODIS $fsc$ datasets. This model suite was designed to quantify the impact of assimilating albedo data, relative both to the model run with no assimilation as well as the snow cover assimilation integration. Configuration 4 was targeted to quantify the joint impact of incorporating both surface albedo and snow cover information.

The subsequent sections first provide a comparison of the spectral albedo estimates from Noah-MP against the GLASS albedo retrievals. This is followed by an assessment of the impact of changes in albedo, SWE, and surface water budget partition from surface albedo assimilation relative to the OL. Evaluation of the impact of data assimilation on modeled snow depth, snow cover, and streamflow are presented using the reference datasets described in Section 2.6. Finally, an evaluation of the relative and joint assessment of assimilating surface albedo and snow cover estimates is presented.

3.1. Evaluation of Simulated Spectral Albedo from Noah-MP

As noted in Section 2.4, surface albedo DA is enabled by assimilating the spectral components of surface albedo. We first compare the Noah-MP OL estimates of ${\alpha }^{d,\mathrm{VIS}}$, ${\alpha }^{d,\mathrm{NIR}}$, ${\alpha }^{i,\mathrm{VIS}}$, and ${\alpha }^{i,\mathrm{NIR}}$ against the GLASS data. Figure 3 shows the average RMSE of the spectral albedo components during the 2000–2017 time period. Overall, larger differences between the model and the GLASS retrievals are observed for the visible components, particularly over the northern parts of the domain. Comparatively, the RMSE patterns are more uniform in the NIR comparisons. It is well established that plant canopies and soils reflect a much larger fraction in the NIR component than in the VIS [68]. This is evident in the NIR RMSE albedo maps that demonstrate the strong influence of vegetation types. Patterns of lower RMSE values in NIR albedo are seen over the forested regions (parts of Northwest, Northeast, Appalachian Mountains) and the croplands (the Mississippi basin, and parts of the Central US), strongly correlated with the landcover map of Figure 1. In contrast, the influence of vegetation is less prominent in the VIS albedo RMSE maps. The larger influence of snow on the VIS albedo is seen in these comparisons, as high RMSE values in ${\alpha }^{d,\mathrm{VIS}}$ and ${\alpha }^{i,\mathrm{VIS}}$ are observed over the High Plains, parts of Midwest, and the Rocky mountains, areas with seasonal snowpacks.

3.2. Impact of Albedo Assimilation Schemes

Figure 4 shows the aggregate impact of the changes in albedo, SWE, runoff and ET as a result of the surface albedo assimilation across the 2000–2017 time period. The general impact of surface albedo DA is to increase albedo over most regions except over the Northeast, parts of the Northwest, and Midwest. The pattern of increased albedo over the High Plains mirrors the pattern of the larger errors in the VIS albedo components over this region, as shown in Figure 3. The increase in albedo allows for reduced net radiation at the surface, which leads to reduced snowmelt and increased snow accumulation, reflected in the map that shows the changes in SWE. For example, over the High Plains where assimilation increases albedo, larger snow accumulation is seen as a result of the reduction in the net radiation available for snowmelt. Over regions where albedo is reduced as a result of DA (such as the Northeast part of the domain), there is increased net radiation and reduced snow accumulation. The aggregate impact of albedo DA configurations is generally small on runoff and ET. The increase in albedo generally leads to increased snow and runoff and corresponding reductions in ET.

3.2.1. Snow Depth

Figure 5 shows maps of changes in RMSE (expressed as RMSE (Open Loop; OL) minus RMSE (DA)) in simulated snow depth from surface albedo assimilation using GHCN, SNODAS, and UA estimates as the reference. These comparisons indicate that the impact of albedo assimilation is generally positive, with improvements in simulated snow depth observed in several areas of the modeled domain. In all three comparisons, the snow depth estimates are improved over the High Plains, parts of the Midwest and the Rocky Mountains, whereas assimilation leads to degradations over the Northeast, parts of the western US, and areas near the Great Lakes. It is also notable that the patterns of improvements and degradations in the simulated snow depth fields are generally consistent in all three comparisons, though the magnitudes of these changes vary.

The propagation of the impact of assimilating the spectral albedo components are shown in the time series comparisons at two representative locations (A and B, shown in Figure 1 in Figure 6 and Figure 7, during the water years of 2014 and 2015. Figure 6 and Figure 7 show the time series of the biases in the Noah-MP OL spectral albedo components, the changes in total surface albedo as a result of DA, the changes in the reflected radiation components from DA, and the time series of snow depth from OL, DA, SNODAS, and GHCN.

As shown in Figure 6, at location A, the model simulated spectral albedo components are negatively biased relative to the GLASS estimates. The biases are largest in the winter months and small in the non-winter (April to October) time periods. Relative to the OL, the spectral albedo components, the total surface albedo, and the reflected solar radiation increase with data assimilation, particularly during the winter time periods. Mirroring the biases in the visible albedo components, DA leads to larger changes in the visible component of the reflected shortwave. The increase in the reflected solar radiation (and the reduction in net radiation) leads to larger snowpacks with DA, and stronger agreements with SNODAS and GHCN at this location. In particular, the early snowmelt in the model OL is significantly improved by surface albedo assimilation.

Similar changes in albedo, net radiation, and snow depth are observed at location B, where DA improves the snow depth simulations, particularly in the water year 2015. The snowpack time series at this location shows the historically dry winter in water year 2014 and the heavy precipitation that led to the reduction of drought in the water year 2015. The reference datasets indicate a large seasonal snowpack in 2015, which is significantly underestimated in the model simulations. The assimilation of albedo datasets helps in improving the early snowmelt in the model, by increasing the reflected solar radiation. There are larger differences in the albedo components, and the reflected radiation in 2015 compared to that of 2014. Similar to location A, the VIS albedo components dominate the changes in total albedo and the reflected solar radiation, which increase as a result of DA.

Though Figure 5 shows similar geographical patterns of improvements in the GHCN, SNODAS and UA comparisons, it should be noted that there are significant disagreements between these reference datasets at several locations, as seen in Figure 6 and Figure 7. In many of these sites, the spread across the reference datasets is larger than the changes in snow depth introduced by DA. Note also that the changes in albedo estimates from assimilation are not always in one direction alone. While there are certain locations and times when the reduction in albedo improves the snow simulation (not shown), the delay in snowmelt through an increase in albedo is generally beneficial for obtaining better agreements with the reference measurements. Conversely, the reduction in albedo, particularly over areas in the Northeast US, leads to reduced agreements with the reference measurements.

The impact of albedo data assimilation on the accuracy of simulated snow depth is examined further in Figure 8, where the RMSE of snow depth is stratified regionally and by the accumulation (SON), peak (DJF) and ablation (MAM) seasons. Overall, the largest errors are seen over the Northwest domain, where DA leads to marginal improvements in the GHCN comparisons and degradations in the SNODAS and UA comparisons (though these domain averaged changes are not statistically significant). In all regional comparisons, it can also be noted that the largest errors occur in the melt season, except over the Southeast domain. The largest statistically significant impact from DA is over the Northeast (degradation) and Great Plains and Midwest (improvements). In the GHCN comparisons, the impact of albedo assimilation is generally positive over all regions except the Northeast. Similarly, the SNODAS and UA comparisons also indicate positive impacts from albedo DA in all regions except the Northwest and Northeast. In most regions and in comparisons with all three reference datasets, DA provides improvements in the peak season (DJF), whereas most degradations and disagreements with the reference datasets occur in the melt (MAM) season. For example, over the Northwest, the overall degradation (in the SNODAS and UA comparisons) is a result of the disagreements in the MAM time period, where DA leads to increased RMSE in snow depth.

The influence of topography on the performance of the albedo DA schemes is evaluated by stratifying the RMSE estimates by elevation, as shown in Figure 9. The SNODAS and UA evaluations indicate that errors in the assimilation estimates increase at higher elevations, though the GHCN evaluations indicate the opposite. All three comparisons suggest that assimilation provides small, but beneficial improvements at lower elevations (<1500 m). As also shown in Figure 9, the coverage of the GHCN stations and SNODAS/UA grid points is highly non-uniform across these elevation ranges, with more than 80% of the GHCN stations at elevations less than 1500 m. The larger errors in the assimilation estimates at higher elevations are likely due to the increased uncertainties in the albedo retrievals from cloud cover. Note that vegetation also plays a significant role in determining the level of improvements obtained from albedo DA. Since a vegetation-based screening strategy is already used in the assimilation (Section 2.5, a similar stratified analysis based on vegetation type is omitted.

3.2.2. Snowcover

Figure 10 shows maps of the changes in POD and FAR from the surface albedo assimilation relative to the OL. The spatial patterns of improvements and degradations in POD are generally consistent with the patterns seen in the comparisons of snow depth relative to GHCN, SNODAS and UA. Most improvements from surface albedo assimilation are observed over the High Plains, and parts of the Rocky Mountains. Reduction in the POD of snow cover fraction estimates is observed over the Northeast, parts of Cascades, Sierra Nevada and Rocky Mountain ranges. Comparatively, the impacts on FAR from DA are smaller, particularly over areas with large snow evolution. Note that this evaluation against MOD10A1 is not truly independent, since the GLASS albedo products are also derived based on MODIS-based reflectances.

3.2.3. Streamflow

Snow conditions have a significant influence on the river runoff, particularly over the mid-latitude and high-latitude regions. For example, studies such as [69] have quantified that snowpack in the western US is the largest component of water storage. Since albedo assimilation has a more direct impact on snow accumulation, here we evaluate the subsequent impact on the simulated streamflow.

Figure 11 shows the $NI{C}_{NSE}$ and $NI{C}_{RMSE}$ from the two assimilation configurations. Overall, the assimilation of surface albedo provides minimal impact on streamflow, as the changes in NIC values are not statistically significant over majority of the USGS stations. Both assimilation configurations lead to degradations in streamflow over the Northwest and Northeast US, mirroring the degradations observed in the snowpack estimates. On the other hand, the surface albedo assimilation improves the streamflow simulations in the Missouri and upper Mississippi basins, likely a result of the improvements in snow over those regions.

3.3. Relative and Joint Impact of Snow Cover Assimilation

As noted earlier, since most studies assimilating optical sensor-based snow measurements are focused on the use of snow cover retrievals, here we evaluate the relative impact of albedo and $fsc$ assimilation. Using the snow cover assimilation strategy described in Section 2.4, DA integrations are conducted during the same time period as that of the albedo DA. The impact of snow cover assimilation on simulated snow depth is shown in Figure 12, which presents the changes in snow depth RMSE from DA (using GHCN, SNODAS, and UA datasets as the reference). Overall, snow cover assimilation has a positive impact on simulated snow depth, though the improvements are generally small relative to those from surface albedo assimilation. For example, while surface albedo assimilation has a strong feature of improvements over the High Plains, the snow cover assimilation has a more muted, yet beneficial impact over the same region, particularly in the GHCN and SNODAS comparisons. Relative to UA, there are little changes in the model skill from $fsc$ assimilation over the High Plains. Over the Northeast US, surface albedo assimilation reduces the agreement of simulated snow depth with SNODAS and GHCN, whereas snow cover assimilation has a neutral influence in this region. Over the mountainous Western US terrain, there are both positive and negative patterns of RMSE from the assimilation of surface albedo. The degradations from snow cover assimilation are small over these areas, though the improvements are also small.

As shown in Equations (3) and (4), fractional snow cover has a direct influence on the parameterization of the ground surface albedo (and the total direct and diffuse albedo) values. If the assimilation of snow cover data improves the representation of $fsc$, then it would also impact the characterization of the albedo, and the calculation of snow mass. This is examined by jointly assimilating both snow cover and surface albedo datasets within a multivariate DA configuration. The univariate assimilation schemes for MOD10A1 snow cover and GLASS surface albedo datasets are combined within a single DA instance. The impact of the joint assimilation on simulated snow depth is shown in Figure 13, evaluated using the three reference datasets. Figure 13 shows that, relative to the assimilation of surface albedo data, the joint assimilation configuration provides added improvements in snow depth simulation. There are larger improvements in the simulated snow depth over the High Plains and parts of the Rocky Mountains. In addition, there are marginal improvements in the disagreements over the Midwest, Western and the Northeast US between modeled and reference snow depth datasets. Some patterns of increased degradations in the joint assimilation configuration are also observed. For example, relative to the UA dataset, the joint assimilation increases the RMSE in snow depth over the Pacific Northwest.

A representative time series example at location B for the water year 2015 is shown in Figure 14, which shows the snow evolution from the snow cover assimilation and the joint assimilation configurations along with other model integrations and reference datasets. The snow cover assimilation marginally improves the OL simulation, particularly during January and February, though the univariate albedo DA configuration provides larger improvements. The joint assimilation configuration is consistently better than either of the univariate DA configurations. Even if the changes to the fractional snow cover from snow cover DA is small, it helps to provide a better background for the albedo assimilation and to generate larger improvements in simulated snow depth, as shown in Figure 14.

3.4. Future Enhancements to Albedo DA

There are several possible reasons for the degradations seen in the snow estimates over certain areas of the model domain from the assimilation of albedo. Over such locations, the general trend is a decreased snowpack due to a reduction in albedo. The decreased snowpack occurs as the net radiation available from the updated albedo is large enough to cause an earlier melt. The potential uncertainty in the albedo estimates, particularly under cloudy conditions and over forested regions, could be a factor resulting in spurious albedo updates. Assimilation configurations that allow the consideration of observational errors and uncertainty could be used to partly mitigate these issues. However, the diagnostic nature of the albedo formulation in Noah-MP prevents the use of such DA methods.

As the albedo assimilation schemes rely on changing the net radiation to generate changes in simulated snowfields, the accuracy of the model energy budget formulation is a significant factor in the successful simulation of snowfields. Recent evaluations of reanalysis products [70,71] have established that the uncertainties in the model physics, and not the errors in the driving meteorology, are the key factors in snow simulation errors. The uncertainties in the net radiation impact the success of the albedo assimilation schemes, as demonstrated in the time series comparisons presented in the article. For example, if systematic positive biases in the net radiation are present in the model at a location where snow is underestimated, the reduction of albedo from assimilation will always lead to degraded snow simulations. The model enhancement efforts to improve the energy budgets in Noah-MP should also ensure consistency with respect to the realism of albedo estimates.

The analysis in the article shows that the physical realism of the Noah-MP albedo estimates needs to be improved, particularly in the presence of snow. It must be noted that during the non-snow time periods, the bias differences between the model simulated albedo and assimilated observations are small, as evidenced in the time series examples in Figure 6 and Figure 7. Calibration of the relevant model parameters with reliable measurements of albedo could be used to improve the physical realism of the simulated albedo values. For larger snowpacks, the decrease in reflectance as a result of the increase in the snow grain size and impurities are likely captured by the albedo observations, though processes such as snow aging effects on albedo are not modeled in Noah-MP. These limitations in the snow physics could be a significant source of errors while assimilating the albedo measurements. We expect that the results of this study will serve as a motivation for developing such future improvements in the LSM. Further, this study also functions as a benchmark to evaluate the impacts of such model enhancements within data assimilation configurations.

4. Summary

This article examines the utility of assimilating remotely sensed surface albedo retrievals for improving snow and terrestrial water budget estimates. MODIS-based black-sky and white-sky albedo estimates from the GLASS project are assimilated within the Noah-MP LSM in the NLDAS2 configuration. The impact of assimilation on snow depth, snow cover, and streamflow is systematically evaluated through comparisons against a number of independent reference data products.

Overall, the assimilation of surface albedo retrievals leads to increases in the model simulated albedo in several parts of the domain, particularly over the High Plains and parts of the Rocky Mountains. The increase in albedo generally contributes to reduced availability of net radiation, delayed snowmelt and larger snowpacks. Conversely, smaller snow evolution is noticed over regions where DA leads to a decrease in surface albedo. An analysis of the spectral components of albedo from Noah-MP indicates that there are larger errors in the visible albedo components, which dominate the changes in the reflected radiative fluxes and snow simulations. The impact on runoff and ET from albedo assimilation is found to be small.

Quantitative evaluation of the simulated snow depth against independent estimates from GHCN, SNODAS and UA datasets show that the increase in snowpack largely contributes to improvements in the estimates from DA. Over the High Plains where assimilation leads to larger snowpacks, the RMSE is smaller relative to these reference datasets. Seasonal evaluations indicate that the largest errors are during the melt time periods. It is also notable that there is significant uncertainty across these reference datasets themselves. Though not completely independent, comparisons against the snow cover fraction estimates from MODIS show similar patterns of improvements and degradations in the skill of snow detection. The observed impact of albedo DA on streamflow is marginal, with no statistically significant changes observed over most locations. Some limited patterns of improvements are observed over the upper Missouri and Mississippi basin with degradations over the Northwest. It is notable that these changes are also consistent with the patterns of improvements and degradations observed for snow depth.

As most prior optical sensor-based snow DA studies are focused on the assimilation of fractional snow cover measurements, this study examines the relative (and joint) utility of snow cover DA to albedo assimilation. The snow cover measurements are assimilated using a rule-based direct insertion approach that corrects the mismatches between the model and the observations. Compared to the albedo DA configurations, the changes (both improvements and degradations) in the simulated snowfields are small with snow cover DA. The configuration that jointly assimilates fractional snow cover and albedo data provided improved snow depth fields compared to the albedo only DA configuration. This suggests that the improved snow cover fields from $fsc$ DA is useful in generating additional improvements when assimilating albedo measurements.