Snow cover over land has a significant impact on the surface radiation budget, turbulent energy fluxes to the atmosphere, and local hydrological fluxes. For this reason, inaccuracies in the representation of snow-covered area (SCA) within a land surface model (LSM) can lead to substantial errors in both offline and coupled simulations. Data assimilation algorithms have the potential to address this problem. However, the assimilation of SCA observations is complicated by an information deficit in the observation—SCA indicates only the presence or absence of snow, not snow water equivalent—and by the fact that assimilated SCA observations can introduce inconsistencies with atmospheric forcing data, leading to nonphysical artifacts in the local water balance. In this paper, a novel assimilation algorithm is presented that introduces Moderate Resolution Imaging Spectroradiometer (MODIS) SCA observations to the Noah LSM in global, uncoupled simulations. The algorithm uses observations from up to 72 h ahead of the model simulation to correct against emerging errors in the simulation of snow cover while preserving the local hydrologic balance. This is accomplished by using future snow observations to adjust air temperature and, when necessary, precipitation within the LSM. In global, offline integrations, this new assimilation algorithm provided improved simulation of SCA and snow water equivalent relative to open loop integrations and integrations that used an earlier SCA assimilation algorithm. These improvements, in turn, influenced the simulation of surface water and energy fluxes during the snow season and, in some regions, on into the following spring.
Average monthly snow cover in the Northern Hemisphere varies from 7% of all land area to more than 40%, making snow cover the fastest varying large-scale surface feature on Earth (Chang et al. 1990; Hall 1988). This variability has a dramatic impact on surface moisture and energy fluxes. Snow insulates the ground beneath, moderating soil temperatures during winter (Decker et al. 2003). Because of its high albedo, snow significantly reduces the absorption of radiation at the land surface, restraining turbulent energy fluxes and lowering near-surface air temperature (Baker et al. 1992; Liston 2004). In high-latitude and high-altitude regions, snow plays a dominant role in the local and regional hydrologic balance, with associated impacts on soil moisture, aquifer recharge, and vegetation growth (Flerchinger et al. 1994; Marsh 1999; Ren et al. 2007). Moreover, it is now understood that snow’s influence on local energy fluxes has remote impacts, modifying atmospheric temperature and circulations on the regional scale (Elguindi et al. 2005; Ellis and Leathers 1999; Zaitchik et al. 2007) and climate dynamics on the continental to hemispheric scales (Bamzai and Shukla 1999; Cohen and Entekhabi 2001; Dery et al. 2005).
Given these significant effects, it is critical that models used in seasonal forecasts and retrospective climate studies accurately simulate the snowpack. For this reason, most modern land surface models (LSMs) include sophisticated snow routines designed to capture the accumulation, aging, and melt of snow under a range of weather conditions (Ek et al. 2003; Roesch et al. 2001; Stieglitz et al. 2001; Takata et al. 2003). Unfortunately, even the most refined LSMs have difficulty capturing these processes over an extended simulation (Foster et al. 1996; Pitman 2003). In a sense, snow is an inherently difficult quantity to model: for large areas of the planet, accumulation and melt processes are sensitive to small changes in air temperature, radiation fluxes, or surface properties, yet the presence versus absence of snow will itself have a dramatic impact on these variables. This means that small uncertainties in parameterization or forcing data can trigger errors in snow cover that quickly cascade into larger inaccuracies in albedo, soil moisture, energy exchange, and atmospheric conditions. As these inaccuracies compound, the model “drifts,” leading to severe reductions in skill on the seasonal time scale. For hydrologic applications, drift with respect to snow simulation is an even larger problem; in many regions, accurate simulation of the snow water equivalent (SWE) of the snowpack and timing of snowmelt are the most important aspects of hydrological monitoring and prediction systems.
The Moderate Resolution Imaging Spectroradiometer (MODIS) instruments, on the National Aeronautics and Space Administration (NASA) Terra and Aqua satellites, offer twice-daily observations of the snow-covered area (SCA) at 500-m resolution, with near-global coverage (Hall et al. 2002). The reliability of these observations is much greater than that of LSM predictions at a similar resolution, but the information contained in a MODIS observation is much less. The LSM simulation of the snowpack includes information on SWE, snow depth, albedo, and, in some models, time-variable snow density. The simulation is also continuous in space and time, at a user-specified resolution. MODIS contains no information on SWE or snow depth, can provide a maximum of two daylight observations per day, and frequently suffers from data gaps due to cloud cover. Because the MODIS SCA product is derived from information in the visible and near-infrared portions of the electromagnetic spectrum, observations can only be obtained under clear-sky, daytime conditions.
Data assimilation (DA) algorithms attempt to use the information in such discontinuous observations by integrating them into a numerical model. The model provides spatial and temporal continuity as well as physically based schemes that allow an observation of one quantity to inform predictions of other modeled variables. The observation system provides reliable, independent information on a variable that the model simulates imperfectly, preventing model drift and improving the accuracy of the simulation. The MODIS SCA observation is a strong candidate for use in a DA system because it has demonstrated reliability (Hall et al. 2002), but it is discontinuous in time and provides hydrologically incomplete information. Further, because LSMs are frequently employed in a coupled mode with atmospheric models, a MODIS SCA assimilation system can be used to improve the initialization of numerical forecasts, a powerful predictive application of MODIS-derived information.
In this study, MODIS observations of SCA were assimilated into the Noah LSM (Chen et al. 1996; Ek et al. 2003) in global, retrospective simulations at 0.25° resolution. Simulations were performed for two Northern Hemisphere snow seasons, September 2005–June 2006 and September 2006–June 2007. The results of two assimilation algorithms are presented. The first, a rule-based direct insertion scheme, is a refinement of a system first presented by Rodell and Houser (2004, hereafter RH04). In this algorithm, MODIS SCA observations are assigned a SWE equivalent on the basis of land cover, and this value is used to overwrite simulated SWE whenever there is a contradiction between model and observation. The second algorithm is a novel, physically based assimilation update that uses MODIS SCA observations from up to 72 h ahead of the simulation to adjust snow accumulation and melt processes. This allows for smooth updates of snow that attempt to preserve the local hydrologic balance while preventing unrealistic add/melt and remove/accumulate cycles that occur in standard snow assimilation schemes whenever observations are in conflict with atmospheric forcing data.
The remainder of the paper is structured as follows: section 2 provides background on the Noah LSM, satellite-derived snow observations, and data assimilation. Section 3 gives details of image processing, validation data, and the assimilation scheme. Simulation results and evaluation are presented in section 4. Section 5 contains conclusions.
a. Noah LSM
The Noah land surface model is an advanced 1D column model that simulates soil temperature, surface temperature, soil moisture (frozen and liquid, in four soil layers), snow states, surface energy and hydrologic fluxes, and subsurface runoff (Chen et al. 1996). The model is currently supported by the National Centers for Environmental Prediction (NCEP) and is in active development (e.g., Ek et al. 2003). It is widely used in operational and research applications and can be run in either the coupled or the uncoupled mode. Noah employs finite-difference spatial discretization methods and a Crank–Nicholson time-integration scheme to solve the governing equations of the physical processes of the soil-vegetation-snowpack medium.
Noah contains a one-layer snow model that simulates snow water equivalent (SWE) as the residual of snowfall (snowmelt + sublimation). Snowfall occurs whenever there is nonzero precipitation and 2-m air temperature is ≤ 0°C. Melt occurs whenever there is net available energy within the snowpack. When melt occurs, some water is stored within the snowpack and excess water is removed as runoff. Fractional snow-covered area is diagnosed as a function of SWE using a generalized snow depletion curve:
where τ is the snow distribution shape parameter that relates the total amount of snow in the grid cell to the percent snow cover within the grid cell (currently set to 4.0, globally) and SWESCA=1.0 is a land cover specific parameter that defines the minimum SWE required for full snow cover (ranging from 13 mm for bare soil to 40 mm for forests). Snow depth is solved on the basis of SWE and snow density. Snow density is itself a predicted field and is a function of total SWE in the snowpack, a temperature and SWE-dependent compaction scheme, and air temperature for freshly falling snow (Koren et al. 1999).
b. Satellite-derived snow observations
Satellite-based snow detection algorithms can be divided into two general categories: those that attempt a direct measurement of SWE and those that produce estimates of SCA. SWE algorithms take advantage of the fact that the diameter of snow grains falls within the same size range as microwave radiation (Chang et al. 1987, 2005). This means that snow scatters microwaves of different wavelengths with widely differing efficiency and that the difference in scattering between two or more microwave channels can be used to estimate SWE. Algorithms founded on this principle have been used to retrieve SWE from the Special Sensor Microwave Imager (SSM/I), the Advanced Microwave Scanning Radiometer for Earth Observing System (AMSR-E), and other passive microwave sensors. All such algorithms, however, are plagued by the fact that microwave scattering is sensitive to snow grain size, phase changes within the snowpack, vegetation, ice lenses, and the presence of liquid water on the snow surface. While there has been considerable progress on the theory of SWE retrieval from passive microwave sensors (e.g., Durand and Margulis 2006), these data products are not yet ready for practical use in DA (Andreadis and Lettenmaier 2006), particularly in global applications (Cordisco et al. 2006).
Detecting the presence or absence of snow is a simpler problem for remote sensing, and several sensors currently offer acceptably accurate estimates of SCA over the globe (Bitner et al. 2002; Brubaker et al. 2005). With respect to the two general goals of assimilating observations of snow—to improve representation of (i) the surface energy balance and (ii) snow water reserves—the assimilation of estimated SCA offers the potential to improve (i) in all regions of the globe and (ii) only in locations with partial or sporadic snow cover, or areas where there is substantial model error in the timing of the initiation and final melt of the seasonal snowpack. This expectation follows from analyses of in situ snow observing stations, which show strong correlations between SCA and SWE during periods of accumulation and ablation but weak correlations in the middle of the snow season, when SCA nears 100% in a snowy region (Gong et al. 2007).
SCA can be retrieved by passive microwave sensors, such as SSMI and AMSR-E, that provide measurements under both clear and cloudy conditions but have coarse spatial resolution (25 km or more) and by sensors that use the visible and near-infrared (VNIR) portion of the spectrum. VNIR sensors cannot penetrate clouds, but they offer superior spatial resolution. MODIS, for example, can detect snow presence/absence at 500-m resolution, which is more than adequate for most climate and many hydrological applications. MODIS was selected for this study on account of its high spatial resolution and general reliability. It is our expectation that DA routines can effectively interpolate between discontinuous measurements in time and space but that it is more difficult to disaggregate a spatially coarse measurement or to extract useful information from an error-prone observation. Because our goal is to design an assimilation scheme that can be applied globally but also at resolutions finer than 25 km, MODIS was identified as the best available data source.
c. Data assimilation
The guiding principle of DA is that an observation and a prediction can be combined such that the result contains less uncertainty than either input. Here, the predictions are provided by a physically based numerical model of land surface processes. The model contributes a high-resolution, spatially and temporally continuous, first-guess approximation of the assimilated state (in this case, snow water equivalent), based on numerical representations of physical processes and other relevant data. The observations provide an external check on model predictions, preventing drift and correcting the model with respect to essential state variables. Hence, the data assimilation scheme enables both scientific understanding and observations to inform estimates of a given variable, as well as associated states and fluxes, while addressing the spatial and temporal limitations of the observing system.
The assimilation of a SCA observation presents some challenge because it provides only partial information on the model field that it is used to inform. A SCA observation contains no information on snow depth, snow density, or snow water equivalent, which will all be influenced by the assimilation update, and the physical link between SCA and these variables is inconsistent in space and time. An early season snowfall, for example, might coat 100% of the observation footprint with a very thin layer of snow, while springtime melt might leave behind deep piles of dense snow that cover only 30% of the footprint.
The relationship between SCA and SWE at a given location can be described using a contextual snow depletion curve (SDC), which relates SCA to SWE on the basis of prior information and ground measurements (e.g., Anderson 1973). Ideally, such a geographically specific SDC is informed by local data on typical patterns of snow-depth distribution (e.g., Andreadis and Lettenmaier 2006; Kolberg et al. 2006). Because such data are not available globally, Liston (2004) proposed a dichotomous key that divides the world into nine SDC landscapes, allowing for some refinement to globally uniform SDCs in the absence of local data. The key lacks the local specificity and spatial resolution of some regionally specific SDCs, but it was shown to have a significant impact on the simulation of snow processes in a regional model. The Subgrid Snow Distribution (SSNOWD) model derived from this key could be implemented within the framework of many existing LSMs.
Another approach to SCA assimilation is to map SCA observations onto predicted variables in the LSM on the basis of cross-state covariance estimates. Andreadis and Lettenmaier (2006) applied an ensemble Kalman filter (EnKF) assimilation scheme that used SCA observations to directly update SWE, while Clark et al. (2006) demonstrated a more extensive EnKF scheme that mapped SCA observations to update SWE, soil moisture, and groundwater. Both studies showed promise in a watershed scale, but the requirement for reliable cross-state covariance estimates and statistically robust observation error estimates currently preclude global application of the techniques (Clark et al. 2006).
Acknowledging the information deficit encountered by both regionally specific SDC assimilation schemes and EnKF routines, RH04 proposed a simple, conservative rule-based assimilation scheme to link SCA observations to SWE. In their assimilation of MODIS SCA into the Mosaic LSM, they implemented a rule in which a thin layer of SWE (5 mm) was applied to any model grid cell for which MODIS reported SCA was >0.4 and the LSM had no snow cover. For cells with MODIS SCA <0.1, all snow was removed from the grid cell. This simple approach can be applied globally with reasonable confidence, and it requires no modification to the basic routines of the LSM. As our goal is a globally applicable assimilation routine that can be applied to existing operational versions of the Noah LSM, we use the RH04 algorithm as the starting point for the algorithm developed in this study.
The RH04 algorithm, like most SCA assimilation methods, suffers from the inherent sensitivity of a snow update to any inconsistencies that might exist between the SCA observation and forcing fields of air temperature and precipitation. This leads to two problems. First, SCA assimilation tends to be much more effective at removing snow than adding snow. This is because model errors that increase SCA are the product of episodic snowfall events, and these events can be quickly counteracted by a single snow-free observation. Snowmelt, on the other hand, can occur continuously as a function of air temperature. The assimilation of a SCA observation that adds SWE to a grid cell is subject to immediate melt if air temperature is above the freezing point. If added snow cover is melted away before the next observation becomes available, then a primary goal of SCA assimilation—to correct SCA and surface albedo, preventing errors in surface energy fluxes to the atmosphere—is compromised, and, if it happens day after day (an “add/melt cycle”), soil moisture and runoff are unduly amplified. In RH04 this problem was observed at a few isolated points, but on the whole assimilation was effective, in large part because the LSM used in that study, Mosaic, suffered from overaccumulation rather than underaccumulation of snow. For an LSM such as Noah, which is prone to underestimating snow, the problem of preventing inaccurate snowmelt is more critical.
The second related difficulty in updating snow fields in an LSM is that most assimilation schemes upset the water budget without consideration of the underlying errors in the simulated precipitation and melt processes. Such introduced imbalances can compound over time and manifest themselves as soil moisture and runoff biases. As mentioned above, the advanced DA schemes capable of simultaneously updating all snow-related hydrological states at once (e.g., Clark et al. 2006) are not yet ready for global applications.
In an effort to mitigate hydrologic imbalance caused by simply adding or removing snow, we present a novel assimilation algorithm that uses MODIS observations from up to 72 h ahead of simulation time. These “future” data are used to nudge the forcing air temperature and, when necessary, precipitation toward likely precursors of the observed SCA. When it is successful, this approach minimizes disruption of the local water balance and provides a smooth simulation of snowpack as informed by MODIS observations. We feel that the use of SCA observations to nudge forcing data is justified by the reliability of MODIS observations. Furthermore, the actual impact on forcing fields typically is small: snow accumulation and melt are sensitive to small changes in air temperature in the range of 0°C, so it is rarely necessary to alter forcing air temperature by more than a few degrees Celsius. In a coupled simulation, the modification would not directly affect the atmospheric model, as air temperature within the LSM is not communicated back to the atmosphere. The modification of air temperature could indirectly affect atmospheric processes, via its influence on surface energy fluxes, but that influence was found to be modest relative to the energetic impacts of snow cover itself (see section 4).
In concept, this algorithm bears some resemblance to the EnKF update presented by Clark et al. (2006), in which snow accumulation and melt were used to update SWE indirectly on the basis of a SCA observation, and soil moisture and groundwater were updated statistically on the basis of locally determined coefficients of variation. Beyond this the methods diverge, as our scheme employs model physics rather than statistical relationships to update hydrological fields. This avoids the need for detailed local cross-state covariance estimates, allowing for global applications.
a. MODIS data
In this study we use the daily, 0.05°-resolution MODIS climate-modeling grid-level-3 product (MOD10C1), which is based on 500-m Terra MODIS observations (Hall et al. 2002). MOD10C1 specifies the fraction of each 0.05° grid cell that was observed to be snow covered, the fraction that was cloud covered, and the fraction (known as the “confidence index”) in which the land surface was visible (i.e., not obscured by clouds, night, or other interference) at the time of the satellite overpass [approximately 1000 local time (LT)]. The MOD10C1 product has been evaluated extensively against independent satellite and ground-based datasets, and it has shown good agreement (Bitner et al. 2002; Maurer et al. 2003). Nonetheless, there are acknowledged limitations to MOD10C1, including the binary nature of the 500-m observations that are used to derive percent cover and the difficulty of distinguishing patchy snow from clouds in some regions.
Following RH04, the reliability of a MODIS snow-cover observation within a given grid square was determined using the MOD10C1 confidence index. MODIS is not cloud penetrating, but clouds are often pervasive where snow exists, so it is essential to make prudent use of data from grid squares that are partially cloud covered, lest useful information be neglected. Taking this reasoning into account and based on a visual assessment of the credibility of the observed snow-cover state at varying levels of the confidence index, it was decided that 6% is the minimum visibility for which a 0.25° aggregation of observations is useful. This parameter can easily be adjusted if a more appropriate value is later identified. For 0.25° tiles that achieved this confidence threshold, the percentage of visible, snow-covered 500-m pixels is divided by the confidence index to establish the fraction of ground-visible (cloud free) pixels that were snow covered at the time of the MODIS observation.
We selected the relatively coarse spatial resolution of 0.25° to accomplish lengthy global simulations, but MODIS would allow the application of the same algorithm at significantly higher resolution. Recently, Salomonson and Appel (2004) have extracted fractional SCA from MODIS at 500-m resolution. High-resolution assimilation experiments using this product are the subject of current investigation.
b. Assimilation procedure
We implemented Noah LSM (version 2.7.1) for global, uncoupled simulations at 0.25° resolution for the periods September 2005–June 2006 and September 2006–June 2007. Experiments were initialized using Global Land Data Assimilation System (GLDAS) Noah restart fields from a multidecadal simulation parameterized and forced by a hybrid dataset (Rodell et al. 2004). Three-hourly atmospheric analyses from the Global Data Assimilation System (GDAS), the operational atmospheric DA system of NCEP (Derber et al. 1991), forced the experimental simulations. GDAS runs on a thinned Gaussian grid, with a resolution of about 0.47°. The operational 2.5° 5-day Climate Prediction Center (CPC) Merged Analysis of Precipitation (CMAP) was used to bias correct the higher-resolution GDAS precipitation fields (Rodell et al. 2004). CMAP estimates are based on a blending of satellite data (microwave and IR) and gauge observations (Xie and Arkin 1997). The combination of GDAS forcing with CMAP bias correction for precipitation was selected to match the standard set of forcing fields that drive GLDAS, as the DA algorithms developed in this study will eventually be used in GLDAS.
All Noah simulations were performed using the NASA Goddard Space Flight Center (GSFC) Land Information System (LIS), version 5.0. LIS is a software framework that supports the integration of observational data with coupled or uncoupled LSMs using a range of DA techniques (Kumar et al. 2006). Implementation of the assimilation algorithm within LIS allows for simulations with a number of different atmospheric forcing datasets and for application to coupled land–atmosphere modeling systems. The LIS framework also facilitates transfer of the assimilation scheme to forthcoming versions of Noah or to other similar LSMs. In LIS–Noah, subgrid variability is captured by simulating each major land-cover class within the grid cell as an independent, 1D “tile.” Gridded output is returned as the area-weighted sum of tiles within a grid cell (Rodell et al. 2004).
As discussed in section 1, we used two assimilation procedures in this study, the push and pull algorithms (Fig. 1). The push algorithm consists of a single update routine, invoked once daily for each LSM tile. It is largely analogous to the scheme used in RH04, slightly modified for Noah’s snow routines.
1) The push algorithm
At 1000 LT, Noah SCA is checked against the present MODIS observation. At this time, if Noah differs from MODIS (LSMSCA ≠ OBSSCA), then
(i) if MODIS SCA > 0.5 and Noah SCA < 0.2, Eq. (1) is inverted to calculate the minimum depth of SWE that must be added to the tile to bring Noah into agreement with MODIS; this quantity is added directly to SWE, and snow depth is also updated, on the basis of added SWE and Noah’s snow density routine for incoming snowfall;
(ii) for MODIS SCA < 0.2 and Noah SCA > 0.5, Eq. (1) is inverted and Noah SWE is replaced with the value indicated by MODIS.
These thresholds were chosen based on the properties of the MODIS observation. The MODIS SCA product reports a “bird’s-eye view” of snow cover rather than an estimate of actual snow cover on the ground (as is produced by most microwave-based snow-cover products). This means that in areas with tall vegetation or steep topography, MODIS frequently reports SCA below 100%, even when the ground below canopy is completely snow covered. For this reason we take a conservative approach, removing snow only when MODIS SCA is particularly low (<20%) and, in adding snow, only adding to nearly snow-free tiles for which MODIS reports considerable (>50%) snow cover. For all other tiles, we accept that LSMSCA ≈ OBSSCA and do not update the model. The push algorithm differs from that of RH04 in that the volume of SWE added/removed is determined by inverting the internal Noah snow-depletion curve rather than setting it as a constant for the entire simulation. Use of the SDC allows for land-use specificity in the assimilation update and allows us to apply the update as a continuous function of SCA. The Noah SDC is highly generalized, however, and is in no way equivalent to the detailed, geographically specific SDCs applied in watershed-scale studies (e.g., Andreadis and Lettenmaier 2006; Kolberg and Gottschalk 2006).
2) The pull algorithm
The pull algorithm comprises two update routines, the first invoked for every tile at each LSM time step, and the second invoked only once daily for each tile.
At each time step, the assimilation scheme checks the simulated SCA for the tile [calculated using Eq. (1)] against the most proximal future MODIS SCA observation. This is the same-day MODIS observation when local time is 10:00 A.M. or earlier and it is the next-day MODIS observation when local time is later than 1000 LT. When the most proximal future observation is unavailable, the algorithm looks up to 3 days (72 h) ahead to find a usable observation. If current Noah SCA differs from the future MODIS observation, then forcing is adjusted as follows.
(i) If MODIS SCA > 0.5 and Noah SCA < 0.2, forcing air temperature is set to −3°C for the time step—this inhibits melt and causes any incoming precipitation to be treated as snowfall.
(ii) If MODIS SCA < 0.2 and Noah SCA > 0.5, forcing air temperature is set to +3°C—this promotes melt and prevents the accumulation of any new snowfall.
Whenever local time is equal to 0945 LT—that is, the time step just prior to the time of MODIS observation—the precipitation field is modified directly.
(i) For MODIS SCA > 0.5 and Noah SCA < 0.2, if precipitation in the forcing dataset is insufficient to produce the MODIS-observed SCA by the end of the time step, the precipitation rate is increased to the minimum rate required to produce SCA that is consistent with MODIS. Routine (1) guarantees that this added precipitation will fall as snow.
(ii) For MODIS SCA < 0.2 and Noah SCA > 0.5, the remaining snowpack is forcibly melted, with meltwater converted to soil moisture and, if the soil saturates, to surface runoff.
All modifications to the forcing fields are recorded throughout the simulation to track the impact of assimilation on simulated air temperature and precipitation. Modifications to forcing imposed through assimilation are relevant only to the LSM, not to any coupled atmospheric model; in coupled simulations, precipitation and air temperature are inputs to the LSM, which are not reported back to the atmosphere. The modifications will, however, influence surface fluxes of energy and water. Effects on these variables are presented in section 4.
Forcible melt of a recalcitrant snowpack [pull, step 2(ii)] was included in the algorithm because Noah determines the available energy for snowmelt as a function of both air temperature and soil surface temperature (Koren et al. 1999). In locations where the LSM has erroneously generated a deep snowpack, the insulating properties of the snow inhibit melt from the top. This difficulty in melting a thick snowpack from the top is inherent to the one-layer snow model in Noah 2.7.1, and would be of less concern for LSMs that contain multiple-layer snow models. In practice, however, the last rule 2(ii) was rarely invoked in the global 0.25° Noah simulations reported in this paper.
c. Evaluation data
Model predictions of SCA and SWE were evaluated using ground-based observations wherever possible. For the United States, snow-depth observations were provided by the National Weather Service Cooperative observer program (Co-op). The Co-op network does not provide direct measurements of SWE, but it is nonetheless preferable to other snow observation networks (e.g., SNOTEL), both because it offers extensive coverage over the United States and because it is not specifically designed as a snow observation network. Observation stations designed for snow monitoring are intentionally sited in high-altitude locations with deep-snow accumulation, as these are the areas of greatest relevance for most snow studies. Because SCA updating is most effective in areas that have marginal snow cover or an ephemeral seasonal snowpack, Co-op stations, which are distributed randomly with respect to snow-accumulation zones, are a better source of evaluation data. Snow water equivalent was estimated as 8.0% of observed depth at Co-op stations, representing an average for the dynamic range of snow density in the Noah LSM.
For locations outside of the United States, observations of snow depth were extracted from the National Climatic Data Center (NCDC) Global Summary of the Day (GSOD) dataset. As with Co-op stations, GSOD observations provide information on snow presence/absence and snow depth. Unfortunately, many GSOD stations report snow observations only when snow depth is greater than zero. It is still possible to evaluate the distribution of snow cover predicted in Noah simulations against the distribution reported by GSOD stations, but with the understanding that snow-free station locations are underrepresented.
In all cases, our analysis is limited by the difficulty of using point location in situ data to evaluate a relatively coarse-scale gridded model. Wherever possible, multiple observations within each 0.25° grid cell were averaged to provide more representative evaluation data for model results. In the majority of cases, however, only a single observation was available for each grid cell, so this single data point was assumed to be representative of the 0.25° area. For assessment of snow presence/absence, snow was deemed to be “present” whenever reported depth exceeded the snow depth required for 50% grid cell coverage in Noah LSM.
4. Results and discussion
a. SCA and SWE
Over most regions of marginal snow cover, data assimilation using the push or pull algorithm led to an increase in total snow-covered area relative to open-loop simulations (Fig. 2). This included earlier onset of the snow season in northern Europe and across most of North America, larger total extent of winter snow cover in the western United States and central Asia—but with markedly less winter snow in Tibet—and later spring melt for much of Canada and, with the exception of Tibet, Eurasia.
Comparisons with in situ reports of snow presence/absence indicate that DA-informed simulations of snow cover were substantially more accurate than the open loop simulation. First, open-loop simulations with Noah systematically underestimated snow presence at reporting Co-op stations in the contiguous United States (CONUS; Fig. 3a). This was true for test regions outside of the United States as well, with the largest discrepancy occurring during springtime melt in high-latitude regions (central Canada and Siberia, Figs. 3b,c) and throughout the snow season in midlatitude regions (Mongolia, Fig. 3d; CONUS). In central Canada, the springtime (March–May) hit rate for observed snow cover was 68% for open loop versus 86% in pull in 2005–06 (n = 1071) and 65% open loop versus 83% pull in 2006–07 (n = 975). In Siberia, 2005–06 springtime hit rates were 75% and 88% (n = 2484) and 2006–07 hit rates were 74% and 89% (n = 1889) in open loop and pull, respectively. Summed over the entire snow season, the ratio of simulated snow cover to observed snow cover (ρ, Table 1) indicates that both push and, more dramatically, pull simulations corrected the low snow bias for the CONUS, Siberia, and Mongolia during both the 2005–06 and 2006–07 snow seasons. In central Canada, only pull led to a reduction in the low snow bias. Reduction in bias was associated with an improvement in overall accuracy for all tested regions in both years, with the exception of central Canada in the push simulation (Table 1). This improvement can be quantified using Cohen’s Kappa statistic, which accounts for errors of omission and commission in the assessment of accuracy (Congalton and Green 1999). Kappa analysis indicates that improvement over open loop was statistically significant (α <0.05) for both push and pull as measured against Co-op stations for the CONUS in 2005–06 and 2006–07. In both cases, both assimilation runs were significantly improved relative to open loop, and pull performed significantly better than push. Outside of the CONUS statistical measures of improvement were weaker, largely due to underreporting of snow-free conditions at GSOD stations (see section 3). Nevertheless, statistically significant improvement was observed for the pull simulation in Mongolia and central Canada in 2005–06. In Siberia, and in all three non-U.S. regions in 2006–07, improvements were not statistically significant in this analysis.
Unfortunately, station reports from Tibet were too limited to allow for a full station-based evaluation of DA results in this interesting region. However, the general tendency of DA to produce patchy, rather than complete, snow cover for Tibet is consistent with the known character of snow cover on the Tibetan plateau, as established in field survey and independent remote sensing analyses (Qin et al. 2006). MODIS SCA observations, for their part, have been validated for Tibet in a focused regional-scale study (Pu et al. 2007), further increasing our confidence in DA results for the plateau.
With respect to snow volume, DA had a positive influence on the simulation of SWE in regions with partial or transient snow cover, including the U.S. Southwest, where deep-snow accumulation occurs only in elevated areas and snow cover exhibits strong interannual variability, the U.S. high plains, and the Mongolian steppe (Fig. 4). Following the pattern observed for SCA, DA tended to increase SWE relative to open loop simulations in these regions. In regions with deep seasonal snowpack, assimilation was only able to provide a marginal improvement in the simulation of SWE. Such was the case for the West Coast focus region, which included the Sierra Nevada and Cascade ranges. DA increased SCA in these ranges (see Fig. 2) and did provide a nominally statistically significant improvement in SWE for certain months (not shown). However, DA could not improve model simulation of SWE in a hydrologically meaningful way for the deep-snow conditions observed in 2005–06 (Fig. 4), as MODIS provides no information on snow mass or depth.
b. Hydrological fluxes
As the SCA and SWE results indicate, the pull approach applied MODIS information more effectively than the push approach, primarily because it offset discrepancies between the satellite observation and atmospheric forcing data. By nudging the atmospheric forcing into consistency with the MODIS observation, the algorithm updated snow fields in a more robust manner, reducing the incidence of daily add/melt cycles and related artifacts. Thus pull maintained a realistic local hydrological balance throughout the simulation. In contrast, the simple addition or removal of snow in the push simulation overwhelmed the water balance in certain locations. This is demonstrated in Fig. 5. For the CONUS, the sum of daily residuals in the local water balance in pull closely resembled those of the open loop simulation: in both cases there is some variability around zero, and pull exhibits a small positive bias due to deep-snow corrections early in the simulation and, to a lesser degree, during the melt season. In contrast, push yields consistent negative departures from balance, and these become highly significant in late winter and spring. This is inconsistent with model physics and is a product of the repeated add/melt cycles in the push simulation.
Figure 6 shows an example of the contrasts between push and pull for a location in Colorado where MODIS observed nearly complete cover throughout the winter months, while open loop Noah simulations exhibited sporadic coverage (Fig. 6a). Both the push and pull integrations effectively corrected the low snow bias, adding SWE to the simulation, and thus bringing the model into better agreement with Co-op observations of snow depth (Fig. 6b). Push, however, accomplished this by the direct addition (removal) of snow to (from) the land surface, which caused nonphysical imbalances in the local water balance (Fig. 6c). Pull maintained the local hydrological balance as effectively as did open-loop simulations by anticipating the presence of snow and reducing air temperature (Fig. 6g). The relatively modest reduction in air temperature was sufficient to change the phase of available incoming precipitation to snow (Fig. 6e) from rain (Fig. 6f), and additional snowfall was added where precipitation was inadequate to provide snow cover. The physically consistent approach to snow updates also had an influence on snowmelt (Fig. 6h) and, by extension, soil moisture. In this case, runoff (Fig. 6d) was small for all simulations. This example also demonstrates some of the limitations of MODIS SCA assimilation: neither push nor pull can remedy the model’s low snow bias in response to the precipitation event on 20 December, as the model achieved 100% snow cover, and SCA observations cannot provide any further information. Also, the accumulated snowpack melted from DA assimilations during a gap in MODIS data in February. Because there is a corresponding gap in Co-op data, we don’t know whether this melt was accurate, but it was unavoidable due to the MODIS data gap.
Figure 7 demonstrates the advantages of pull for a location where DA had to remove snow from the model. In this example, from a location in Wisconsin, the open loop Noah simulation erroneously overestimates the magnitude of a late-November snow event, creating a snow cover error that persists for the first half of December. In the push simulation, low-snow MODIS observations allow for a rapid correction of the erroneous snow accumulation, but this is not achieved until after the initial error has occurred (Figs. 7a,b), and the correction imposes a nonphysical departure from local water balance (Fig. 7c). In pull, MODIS advance observations were utilized to prevent false snow accumulation from occurring. This was accomplished through a slight increase in air temperature (Fig. 7g) that caused precipitation to fall as rain rather than snow (Figs. 7e,f). As a result, the pull simulation was closer to Co-op observations (Fig. 7b) and exhibited less snowmelt (Fig. 7h) and runoff (Fig. 7d) over this period.
The secondary pull add-snow procedure—increasing the local precipitation rate when air temperature modifications failed to produce an adequate snow cover—preserves the surface water balance at the expense of a possible precipitation bias. That consequence should be examined closely in any application. For the simulations performed in this study, the influence of pull on precipitation was substantial in some locations, but differences between pull and open loop precipitation were generally within the uncertainty level, that is, the range of estimates from various global precipitation datasets. In fact, the pull algorithm often brought the open loop precipitation (disaggregated CMAP) into better agreement with alternative forcing datasets (Table 2). One exception is for high-latitude regions in the spring. In these regions, the DA-induced addition of precipitation due to the pull algorithm did sometimes impose a wet bias beyond the range of common precipitation datasets. This assimilation artifact reflects a limitation in the MODIS SCA data. Because MODIS depends on visible and near infrared data to determine snow cover, the sensor cannot provide SCA estimates at high latitudes during winter. During this extended data gap, the model has the opportunity to drift to lower SCA values (e.g., Fig. 3, central Canada). When the MODIS record for these regions resumes in spring, a relatively large volume of snow must be added to the simulation to regain accurate snow cover. To address this issue, we have begun to investigate the use of multisensor snow-cover products (Helfrich et al. 2007; Romanov et al. 2003) that use high-resolution visible/near-infrared observations whenever available but use microwave data to fill in data gaps at high latitudes. It is anticipated that such multisensor products will correct for any artificial wet bias that pull produces at high latitudes.
Pull generally produced more realistic snowmelt, soil moisture, and runoff than push. This was most obvious in locations where push repeatedly caused add/melt cycles over the course of the snow season (e.g., Fig. 7). These cycles resulted from inconsistencies between observation (SCA > 0.5) and atmospheric forcing (Tair > 0°C). By making small adjustments to air temperature (Fig. 7g), pull made it possible for the LSM to retain an applied positive snow increment rather than melt it away (Figs. 7b,c,h). Even under less dramatic circumstances, differences between pull and push can accumulate over the season to have a substantial impact on total snowmelt, as well as the associated fields of soil moisture and runoff.
c. Energy fluxes
As expected, data assimilation had a substantial impact on the surface energy balance and on land–atmosphere energy fluxes. Changes to SCA altered surface albedo, which in turn influenced radiative and turbulent fluxes. These effects were most pronounced at midlatitudes during winter and at high latitudes in fall and spring, and they were quite substantial for some regions. On the West Coast of the United States, for example, DA increased snow- covered area and albedo, causing a significant reduction in net radiation and sensible heat flux in winter (Table 3). The effects of the SCA update persisted into spring, when snowmelt in the DA simulations led to significantly enhanced latent heat flux, along with a nonsignificant increase in soil moisture and decrease in surface temperature. The influence of DA on energy fluxes was also observed in midlatitude Eurasia. In Mongolia, DA produced simulations with increased wintertime SCA and albedo and significantly decreased sensible heat flux and net radiation (net radiation effect significant for pull only). The persistence of snow memory in the form of soil moisture was smaller in this region than it was on the West Coast.
At high latitudes, where winter snow cover was not ephemeral, DA had a smaller impact on energy fluxes during winter, but it had consistent, significant impacts in spring (Table 3: central Canada and Siberia). Both push and pull increased springtime snow cover and surface albedo in these regions, which in turn reduced net shortwave radiation at the surface (statistically significant for pull only). Outgoing longwave radiation also decreased, owing to reduced surface temperature. The reduction in surface temperature was caused by increased snow cover and, in the case of the pull simulation, reduced air temperature. Sensible heat flux was diminished significantly relative to open loop because of the combined effects of reduced net surface radiation and enhanced latent heat flux.
These results demonstrate the value of SCA observations for coupled land–atmosphere simulations. They also allay a concern associated with the pull algorithm— that manipulation of air temperature might lead to unrealistic radiation and energy fluxes. For example, in locations where the pull algorithm reduces air temperature to increase SCA, it is possible that the enhanced surface-to-air temperature gradient might lead to increased surface sensible heat flux, even though the addition of snow would be expected to reduce it. In this application, it was found that pull imposed only small changes in the air temperature forcing field, and any influence that this had on turbulent fluxes was small relative to the influence of improving the model’s representation of SCA. This is clear from the fact that push and pull had similar influences on energy fluxes in all regions: increased SCA led to substantial reductions in sensible heat flux, and any mitigation of this due to enhanced surface-to-air temperature gradient in pull was negligible both in the monthly average (as shown in Table 3) and on daily and 3-hourly time scales. Similarly, the influence of DA on radiative fluxes was dominated by the impact of SCA on surface properties rather than the influence of pull on atmospheric forcing fields. Even in high-latitude regions such as central Canada, where pull reduced surface temperatures and net longwave radiation relative to push, the energetic differences between simulations is attributed to the fact that pull more effectively maintained snow cover, leading to larger values of SCA. The direct manipulation of air temperature (Table 2) was small relative to simulated differences in surface temperature (Table 3) and would be expected to have the opposite influence on longwave radiation than was observed in the simulation: the magnitude of net outgoing longwave radiation at the surface was reduced in pull relative to other simulations, even though the manipulation of air temperature decreased incoming radiation, which would, on its own, have increased outgoing net radiation.
The assimilation of MODIS SCA observations to the Noah LSM yielded improved estimates of snow cover in offline global simulations. This had significant impacts on the surface radiation balance, indicating that the technique would be of value in coupled land–atmosphere simulations. Regions where assimilation had a significant impact include the western United States and Mongolia, in winter, and Siberia and central Canada in spring. Both assimilation schemes tested in this paper—the rule-based push assimilation scheme and the pull assimilation scheme that updated air temperature and precipitation fields—improved the simulation of SCA and SWE relative to open loop simulations. Pull tended to produce more robust results, in large part because positive snow increments in push often melted away rapidly. Any unintended influence that pull might have had on surface energy fluxes as a by-product of air temperature manipulations was small compared to the effect of updating SCA. Improved realism in the simulation of SCA and associated surface energy fluxes is important for retrospective climate analyses and, in the context of coupled models, for seasonal forecasts. This is particularly true for regions like Siberia, where variability in SCA is known to have a significant impact on atmospheric circulations and climate throughout the Northern Hemisphere (Gong et al. 2003), and for midlatitude regions such as the U.S. Great Plains, where snow influences springtime soil moisture and vegetation, both of which are involved in important land–atmosphere feedbacks on climate (Koster et al. 2004). The degree to which local improvements in the simulation of SCA lead to improved seasonal predictions on the regional scale will, of course, depend on the accuracy of model physics and parameterizations (Schlosser and Mocko 2003).
With respect to hydrological variables, the benefits of data assimilation were more modest. The assimilation of SCA is inherently limited in this regard, as most snow is stored in regions of deep snow, where SCA observations have little information to provide beyond the first few weeks and last few weeks of the snow season. In evaluation against snow-depth measurements from available surface stations, data assimilation did appear to improve simulation of SWE in some regions. Not surprisingly, these were regions of ephemeral snow cover, where SCA observations are most indicative of SWE. The pull algorithm did produce an increase in precipitation in some regions. However, the CMAP-based forcing data used in the experiments had a low bias compared with other precipitation products, and the resulting total precipitation in the pull simulation was somewhere between CMAP and the others. Consequently, pull also tended to generate more runoff, latent heat flux, and soil moisture than open loop. Push, in contrast, had mixed hydrological effects, as the net impact of add/melt and remove/accumulate cycles somewhat offset each other when averaged across focus regions. These results could not be evaluated reliably against observations, because of the many possible sources of biases in modeled runoff and soil moisture, and hence the difficulty determining an improvement or degradation.
In high-latitude regions the hydrological benefits of data assimilation were limited by the fact that MODIS cannot observe surface conditions without daylight or during periods of extended cloud cover. This problem will be mitigated by the use of multisensor snow-cover products in future studies. Nonetheless, MODIS observations alone do provide significant information for many snow-covered regions of the world.
In summary, the shortcomings of existing snow-cover assimilation schemes motivated the development of the new pull algorithm. In particular, the algorithm overcomes a key limitation: in cases where the atmospheric forcing fields contradict the assimilated observation, updates to SWE can be quickly nullified by subsequent snowmelt or (less frequently) accumulation. The new scheme nudges the air temperature and precipitation into consistency with the near-future SCA observation, thus preserving the update. In our experimental simulation, this led to simulated SCA values that were higher than push and the open loop and that were more consistent with ground-based observations. The second shortcoming we addressed was the tendency of assimilation schemes to disrupt the local water balance. Pull inherently prevents nonphysical hydrological imbalance and limits the add/melt and remove/accumulate cycles that produce undesirable artifacts in the simulated snowmelt, soil moisture, and runoff fields. Nevertheless, there is still room for improvement in the application of data assimilation to snow hydrology on the global scale. More reliable satellite retrievals of SWE are needed to update snow fields in areas of 100% cover. In the interim, the application of multisensor SCA products will fill gaps in the MODIS product, allowing for smoother and more complete updates in high-latitude and cloud-covered regions.
This research was supported by NASA’s Energy and Water cycle Study (NEWS) program. The authors thank Dorothy Hall and George Riggs for helpful discussions on the MODIS SCA product, Sujay Kumar and Jim Geiger for assistance with LIS, and Hiroko Kato for advice on GLDAS simulations.
Corresponding author address: Benjamin F. Zaitchik, Dept. of Earth and Planetary Sciences, Johns Hopkins University, 3400 N. Charles St., Baltimore, MD 21218. Email: firstname.lastname@example.org