The benefit of using solar and longwave surface irradiance data from NASA’s Clouds and the Earth’s Radiant Energy System (CERES) synoptic (SYN) satellite product in simulations of snowmelt has been examined. The accuracy of the SYN downwelling solar and longwave irradiances was first assessed by comparison to measurements at NOAA’s Surface Radiation Network (SURFRAD) reference stations and to remote mountain observations. Typical shortwave (longwave) biases had magnitudes less than 30 (10) W m−2, with most standard deviations below 140 (30) W m−2. The performance of a range of snow models of varying complexity when using SYN irradiances as forcing data was then evaluated. Simulated snow water equivalent and runoff from cases using SYN data fell in the range of those from simulations forced with irradiances from well-maintained surface observation sites as well as empirical methods that have been shown to perform well in mountainous terrain. The SYN irradiances are therefore judged to be suitable for use in snowmelt modeling. It is also noted that the SYN upwelling shortwave irradiances, and hence albedos derived from them, are frequently not representative of individual monitoring stations because of the high spatial variability of snow cover and other surface properties in mountainous regions. In addition, adjusting the SYN downwelling longwave irradiances to reflect the exact elevation of the point of interest relative to the mean altitude of the satellite grid box is recommended.
Snowmelt supplies as much as 75% of surface water in basins of the western United States (Beniston 2006). However, as temperatures rise, more winter precipitation falls as rain rather than snow (Knowles et al. 2006), and snow-fed streamflow shifts earlier in the year (Stewart et al. 2005; Regonda et al. 2005). In the face of rapid warming, resource managers are now moving away from the use of historic records and toward hydroclimatic predictions for making decisions regarding reservoir management, city planning, and environmental restoration (e.g., Battin et al. 2007; Wiley and Palmer 2008). Basin-specific information regarding how local snow processes will respond to climate forcing is needed to inform these planning and management activities. In turn, accurate modeling of snow processes, particularly snow cover duration, is needed to treat boundary layer feedbacks in climate models.
Local snow water equivalent (SWE) depends on both accumulation and melt. Accurate simulation of accumulation requires precise precipitation estimates, including both amount and phase. Factors that affect the rate of snowmelt include incoming longwave (LW) and shortwave (SW or solar) irradiance, surface albedo, snow emissivity, snow surface temperature, sensible and latent heat fluxes, ground heat flux, and energy transferred to the snowpack from deposited snow or rain (Gray and Prowse 1993; Pomeroy et al. 2003). The net radiative flux generally constitutes about 80% of the energy balance (Male and Granger 1981; Marks and Dozier 1992; Cline 1997). Thus, the greatest potential sources of error in simulating snowmelt rates and timing are inaccurate solar and longwave radiative inputs.
In situ irradiance measurements, particularly of longwave irradiance, are rare in areas of seasonal mountain snowpack, and attenuation of solar (or shortwave) radiation by clouds is especially difficult to account for in snowmelt models (Male and Granger 1981). Most hydrologic models calculate the theoretical potential solar flux and estimate transmissivity from the local diurnal temperature range (e.g., Bristow and Campbell 1984; Hungerford et al. 1989; Thornton and Running 1999). Air temperature and relative humidity are often used to empirically estimate incoming longwave irradiance, with some methods assuming clear sky conditions and others using a measurement or approximation of shortwave irradiance to estimate a cloud cover fraction, which is used to adjust the longwave irradiance for clouds [see review of empirical longwave methods in Flerchinger et al. (2009)]. However, using these methods to compute net radiative inputs to a snow model can lead to overestimation errors of up to 50% in the snowmelt rates (Raleigh and Lundquist 2012; Raleigh 2013).
Further progress in snowmelt modeling thus requires the development of better sources of surface irradiance data. Satellite products have great potential to fill this gap because they provide measurement-based flux estimates over entire river basins, including locations where surface observation sites have not been installed. They also frequently incorporate longwave fluxes in addition to shortwave values. In this study, we examine the performance of satellite-based surface irradiance estimates from NASA’s Clouds and the Earth’s Radiant Energy System (CERES) synoptic cloud and radiation (SYN) product as forcing for snowmelt in both point simulations and basinwide hydrologic modeling.
We begin by describing the SYN dataset, which includes long- and shortwave up- and downwelling irradiances at the surface. We then present comparisons of SYN data with observations from reference and mountain measurement stations. We next apply irradiance time series from SYN and other sources as forcing data in simulations run with a range of models at several locations. Finally, we examine the outcomes of the model runs to assess whether SYN fluxes perform as well as fluxes from common empirical methods in snow modeling.
2. The CERES SYN data product
a. Product description
The CERES (Wielicki et al. 1996) synoptic cloud and radiation data products (Doelling et al. 2013; Rutan et al. 2015) are designed for use in studies of climate and the global or regional surface energy budget. They are available at three different temporal resolutions (monthly, diurnal averaged over a month, and 3 hourly) on a global 1° × 1° grid. In this study, we used solar and longwave downwelling surface fluxes from the Edition 3a CERES SYN1deg_3Hour dataset (hereafter referred to simply as SYN), which currently extends from March 2000 through July 2014.
SYN is based on measurements from instruments on board the NASA satellites Terra and Aqua. Because satellite instruments cannot directly measure irradiance at Earth’s surface from their position in the sky, the surface radiative fluxes in SYN are computed by applying a radiative transfer model to multiple inputs. For this reason, we do not refer to them as “observations” in this paper. Model inputs include MODIS cloud property data, surface snow and ice cover from NSIDC, and atmospheric profiles from the Goddard Earth Observing System model. The computed fluxes are then compared to top-of-atmosphere radiative fluxes from the CERES instruments and the input data tuned to yield matching results. Both the tuned and untuned fluxes are included in the SYN dataset. We used the untuned fluxes because tuning can introduce inconsistencies in the surface fluxes (D. Rutan 2012, personal communication).
Terra and Aqua are polar-orbiting satellites, which means that they each pass over a given location only twice per day, for a total of four daily observations. These overpasses are nominally at 1030 and 1330 local time, respectively, and the corresponding times at night. To create the 3-hourly product, measurements from Terra and Aqua are supplemented by radiances and cloud variables obtained at 3-hourly intervals from weather satellites around the world, including NOAA’s Geostationary Operational Environmental Satellites (GOES). Radiances from these satellites are cross-calibrated with coincident MODIS measurements so that the cloud properties retrieved from instruments on GOES and Terra and Aqua are comparable (Doelling et al. 2013).
Note that the 3-hourly values in SYN have been adjusted so that they correspond to averages over the 3-h time periods, such that averaging eight consecutive values will yield the correct total energy flux over a 24-h day. The times listed in the files indicate the beginning of the averaging periods. Although the underlying computations are performed for a quasi-equal area grid, the product is provided on a 1° × 1° grid. The increased spatial resolution at higher latitudes is provided by repeating the larger-scale irradiances in the corresponding smaller grid cells, not interpolation. SYN data can be downloaded in netCDF (from http://ceres.larc.nasa.gov/order_data.php) and in Hierarchical Data Format (HDF; from https://eosweb.larc.nasa.gov/project/ceres/ceres_table), although elevation information is only available in the HDF files at the time of this writing.
b. Corrections for elevation
Downwelling longwave irradiance is primarily determined by the vertical temperature and humidity profiles in the atmosphere above a station, which are a strong function of elevation (e.g., Ohmura 2001). Because the SYN grid cells cover an area of approximately 110 km × 110 km at midlatitudes, the average elevation over a grid box covering a mountainous area rarely matches the elevation of an individual point within this area. The SYN LW irradiances must thus be adjusted to make them representative of a particular station.
Accurate correction of LW irradiances for elevation requires knowledge of the specific meteorological conditions at the site of interest, which is rarely available except from atmospheric models. Thus, in this work we apply the constant irradiance gradient of −29 W m−2 km−1 estimated by Marty et al. (2002) using annual mean data collected over the period 1995–98 at 11 stations in the Swiss Alps. While this value is a climatological approximation for a specific geographic region, the gradient is sufficiently large that even a rough correction should improve the representativeness of our data. To give an idea of the magnitude of the corrections, we list the elevations of several well-studied mountain sites in the western United States, the mean elevation of the SYN grid boxes in which they fall, and the corresponding average longwave irradiance adjustments in Table 1. Snow study sites are often situated in complex terrain where the local elevation is greater than the gridcell average. For example, Reynolds Mountain East, Dana Meadows, and Lost Horse are all about 800 m higher than the average in the corresponding SYN grid box and are thus expected to receive 20–25 W m−2 less longwave irradiance than the grid box as a whole. In the Rocky Mountains of Colorado, differences greater than 1200 m occur at Niwot Ridge and the Senator Beck study plot, leading to irradiance offsets greater than 35 W m−2. On the other hand, the Central Sierra Snow Laboratory is 170 m below the gridbox mean, requiring a positive adjustment of about 5 W m−2.
Note that there is also a positive gradient of shortwave irradiance with elevation due to decreased scattering as the optical path through the atmosphere is reduced. Under clear skies, this gradient has been found to range between −7 and −26 W m−2 km−1 (Marty et al. 2002). However, in practice, it is modified significantly by local cloud, aerosol, and surface albedo conditions. When cloudy periods were included, Marty et al. (2002) found both positive and negative gradients as a function of elevation in their measurements, with a mean annual value of −13 W m−2 km−1. Because this value is location dependent and small relative to typical midday irradiance values, we have elected to neglect elevation correction for solar irradiance. In addition, prior work has shown that snow models are more sensitive to biases in LW than SW irradiance forcing data (Lapo et al. 2015a).
c. Data quality
Before applying the SYN fluxes as model inputs, it is useful to examine their accuracy relative to direct measurements. The CERES team has evaluated the SYN solar and longwave downwelling surface fluxes against measurements at 85 land and ocean sites around the world over the period from March 2000 to December 2007 (Rutan et al. 2015). The sites were selected based on quality and length of record and include locations in a variety of climates. The evaluation revealed a monthly mean bias of 2.2 W m−2 (1.1%) for the shortwave and −4.1 W m−2 (−1.2%) for the longwave irradiances. Standard deviations at the 3-hourly time scale were 55.0 W m−2 (28.4%) and 20.9 W m−2 (6.3%), respectively. Note that these comparisons included nighttime data, which reduces the average difference between satellite and surface measurements of solar irradiance.
For our application to snow simulation, it is particularly important to evaluate the satellite fluxes in mountainous areas and for snow-covered surfaces. In the study mentioned above, most sites were located in flat terrain and sites with high albedos were not broken out from the rest, so the specific conditions of interest were not addressed. It is well known that satellite data processing algorithms can incorrectly identify bright surfaces, such as deserts, ice, or snow, as clouds or vice versa (Rossow and Garder 1993; Sun and Cess 2005; Radkevich et al. 2013). The former error is more common (Q. Trepte 2015, personal communication) and can lead to a low (high) bias in downwelling shortwave (longwave) irradiance. In addition, low-resolution satellite products have difficulty representing highly variable terrain. Both of these factors suggest that SYN irradiances over snow-covered mountains may be of lower accuracy.
We performed additional comparisons at mountain stations and in the presence of snow to augment the CERES team’s evaluation of the SYN surface flux estimates. It must be noted that the quality of irradiance measurements is likely to be suboptimal when the sites are in remote locations and exposed to extreme environments. Inaccessibility makes it difficult to perform recommended maintenance and calibration procedures, such as the daily cleaning and inspection and annual calibration of instruments called for by WMO guidelines (WMO 2008). As a result, contaminants can build up on radiometer domes or sensitivity may decline over time (Brotzge and Duchon 2000). In addition, it is not uncommon for radiometers installed at mountain sites to be shaded by surrounding terrain or vegetation in the morning or evening, or for snow to accumulate on radiometer domes, particularly when lack of electrical power prevents the use of heating and ventilation systems (Betts and Ball 1997; Michel et al. 2008; Landry et al. 2014). Each of these problems causes negative biases in the recorded SW irradiance. (Snow and contaminants on the dome cause positive biases in the longwave.) Achieving horizontal alignment of instruments on towers is also difficult (Betts and Ball 1997). For these reasons, we were selective about the sites chosen for these comparisons and included some comparisons at climate reference stations. We also developed our own quality-control procedures for the mountain datasets used in this evaluation [see section 2c(2) for details]. Nevertheless, it is important to remember that differences between observations and satellite-retrieved irradiance values can originate with either of the two datasets [Zhang et al. (2010) discuss this issue further].
In addition to data quality issues at measurement sites, comparisons between satellite and surface irradiance data are confounded by the difference in spatial representativeness of the two data types. Surface observations record the irradiance received at an individual point, and although this is influenced by conditions over as much as several kilometers, it is clearly not the same as the mean irradiance over the ~110 km grid boxes of the SYN product. The standard approach to evaluating satellite data thus employs a temporal average of the surface observations, based on the assumption that movement of the clouds and other atmospheric constituents by the wind will allow the surface instrument to sample a representative portion of the satellite’s view. While performing such averaging has been empirically found to improve satellite–surface irradiance comparisons (Raphael and Hay 1984; Li et al. 2005; Su et al. 2008), no optimal relationship between spatial and temporal averages has been determined from theory. We skirt consideration of this issue by assuming that we are interested only in the practical performance of the satellite product. That is, would the available satellite data yield the same energy inputs as surface measurements when used as forcing data in a snow model? We thus compared the satellite-based irradiances to averages of the surface observations over the corresponding 3-hourly periods.
1) Evaluation with respect to high-quality reference sites
We began the evaluation using data collected at seven stations maintained by NOAA as climate reference sites. These stations make up the Surface Radiation Network (SURFRAD; Augustine et al. 2000, 2005) but are also part of the international Baseline Surface Radiation Network (BSRN; Ohmura et al. 1998). BSRN member stations are required to meet strict operating standards, including the use of high-quality instruments that are frequently calibrated and well maintained (McArthur 2005). At the SURFRAD sites, detailed in Table 1 and Fig. 1, irradiance is collected at 1-s intervals and provided as 3- or 1-min averages. Basic quality-control procedures are applied before the data are released. This dataset is available through NOAA’s Global Monitoring Division (http://www.esrl.noaa.gov/gmd/grad/surfrad/).
SURFRAD SW and LW irradiance data were averaged to 3-hourly periods matching the SYN intervals. Average irradiances below 10 W m−2 were excluded because this approaches the limit of radiometer instrument accuracy (Dutton et al. 2012). This also eliminated the comparison of shortwave values at night, when differences are by definition small. The remaining values were then compared to the SYN irradiances, with differences defined as SYN minus observation. Separate comparisons were performed for periods when snow was likely present at the measurement stations. These conditions were assumed to occur when the albedo determined from measured up- and downwelling shortwave irradiance was 0.5 or greater. Meeting this criterion does not imply that there was snow over the entire corresponding SYN grid cell, but it does indicate that there was some snow in the area. Metrics used in these evaluations were the mean and standard deviation of the differences as well as the mean absolute difference (MAD) and Pearson correlation coefficient r (Wilks 1995).
Results of these comparisons are shown in Table 2 and in Fig. S1 in the supplemental material. Here we see that the SYN shortwave irradiances are generally higher than the surface measurements, with biases of 3–23 W m−2 in magnitude. The standard deviations of the differences are much greater, at 75–125 W m−2. This reflects the fact that cloud conditions can vary greatly across an area the size of the SYN grid cells. Mean absolute differences are 50–85 W m−2, depending on the site. Correlation coefficients are all above 0.85, indicating that SYN captures most of the irradiance variability. We note that Desert Rock, Nevada, is the only site where SYN has a large negative bias. This is likely because of the bright surface or high amounts of reflective aerosols being interpreted as cloud at some times.
In Table 3 and in Fig. S1 in the supplemental material, the elevation-corrected SYN longwave data show all-around better agreement with the surface measurements than the shortwave irradiances do. The biases have magnitudes below 6 W m−2 and can be positive or negative. The standard deviations run from 17 to 26 W m−2, with MADs of only 13–20 W m−2, while correlations range between 0.90 and 0.94. In general, better agreement is expected in the LW because the controlling factors, such as temperature and water vapor profiles, are relatively uniform over large areas.
Under high-albedo (snow) conditions, SW biases all decrease or become more negative. This is consistent with snow being identified as clouds: both are highly reflective, but only clouds will reduce the SW flux received at the surface. The largest bias is −52 W m−2, while the others are all below 30 W m−2 in magnitude. The standard deviations and MADs are lower than those for the full dataset, but this is in part because SW irradiances are smaller during the winter—both statistics are larger when considered in relative terms. LW biases are generally larger or more positive when albedo is high, running from 1 to 19 W m−2. This may also be due to erroneous inclusion of clouds in the radiative calculations, since longwave emission from clouds is greater than that from clear air. Both the standard deviation and MAD are larger when the high-albedo times are broken out, with values of 27–31 and 21–24 W m−2, respectively. Note that the high-albedo samples are a tiny subset of the overall number of points because of limited snowfall at these locations and fewer daylight hours during the winter.
2) Evaluation against measurements from mountain sites
This portion of the evaluation utilized irradiance measurements from stations in mountainous areas of the western United States. Shortwave (longwave) irradiances were taken from 17 (10) different sites spread over the western United States and operated by a number of different institutions (see Table 1 and Fig. 1). As mentioned previously, quality-control procedures were developed for SW irradiances. These tested for nonphysical values, shading from surrounding terrain and vegetation, and the likelihood of snow buildup on the radiometer dome based on recent occurrence of snowfall and lack of sufficient wind or temperature to remove the snow (Lapo et al. 2015b). Samples failing these tests were not used in the evaluations. No special quality control was done for the LW irradiances because it is difficult to identify erroneous LW values unless they are extreme. However, LW data from periods when the SW radiometers were likely snow covered were eliminated from the LW comparisons. Irradiance data from the mountain stations were available as half-hourly or hourly means, which were averaged to produce 3-hourly means. Again, average values of 10 W m−2 or less were eliminated.
Results for the SW irradiance comparisons are given in Table 4 and in Fig. S2 in the supplemental material. While the statistics at many locations are comparable to those from the SURFRAD stations, they are much higher at a few sites. For example, the biases at four of the sites are 40 W m−2 or more. The range of standard deviations at the mountain stations is 93–157 W m−2, while the mean absolute differences run from 60 to 115 W m−2. These ranges overlap but also exceed those from the SURFRAD comparisons. The correlations are only slightly lower, at 0.84–0.94. The fact that the agreement between SYN and surface irradiance measurements is worse at the mountain sites could reflect lower measurement quality at these locations and/or the difficulties inherent to computing irradiance over a large area consisting of complex terrain. In addition, while removal of observations at times when radiometers are shaded can improve computed biases, it also tends to increase the average magnitude of random differences. This is because shading typically occurs in the morning or evening when irradiances, and hence possible differences, are small.
As for the SURFRAD comparisons, the agreement between observed and satellite-based LW irradiances (Table 5 and Fig. S2 in the supplemental material) is better than for the SW irradiances. However, it is still poorer at the mountain stations than at the SURFRAD sites. The biases at seven of the stations are in the single digits, like the SURFRAD data, but are −17, 18, and 31 W m−2 at the three other stations. The larger differences could be caused by poor representation of atmospheric conditions over variable terrain in the SYN LW computations or instrument error. The standard deviations (23–30 W m−2) and MADs (18–25 W m−2) fall within 5 W m−2 of the ranges found for SURFRAD except for standard deviations at SASP and SBSP, at 33 W m−2. Correlations are between 0.75 and 0.89.
3. Application to snowmelt models
Having investigated the basic accuracy of the SYN product, we assessed its performance in a range of snowmelt modeling scenarios. The approach taken here was to create multiple short- and longwave irradiance time series based on observations, empirical methods, and satellite remote sensing. These time series were then used as input to several snow models and the outcomes assessed using metrics such as snow disappearance date and melt rate. We first used two point models of differing complexity and sensitivity at test locations in a range of climates, in order to ensure that the behavior was not model or location specific. We then ran a distributed model for one river basin to examine how well the SYN irradiances represented a larger area.
Point simulations were performed using the Utah Energy Balance Snow Accumulation and Melt Model (UEB; Tarboton and Luce 1996) and the Snow Thermal Model (SNTHERM; Jordan 1991). Both of these models are physically based and have been used successfully in a range of hydrological studies (e.g., Cline 1997; Luce et al. 1998; Schulz and de Jong 2004; Watson et al. 2006). However, the representation of the snowpack and snow processes in these models is significantly different. SNTHERM is an energy balance model that represents snow and soil as a stack of multiple layers and uses an iterative process to compute the heat and energy content of each layer at each time step. SNTHERM explicitly models all of the physical processes that regulate snowpack evolution, including liquid water retention and densification. In contrast, UEB is a composite energy balance model that treats the snow and top layer of soil as a single bulk layer. It uses a linearized energy balance equation to compute energy exchange across the interface between this bulk layer and the atmosphere. While UEB includes a representation of liquid water retention, it assumes a constant snow density. For consistency in the representation of the shortwave energy balance, the U. S. Army Corps of Engineers (USACE; U. S. Army Corps of Engineers 1956) albedo scheme was integrated into both point models. Separate decay parameters were used for the accumulation and ablation seasons. For additional details about the physics of SNTHERM and UEB and how these affect model sensitivity to perturbations in irradiance, see Lapo et al. (2015a).
The Distributed Hydrology Soil Vegetation Model (DHSVM; Wigmosta et al. 1994; Cristea et al. 2014) was used in the basin-scale, distributed simulations. This model is of intermediate complexity between UEB and SNTHERM in that it represents the snowpack as two layers—a base snowpack that serves as a reservoir of mass and energy and can exchange heat and meltwater with the upper layer and ground, and an upper surface layer that interacts with the atmosphere through energy and mass transfer. In addition to computing the local snow energy balance, however, DHSVM simulates both overland and subsurface water flow. Preprocessed spatial distributions of vegetation, elevation, geology, soil depth and type, and terrain shading are required to represent the extended domain. In prior studies, DHSVM was found to represent snow accumulation and melt, streamflow, and interaction with vegetation accurately in complex terrain (e.g., Nijssen et al. 1997; Bewley et al. 2010).
b. Study sites
Test simulations were conducted at four sites in the western United States: Central Sierra Snow Laboratory (CSL; Osterhuber 1999; Taylor et al. 2001) on the crest of the Sierra Nevada, Lost Horse SNOTEL station (LOS; McMillan 1981) on the eastern slope of the Cascade Mountains, the Reynolds Mountain East sheltered site (RSH; Reba et al. 2011) in the Owyhee Mountains of Idaho, and the upper Tuolumne River basin in Yosemite National Park, also in the Sierra Nevada (Cristea et al. 2014). The locations of the CSL, LOS, and RSH sites are shown in Fig. 1, and a diagram of the upper Tuolumne River basin is given in Fig. 2. These sites were selected from a pool of well-studied locations in different climatological regions based on the availability over a common period of high-quality observations of the variables required for forcing and evaluation of the model runs. Air temperature, pressure, relative humidity, precipitation, wind speed and direction, downwelling solar irradiance, and SWE were observed at hourly intervals at CSL, LOS, and RSH, where point simulations were performed. Observations from the DAN, TUM, and TES sites, all within the Upper Tuolumne River basin, were pooled to obtain a full set of variables to force and evaluate DHSVM. (This was considered inappropriate for point modeling but reasonable when modeling an extended area.)
c. Irradiance forcing datasets
We selected several commonly used options for irradiance forcing for the comparisons, for a total of four SW and five LW irradiance series. The SW time series consisted of data from observations, the Bristow and Campbell (1984) empirical model (hereafter BC84) and Thornton and Running (1999) empirical model (hereafter TR99), and the CERES SYN product. BC84 is a widely used SW empirical model that relates incoming irradiance to the diurnal air temperature range, while TR99 is a version of BC84 that has been modified to be less dependent on local climate. As a result, TR99 was recently found to perform well across a wide range of climates (Bohn et al. 2013).
Longwave irradiance observations were not available from any of the selected stations, so three empirical methods were included. The datasets were based on Dilley and O’Brien (1998) with the Kimball et al. (1982) cloud correction (hereafter Dilley), Prata (1996) with the Crawford and Duchon (1999) cloud correction (hereafter Prata), Brutsaert (1975) with the Brutsaert (1982) cloud correction (hereafter Brutsaert), uncorrected CERES SYN, and CERES SYN corrected for local elevation. The three empirical methods were selected because they performed well in the comprehensive comparison of longwave irradiance estimation methods by Flerchinger et al. (2009). While Brutsaert was only the fourth best performing algorithm in that evaluation, the Brutsaert methods are by far the most cited longwave empirical methods. Each method was implemented as described in Flerchinger et al. (2009).
Point simulations of water-year (WY) 2009 (from 1 October 2008 to 30 September 2009) were performed for CSL, LOS, and RSH using both UEB and SNTHERM forced by all possible combinations of SW and LW irradiance data. The measured forcing data were ingested in their native 1-hourly format while the 3-hourly SYN and daily empirical irradiances were interpolated to this resolution. Limited tuning was performed for UEB on parameters related to accumulation, such as adjusting the undercatch factor and the rain versus snow parameters, so that the cumulative solid precipitation better matched observed SWE. All other parameters were unaltered from the default values. SNTHERM was used without alteration. Model performance was assessed using metrics derived from SWE observations at the measurement sites.
Distributed modeling runs using DHSVM were performed for the entire period of WY 2003–09 using 3-hourly time steps. Model setup and calibration followed the methods described in Cristea et al. (2014). Shortwave irradiance was distributed accounting for slopes and shadowing, while longwave irradiance was distributed uniformly across the domain. Both basinwide and point outputs of the model were evaluated, using SWE measured at the TUM and DAN study sites and streamflow at the Tuolumne River where it passes under the Highway 120 bridge as references.
Since the general outcomes with respect to the CERES data were consistent across all of the simulations, we concentrate on the distributed modeling here. In addition, we focus on the results from WY 2004 to simplify the presentation. However, averaged streamflow performance metrics for the entire WY 2003–09 period are also provided.
The mean daily values of the entire set of irradiance time series applied at Dana Meadows are shown for WY 2004 in Fig. 3. Here, the SYN shortwave and elevation-corrected longwave irradiances fall within the range of the observed and parameterized values. The relative consistency of this ordering from site to site is demonstrated by the annual mean irradiances given in Fig. 4. In nearly all cases, TR99 yields the highest shortwave irradiances followed by SYN, BC84, and the observations. Likewise, uncorrected SYN presents the highest longwave irradiances, followed by Prata, Dilley, corrected SYN, and Brutsaert. The exception is at CSL, the only location where the SYN LW irradiance required a positive adjustment for elevation. Because the Dilley LW irradiances consistently fall between the Prata and Brutsaert values, we reduced the computational effort by omitting Dilley from the distributed simulations.
Time series of SWE at the Dana and Tuolumne Meadows sites from the DHSVM runs are shown in Fig. 5, and metrics comparing measured and modeled SWE are given in Fig. 6. The SWE plots illustrate that the choice of LW irradiance data has a greater impact on the modeled SWE than the choice of SW data, since the lines are grouped together by LW choice (indicated by color) rather than SW choice (indicated by line style). In every case, the earliest melt [see snow disappearance date (SDD) in Fig. 6] occurs for the uncorrected SYN data, and the latest occurs when Brutsaert is used. The Prata and corrected SYN results fall in between, with Prata causing earlier melt than SYN. In terms of SW irradiances, the TR99 dataset always corresponds to the earliest melt and BC84 to the latest, with the SYN and observed SW irradiance cases in between.
Streamflow from the entire basin is shown in Fig. 7. Here, the plots are separated according to SW data source so that the lines are more easily distinguished. While the shapes of the curves are generally similar, large deviations are evident in the early part of the season (bracketing 1 April) and at the end-of-season tail. Consistent with the snow disappearance dates at DAN and TUM, the streamflow center timing (the date on which fractional cumulative discharge reaches 50%; Stewart et al. 2005), shown by the vertical lines in Fig. 7 and plotted in Fig. 8, is advanced for uncorrected SYN LW, delayed for Brutsaert, and close to the observation for both Prata and the elevation-corrected SYN LW data. Center timing is nearly the same for the observed and SYN SW irradiances, but advanced for TR99 and delayed for BC84.
Differences between measured and modeled melt rates from peak SWE to snow disappearance are illustrated in Fig. 6. In spite of the differences in timing evident in the SDD and streamflow center timing, the melt rates are quite consistent among the various simulations. In combination with the SWE plots, this suggests that most of the timing differences are due to offsets between the dates of peak SWE or in the peak amount of SWE rather than the rate of melt. For instance, differences in response to a brief early season melt event at TUM lead to variations in accumulated SWE that persist through the remainder of the snow season, ultimately affecting melt-out timing.
Wholesale shifts in the snow accumulation or melt onset mean that accuracy of the simulations as measured by the metrics of Figs. 6 and 8 is not simply a function of the ordering of the irradiance time series. For example, the SDDs are generally shifted late at Dana Meadows, with the result that the uncorrected SYN LW performs well, despite having the largest values. In fact, this dataset’s high bias offsets BC84’s negative bias, so that the SDD is closest to the measurement for this combination. The ordering of the results for the different irradiance datasets also varies from metric to metric. For example, for TR99 at TUM, the SYN and Prata SDDs are closest to the observation, but their melt rates compare the worst with the observation.
The SYN shortwave and corrected longwave time series performed as well in the point models as in the distributed simulations (Figs. S3–S5 in the supplemental material). However, the absolute accuracy of each simulation depended on both the location and model used, and, as for the distributed model results, varied from metric to metric. For example, UEB and SNTHERM perform very similarly at LOS but less so at CSL. At CSL, melt rate discrepancies are positive for UEB and negative for SNTHERM while SDD errors are of either sign for both models. Selection of “winners” and “losers” for both the models and the irradiance time series thus depends on which metrics are considered most important. Similarly, the exact snowmelt results over the full time period of the DHSVM runs (Fig. 8) depended on metric. However, the SYN SW and corrected LW metrics consistently fell within the range of the other datasets.
In a related study of the effect of radiative forcing errors on snow models, Lapo et al. (2015a) found that LW biases less than 20 W m−2 and SW biases less than 40 W m−2 led to root-mean-square differences (RMSDs) in SWE below 50 mm under most conditions investigated. Based on our analysis, both SYN SW and corrected LW irradiances meet these criteria at all of the SURFRAD sites (except SW in the winter at Bondville) and most of the mountain stations (13/17 in the SW and 9/10 in the LW). In the same study, random errors in irradiance with standard deviations up to 150 W m−2 were found to have little effect on snow modeling, and all of the standard deviations of the differences in this study were less than 150 W m−2 except for SW irradiance at USNR (157 W m−2). We therefore conclude that the SYN SW and elevation-corrected LW irradiances are of good quality. We note, however, that the above-mentioned criteria were based on a limited number of simulations at a few snow stations and may not translate to the same SWE RMSDs in other cases.
The discussion of modeling results above calls attention to the fact that it is difficult to assign true winners and losers among the irradiance datasets based on modeling results. This is similar to the problem of attempting to select a “best” overall snow model (Slater et al. 2001; Rutter et al. 2009; Essery et al. 2013). First, there are many different metrics of model performance, so that ratings depend on the user’s priorities. Second, the observations used as a reference for these metrics have associated error. Stream gauges require a rating curve to convert depth to discharge, and snow pillows can leak, for example. At the same time, there are uncertainties in the other measurements used to force the models, particularly precipitation, but also temperature and relative humidity. Poor modeling results may be caused by errors in these variables rather than the irradiance data. Alternatively, errors in different inputs or the model algorithms can compensate for each other to give the appearance of good performance. Although minimal model tuning was applied in this study, clearly tuning could affect the way a model performs using different inputs (Xia et al. 2005).
Regardless of the exact outcomes of the model runs, one certain result is that model performance using SYN SW and elevation-corrected SYN LW irradiance is not substantially different than that achieved when typical empirical irradiance formulae or surface observations are used. This can be traced directly to the fact that the SYN irradiances fall within the range of the observed or parameterized irradiances, as shown in Figs. 3 and 4.
Complications arise from the fact that SYN’s resolution is quite coarse by hydrological modeling standards. A SYN irradiance value may not be representative of a particular point falling within its grid cell. Large elevation differences can grossly affect downwelling LW fluxes. This problem can be compensated for by making an elevation-dependent flux adjustment. In this work, we have used a rough LW gradient value from a study in the Swiss Alps. Effecting a more accurate correction would require empirical data from the site of interest or knowledge of the local atmospheric temperature and water vapor profiles over time. Typical elevation corrections to the SYN LW irradiances here ranged from −40 to +6 W m−2 at the mountain stations, depending on their offset from the SYN gridcell mean elevation. The effect of elevation on downwelling SW irradiance is more moderate and was not corrected for here. Spatial variability of cloud cover can also cause differences between regional and local radiative fluxes.
While it may be tempting to use the CERES SYN down- and upwelling solar fluxes to determine albedo, which is also rarely, if ever, observed, we caution against this approach. Albedo is far more dependent on surface conditions than any other radiative energy budget term and is thus more susceptible to errors caused by rapid terrain and groundcover variations. Variability in snow cover, which is far brighter than any other surface type, is a particular problem. As an illustration, albedo computed at CSL during the winter of 2009 using the ratio of up- to downwelling solar irradiance from SYN is compared to albedo from the USACE parameterization in Fig. 9. Unlike the albedo from the USACE method, the SYN albedo never reaches 0.6 but fluctuates between 0.2 and 0.4 all winter. Clearly, values in this range are unrealistic for snow cover and must be the result of averaging over snow-covered and snow-free areas within the SYN grid cell.
Given that the SYN downwelling surface irradiances perform as well as irradiance values from empirical methods or surface observations when used as inputs to snowmelt models, they have several attractive features. First, they are available globally every 3 h from March 2000 to the present. The dataset is nearly spatially and temporally continuous; the few missing values tend to occur over high latitudes. This means that nearly 15 years of data are available at nearly any location of interest. The models used to compute the SYN fluxes are also consistent over time and space. Given radiometer maintenance issues at mountain stations, SYN fluxes may at times be more accurate than the available direct measurements—for example, when a radiometer is tilted or snow covered. At the same time, it is known that the uncertainty in empirically estimated irradiance grows as the distance between the modeled location and the meteorological measurement site increases (Mizukami et al. 2014). This means that the performance of irradiances from empirical models shown in this study, which were based on quality temperature and relative humidity measurements at the modeled location, is likely the best-case scenario. We note that SYN’s areal means may actually be more appropriate for large-scale distributed modeling than forcings based on measurements at a single observation station, given that this point may not be representative of the entire region. For example, it may be shaded by local surface features or vegetation. An additional problem with empirical models is that they are fit to historic conditions in a subset of possible climates and require retuning to obtain good fidelity in other locations. CERES SYN has the advantage of being more physically based and tested under a wide range of surface and atmospheric conditions. Thus, it makes sense to prefer SYN over empirical models for use in “new” or poorly studied regions. In addition, because it is not tied to historical measurements, SYN will include any long-term variations in irradiance that occur because of changes in aerosols, clouds, or other climatological factors. This is important as hydrologic conditions shift away from historical patterns.
The suitability of short- and longwave surface irradiance values from the CERES 3-hourly 1°SYN data product for use as forcing data in snowmelt models has been investigated. Comparisons show that the mean differences between observations and SYN SW irradiance are generally less than 30 W m−2, with the exceptions occurring at mountain sites where the quality of observations is uncertain. LW differences are smaller, typically below 10 W m−2, again with exceptions at certain mountain stations. SYN SW values tend to be smaller under conditions of high surface albedo (snow), likely because of misidentification of the bright surface as clouds. LW values tend to be higher for the same reason. Because of SYN’s low resolution, the LW irradiances should be corrected for differences in elevation between the location of interest and the mean elevation of the SYN grid box. Spatial variability of surface albedo within grid boxes also means that SYN upwelling SW irradiances may not accurately reflect local conditions, especially during the winter.
Performance of the SYN data in snowmelt models was thoroughly assessed using two point models and one distributed model and four different test locations. Runs using irradiance observations and irradiances computed using typical empirical methods were used as references for the SYN simulations. The simulations were run for all combinations of the selected LW and SW datasets and the results evaluated using a range of metrics such as snow disappearance date, seasonal melt rate, and RMSD between measured and modeled SWE. While specific results varied from model to model and run to run, both the SYN SW and LW data products were found to perform within the range of typical empirical methods and SW irradiance observations. We thus recommend their use in hydrological modeling. Given the advantages of complete geographic availability and a physical basis, the satellite products are especially useful for locations where no other irradiance data are available.
We thank David Rutan and David Doelling of the CERES production team at the NASA Langley Research Center for patiently answering numerous questions about the SYN product. Thanks are also due to the Lundquist Mountain Hydrology research group for their suggestions and helpful review of this work as it progressed. In addition, we would like to recognize the efforts of the dedicated staff who maintain the observation stations that provided irradiance and other data for this study. Funding for this project was provided by NASA Grant NNX11AF54G as part of the Science of Terra and Aqua program. K.E.L. also benefitted from support from the University of Washington Program on Climate Change Graduate Fellowship and a NASA Earth and Space Science Fellowship, NASA Grant NNX13AN78H. CERES data were obtained from the NASA Langley Research Center Atmospheric Science Data Center.