[1] Sixteen years of high‐quality surface radiation budget (SRB) measurements over seven U.S. stations are summarized. The network average total surface net radiation increases by +8.2 Wm −2 per decade from 1996 to 2011. A significant upward trend in downwelling shortwave (SW‐down) of +6.6 Wm −2 per decade dominates the total surface net radiation signal. This SW brightening is attributed to a decrease in cloud coverage, and aerosols have only a minor effect. Increasing downwelling longwave (LW‐down) of +1.5 Wm −2 per decade and decreasing upwelling LW (LW‐up) of −0.9 Wm −2 per decade produce a +2.3 Wm −2 per decade increase in surface net‐LW, which dwarfs the expected contribution to LW‐down from the 30 ppm increase of CO 2 during the analysis period. The dramatic surface net radiation excess should have stimulated surface energy fluxes, but, oddly, the temperature trend is flat, and specific humidity decreases. The enigmatic nature of LW‐down, temperature, and moisture may be a chaotic result of their large interannual variations. Interannual variation of the El Niño/Southern Oscillation (ENSO) ONI index is shown to be moderately correlated with temperature, moisture, and LW‐down. Thus, circulations associated with ENSO events may be responsible for manipulating (e.g., by advection or convection) the excess surface energy available from the SRB increase. It is clear that continued monitoring is necessary to separate the SRB's response to long‐term climate processes from natural variability and that collocated surface energy flux measurements at the SRB stations would be beneficial.

1 Introduction [2] The difference between the downward and upward components of shortwave (SW) and thermal infrared longwave (LW) irradiance at the surface defines the surface radiation budget (SRB). The SRB represents absorbed energy at the surface that is partitioned among latent, sensible, and ground heat fluxes and is the bulk of the energy available for atmospheric dynamic processes. Resultant evaporation and atmospheric heating drives circulations at all scales from the small convective eddies of the surface layer to the large gyres of the earth's general circulation and ocean currents that are important to short‐term climate variability and climate change. The direct radiative response to increasing greenhouse gases should be revealed by systematic changes in the SRB, especially in LW‐down. It follows that the success of modern climate models and satellite algorithms that attempt to estimate the SRB, or assess the state of the Earth's climate, hinges on correctly modeling the SRB and surface energy fluxes. Hence, highly accurate and continuous SRB monitoring for satellite and model validation, as well as for a host of other applications, is essential. [3] To fill that niche, accurate and precise SRB networks such as SURFRAD [Augustine et al., 2000; 2005], the Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) Program Southern Great Plains site [Stokes and Schwartz, 1994], and other Baseline Surface Radiation Network (BSRN) stations worldwide [Ohmura et al., 1998] were established in the mid‐1990s. These networks have since served as validation for modern programs such as NASA's Earth Observing System (EOS) [Wang et al., 2009; Liu et al., 2009], NOAA's current and future satellite systems (e.g., GOES‐R, JPSS) [Yu et al., 2009], and modern climate and weather forecast models [e.g., Markovic et al., 2008; Yang et al., 2008; Morcrette, 2002]. They also have been used extensively in hydrology, biological, and renewable energy applications [e.g., Maurer et al., 2002; Hall et al., 2009; Renne et al., 2008]. [4] The seven U.S. SURFRAD stations (Figure 1) are distributed among diverse climates that include the hot, dry desert of the southwestern United States, the moist environment of the eastern United States, and the semiarid northern plains. Sites were chosen to be regionally representative. In this study, long‐term averaging acts to increase the area of influence of the stations. For example, Augustine et al. [2006] used the method of Dutton et al. [2006] to cross‐correlate annual averages of downwelling shortwave (SW‐down) with satellite‐based estimates of surface radiation from the International Satellite Cloud Climatology Project (ISCCP). They found that, collectively, annual averages of SURFRAD station data represent a large part of the southwestern United States, the region from the Front Range of the Rocky Mountains to the Appalachians, and a large part of New England. Long‐term averages also reduce the data volume and act to remove the random error of the measurements. The operational nature of the network as opposed to a campaign approach is beneficial because the surface radiation environments of all types of atmospheric and surface conditions, both expected and unexpected, are continuously sampled and archived, building a high‐quality data set from which climate variability, climate response, and random atmospheric phenomena [e.g., Stone et al., 2011], can be studied. A valuable consequence of the establishment of SURFRAD is that an unprecedented and exceptional multiyear data set of the highest quality has evolved. Never before has the SRB in the United States been monitored so accurately with such spatial and temporal coverage. In this article, we document the variability of the SRB and its components at the U.S. SURFRAD stations from 1996 through 2011. Figure 1 Open in figure viewer PowerPoint The SURFRAD surface radiation budget network. Stations at Bondville, Fort Peck, Goodwin Creek, and Table Mountain began operation in 1995; Penn State and Desert Rock began in 1998; and Sioux Falls was added in 2003.

2 Measurements 2.1 Instruments [5] Substantial effort is directed toward achieving the most accurate and precise measurements possible at SURFRAD stations. The four basic components of the SRB, SW‐down, SW‐up, LW‐down, and LW‐up, are measured independently by using radiometers with thermopile detectors. The advantage of a thermopile detector for SW measurements is that it senses the full solar spectrum from the ultraviolet through the near infrared (280 − 4000 nm), as opposed to more widely used silicon cell detectors that cut off sensitivity at 1100 nm and miss much of the variability caused by water vapor absorption in the near infrared, leading to greater uncertainty. In SURFRAD, SW‐down is derived primarily from independent direct and diffuse SW measurements made by an Eppley Normal Incidence Pyrheliometer model NIP and an Eppley model 8‐48 pyranometer, respectively. Given that the solar tracker on which these instruments are mounted and the instruments themselves are well maintained, the “component sum” of direct + diffuse is preferable to a single pyranometer measurement [Ohmura et al., 1998; Gueymard and Myers, 2008; Bush et al., 2000]. However, a single pyranometer (Spectrosun model SR‐75) measurement is also made as a backup for times when the solar tracker fails. [6] An Eppley model PIR pyrgeometer is used to measure LW‐down and is shaded according to the specifications of the BSRN [McArthur, 2004; Ohmura et al., 1998]. Shading the instrument dome cuts out the small amount of thermal infrared contained in the solar beam and inhibits uneven solar heating of its mirrored surface, thus preventing associated errors [Philipona et al., 1995]. SW‐up and LW‐up are measured at the 8‐m level of a 10‐m tower by an inverted Spectrosun model SR‐75 pyranometer and an inverted Eppley PIR, respectively. To aid research, aerosol optical depth (AOD), cloud coverage, and atmospheric state measurements are also made. AOD is derived from spectral SW measurements made by a MultiFilter Rotating Shadowband Radiometer [Harrison et al., 1994], and cloud coverage is measured with a Total Sky Imager (TSI), both of which are manufactured by Yankee Environmental Systems, Inc. From 1995 through 2008, SURFRAD data were reported as 3‐minute averages of 1‐second samples; 1‐minute averages of 1‐second samples have been compiled since 1 January 2009. 2.2 Calibration, Quality Assurance, and Data Quality Control [7] For long‐term analysis it is important that the measurements do not introduce erroneous results owing to poor maintenance, systematic errors, or calibration drift. To avoid such problems, quality assurance is stressed. Radiometers are cleaned regularly, exchanged annually, and recalibrated before deployment. Shortwave instruments have been calibrated by the National Renewable Energy Laboratory in Golden, Colorado, and the World Meteorological Organization's (WMO) Region 4 Regional Solar Calibration Center/NOAA in Boulder, Colorado. Calibrations from both of these facilities are traceable to the World Radiometric Reference (WRR) [Fröhlich, 1977]. Calibration uncertainty for the minute‐scale measurements varies by instrument and ranges from ±2% to 5% for pyranometers and from ±1% to 2% for pyrheliometers [Michalsky et al., 2011]. Pyrgeometer calibrations are traceable to the World Infrared Standard Group (WISG) at the World Radiation Centre (WRC) in Davos, Switzerland (http://www.pmodwrc.ch/pmod.php?topic=irc). The WISG is made up of two Eppley PIR pyrgeometers and two Kipp & Zonen CG4 pyrgeometers. Absolute calibrations of those instruments are traceable to the IPASRC‐I experiment [Philipona et al., 2001], in which the WISG standards were calibrated using an Absolute Sky‐Scanning Radiometer. For SURFRAD, a standard group of three pyrgeometers is sent to the WRC for calibration on a quasi‐annual basis. Uncertainty of the WRC calibration is ±4.2 Wm−2. Operationally, SURFRAD field pyrgeometers are run alongside the SURFRAD standard group for calibration before deployment. That transfer increases the uncertainty of the field pyrgeometer measurements to about ±5 Wm−2, which equates to ±2% for a 300 Wm−2 signal. Uncertainties of other measurements made at SURFRAD stations have been discussed by Augustine et al. [2000; 2008]. [8] SURFRAD data are quality controlled on a daily basis by automated methods and by a human analyst before being released. The result of careful calibration, annual instrument exchange and maintenance visits, other quality assurance procedures, and data quality control is that SURFRAD has compiled one of the most accurate, precise, and continuous surface radiation data sets in the world.

3 Data Analysis 3.1 Data Reduction [9] Minute‐scale data are superfluous for examining interannual variability. To reduce the data volume, monthly averages were generated from the elemental minute‐scale measurements, and annual means were computed from the monthly data. Because of the diversity of climate across the network and to isolate potential continent‐scale phenomena, annual averages were normalized to zero by computing annual anomalies. Only station years with a full 12 months of data were used. For example, because the Desert Rock, Nevada, station was installed in March 1998, only data from 1999 onward were used. Although the SURFRAD data archive commences in 1995, this analysis begins in 1996 because of multiple startup issues during the first year. The 16‐year data set is fairly complete, and only one station‐year has been excluded. In 2004 the Goodwin Creek station was struck by lightning and remained nonfunctional for 6 weeks. No measurements were available during that period, which precluded any gap‐filling attempts. With 2 months of all measured quantities missing, meaningful annual averages for Goodwin Creek for 2004 could not be computed. [10] To ensure high‐quality annual averages, erroneous measurement biases have to be avoided or corrected. For example, SW‐down measurements from single pyranometers are adversely affected by cosine errors for times when the solar beam illuminates the sensor [Michalsky et al., 1995; Vignloa et al., 2012]. Cosine errors vary as a function of solar zenith angle (SZA) but are ignored in practice because only the responsivity corresponding to 45° SZA is applied. That practice introduces small “cosine” errors in single‐pyranometer SW data at all other SZAs. Also, thermal offset errors caused by infrared cooling of the glass dome and sensor surface from a cold sky are common for pyranometers that have solid black thermopile detectors [Gulbrandsen, 1978; Dutton et al., 2001]. Those errors are nominally on the order of 0 to −10 Wm−2 and manifest themselves as negative values at night but also erroneously reduce the thermopile voltage output slightly during the daytime, especially under clear skies. For this analysis, cosine and thermal offset errors are avoided because the primary source of SW‐down is independent measurements of the direct and diffuse solar components. Pyrheliometers that measure direct solar have a small field of view and are always pointed directly at the sun and, therefore, have neither cosine nor thermal offset problems. By definition, diffuse measurements do not include the solar beam and thus have little or no cosine error, given that they are calibrated at 45º [Vignloa et al., 2012]. After 1999, the Eppley model 8‐48 (also known as “black and white”) pyranometer, which does not have a significant thermal offset error (<1 Wm−2), was used to measure diffuse SW. Thermal offsets in diffuse SW measurements made before 1999 with solid‐black‐detector pyranometers were corrected using the method described by Augustine et al. [2005], which is based on work of Dutton et al. [2001]. SRB calculations for this analysis include only the backup global solar measurement made with a single black‐detector pyranometer for times when the solar trackers were not functional. The upwelling SW measurement does not have a thermal offset because the inverted pyranometer views the ground surface, which is nearly at the same temperature or warmer than the thermopile sensor. After minimizing the thermal offset errors in practice, any residual offsets at night were preemptively set to zero prior to long‐term averaging. [11] Longwave data have no offset issues because the temperature of the pyrgeometer dome is measured and accounted for in the signal processing. SURFRAD network pyrgeometers used for the LW‐down measurements prior to the instrument exchanges in 2000 were not shaded, as required by the BSRN, so those data include the small amount of thermal infrared contained in the solar beam. For an annual mean, that contribution to LW‐down is approximately 1 Wm−2, as projected onto a horizontal surface and accounting for mean cloud effects. For this analysis, that small contribution was removed from the LW‐down annual averages prior to 2001. 3.2 Gap Filling [12] Data gaps were not a large problem because SURFRAD stations have had little down time. Table 1 summarizes data coverage for all measurements relating to the SRB. Well over half of the station‐months have 100% data coverage (column 3, Table 1), and most are represented by more than 95% of the possible data (column 4, Table 1). Station months with 70 − 94% coverage represent only 3 − 5% of the data set. Computed averages for those months were used as is, but typically they had data coverage in the upper 80% to low 90% range. Temporal consistency of the SRB components was checked by comparing monthly averages of total surface net radiation computed by two different approaches. First, all four SRB components, SW‐down, SW‐up, LW‐down, and LW‐up, were summed at the minute‐scale and averaged over monthly periods. Next, monthly means of surface net‐SW and net‐LW were computed separately from the minute‐scale data and added to get the second representation of the total surface net radiation monthly mean. Monthly averages from those two approaches are plotted against each other in Figure 2 for each station. When they differ, the net‐SW and net‐LW components represent different periods, and, when they agree, the SW and LW measurements that contribute to the monthly total surface net radiation represent similar periods. In Figure 2 only a few points deviate from the 1:1 line, indicating that most of the independently measured upwelling and downwelling SW and LW signals are consistent in time and thus sampled the same physical processes on the monthly time scale. Table 1. Percentage of Months With Designated Data Coverage for Specified Surface Net Radiation Quantities and Stations Station Quantity 100% Data Coverage ≥95% Data Coverage 70 − 94% Data Coverage <70% Data Coverage Bondville Net‐SW 83.8 96.4 2.1 1.6 1996 − 2011 Net‐LW 76.6 93.8 5.2 1.0 192 Months Total net 75.0 93.8 4.2 2.1 Fort Peck Net‐SW 77.6 94.3 1.6 4.2 1996 − 2011 Net‐LW 72.9 90.1 6.3 3.6 192 Months Total net 67.2 88.5 6.8 4.7 Goodwin Creek Net‐SW 76.1 92.2 5.6 2.2 1996 − 2011 Net‐LW 73.9 91.7 5.0 3.3 180 Months Total net 70.6 88.9 7.8 3.3 Table Mt. Net‐SW 71.9 97.9 2.1 0.0 1996 − 2011 Net‐LW 64.6 93.2 4.7 2.1 192 Months Total net 56.3 92.7 5.2 2.1 Desert Rock Net‐SW 75.6 93.0 6.4 0.6 1999 − 2011 Net‐LW 67.3 89.2 5.1 5.8 156 Months Total Net 65.4 87.2 6.4 6.4 Penn State Net‐SW 78.9 97.4 1.3 1.3 1999 − 2011 Net‐LW 82.1 96.2 1.3 2.6 156 Months Total net 76.3 94.9 2.6 2.6 Sioux Falls Net‐SW 83.3 99.0 1.0 0.0 2004 − 2011 Net‐LW 56.3 97.9 2.1 0.0 96 Months Total net 55.2 96.9 3.1 0.0 Figure 2 Open in figure viewer PowerPoint Data temporal consistency test for each station that compares monthly averages of the surface radiation budget computed in two ways: (1) the sum of independently computed monthly averages of net shortwave and net longwave (y‐axis), and (2) monthly averages of the minute‐scale four‐component sum (x‐axis). One‐to‐one correspondence indicates temporal consistency among the SW and LW measurements within a month. [13] Station monthly averages represented by less than 70% data coverage, which account for only 1.5 − 3.2% of the data set, were not used. Rather than using simple linear interpolation to fill those gaps, those averages were reconstructed from monthly means of other variables using physically based associations, empirical relationships, or alternative data sets. Making use of physical linkages among measured quantities is justified because monthly means are generally well behaved in a physical sense and are relatively free of random error. For example, missing SW‐up for July and August of 2009 at Fort Peck was filled by using the physical relationship that links SW‐down, SW‐up, and albedo. The monthly average SW‐up for each of those months was estimated by multiplying the monthly mean of SW‐down times the mean surface albedo for June of that year. This procedure was justified because June's average albedo (0.171) and that for September (0.175) were similar, indicating that the missing albedo for July and August should be well approximated by the June value. Missing LW‐down for each of those two months at Fort Peck was estimated empirically from a linear least squares fit between monthly LW‐up and LW‐down from other months of that year (R2 = 0.96). Similarly, missing LW‐down for June and for July 2003 at Goodwin Creek were filled by using a second‐order polynomial fit between the monthly means of air temperature and LW‐down from the other months of 2003 (R2 = 0.98). When possible, missing monthly averages of temperature, pressure, and relative humidity were recovered from collocated data sources (e.g., the USGS SCAN Network data at Goodwin Creek). When such resources were not available, NCEP Reanalysis data were used. As shown in the last column of Table 1, only a small percentage of station month averages were artificially filled by these methods.

4 Results [14] Results are displayed as time series of annual anomalies from the long‐term annual mean. Most time series show the network mean annual anomaly as a thick black curve (except in Fig. 4). Instead of error bars, time series of all contributing stations are shown. The scatter of the individual stations visually portrays the uncertainty of the network mean, as well as the coherence of detected signals among the various stations. Linear fits are applied to the annual anomaly time series only to quantify per‐decade changes and to normalize the results for comparison. Those linear tendencies are not meant to be extrapolated, because any trends, even those that are statistically significant, are likely to be part of longer‐term cycles that cannot be resolved in a 16‐year period. 4.1 Total Net Surface Radiation [15] Annual anomalies of total surface net radiation from 1996 through 2011 are shown in Figure 3. The network mean as well as the means of the individual stations generally increase over the entire period. A linear least squares fit to the network mean reveals a +8.2 Wm−2 per decade increase, which amounts to about +13 Wm−2 over the entire period. That fit has a correlation coefficient of 0.89, accounting for 79% of the variance, and its trend passes the nonparametric Mann‐Kendall (MK) test [Mann, 1945; Kendall, 1975] for significance at the 95% level (Table 2). Figure 3 Open in figure viewer PowerPoint Time series of total surface net radiation annual anomalies. Color curves represent the individual stations, and the heavy black curve is the network mean. The dashed line is a least squares linear fit to the network mean for which the best‐fit equation and coefficient of determination (R2) are given at upper right. Table 2. Mann‐Kendall Trend Significance Results Quantity MK Median Trend 95% Lower OLimit Trend 95% Upper Limit Trend MK Significance Total net radiation +0.86 +0.57 +1.08 Significant Net‐SW +0.53 +0.27 +1.08 Significant Net‐LW +0.20 −0.05 +0.49 Not significant SW‐down +0.72 +0.26 +1.11 Significant SW‐up +0.07 −0.09 +0.29 Not significant LW‐down +0.26 −0.13 +0.52 Not significant LW‐up −0.08 −0.43 +0.30 Not significant [16] To determine the source of the observed increase, components of the total surface net radiation are examined separately. Total surface net radiation is made up of four quantities, SW‐down, SW‐up, LW‐down, and LW‐up, each of which responds to different stimuli. By convention, downward radiation is a source and upward radiation is a loss. Therefore, the total surface net radiation is defined as: (1) To determine the source of the observed increase, components of the total surface net radiation are examined separately. Total surface net radiation is made up of four quantities, SW‐down, SW‐up, LW‐down, and LW‐up, each of which responds to different stimuli. By convention, downward radiation is a source and upward radiation is a loss. Therefore, the total surface net radiation is defined as: [17] Time series of network mean annual anomalies of surface net‐SW (SW‐down – SW‐up) and surface net‐LW (LW‐down – LW‐up) are plotted separately in Figure 4. For clarity, individual station time series, as shown in other figures, are not included. These basic components of the SRB are much more variable than their sum in Figure 3, but that is expected because surface net‐SW and net‐LW are inherently anticorrelated as a result of cloud effects. That behavior is obvious in Figure 4. Greater surface net‐SW generally means more SW input from reduced cloud cover and/or aerosols. Alternately, fewer or optically thinner, clouds lead to less surface net‐LW, i.e., more LW‐up because of a warmer surface from greater SW heating and less LW‐down because of a colder (less cloudy) sky. The opposing nature of surface net‐SW and net‐LW combined with their interannual variability increases scatter, causing linear least squares fits to their annual anomalies to have lower correlation coefficients than their sum (Figure 3). According to the MK results (Table 2), the trend of +5.6 Wm−2 per decade for surface net‐SW is significant at the 95% level, but that for surface net‐LW (+2.3 Wm−2 per decade) is not. The difference in magnitude and statistics of these two slopes indicates that the increase in surface net‐SW dominates the change in the SRB over the United States from 1996 through 2011. Also, that the sum of those two individual slopes (+0.80 Wm−2 per decade) is close to the independently computed slope of the total surface net radiation in Figure 3 (+0.82 Wm−2 per decade) is a testament to the quality and consistency of the data set. Figure 4 Open in figure viewer PowerPoint The blue curve represents annual anomalies of the network mean of surface net‐SW (SW‐down – SW‐up). The red curve represents annual anomalies of the network mean of surface net‐LW (LW‐down – LW‐up). Dashed lines are least squares linear fits to each time series for which the best‐fit equations and coefficients of determination (R2) are given at upper right. 4.2 Surface SW Net Radiation [18] The increase in surface net‐SW shown in Figure 4 is dominated by SW‐down (Figure 5a), which rises relatively steadily at a mean rate of +6.6 Wm−2 per decade (+11 Wm−2 over the 16‐year period) and is statistically significant according to the MK test (Table 2). Figure 5b shows that a slight increase in the network mean of SW‐up of +1.1 Wm−2 per decade, which represents a loss at the surface, slightly tempers the surface net‐SW. At least two sources contributed to the increase of the SW‐up. First, SW‐up inherently escalates with increasing SW‐down, and, second, there was greater mean reflectance in the latter part of the 16‐year period from anomalous snow cover in 2007 and 2010. The largest peak in SW‐up occurs in 2010 when the eastern two‐thirds of the United States experienced unusually snowy weather and snow cover extent [Blunden et al., 2011]. Note in Figure 5b that the eastern U.S. stations contribute most to the 2010 peak. Removing 2007 and 2010 from the SW‐up time series diminishes its slope to +0.3 Wm−2 per decade. Another anomalous feature in Figure 5b is a prolonged minimum of SW‐up at the Fort Peck station from 2005 through 2008, which was caused by four consecutive dry winters there. However, the other stations’ normal or higher than normal reflectance during the later part of the period overwhelm the influence of the extended minimum at Fort Peck. Figure 5 Open in figure viewer PowerPoint Annual anomalies of the shortwave and longwave components of the total surface net radiation for individual stations (colored curves) and the network mean (black solid curves) in time series. The four panels represent SW‐down (a), SW‐up (b), LW‐down (c), and LW‐up (d). Dashed lines are least squares linear fits to the network annual means for which the best‐fit equation and coefficient of determination (R2) are given at upper right in each panel. [19] Decreased aerosols and/or cloud coverage could cause the observed solar brightening. Although Wild [ 2012 Long et al. [ 2009 Long et al. [ 2009 Augustine et al., 2008 Long et al., 2006 aerosols ) using (2) S o is the extraterrestrial downward SW flux at the top of the atmosphere at 40°N (343 Wm−2) as described by Dutton and Cox [ 1995 o is a representative single scattering albedo of continental aerosols (0.97), and g is the typical asymmetry parameter for continental aerosols. The value of g (0.87) was obtained from Paltridge and Platt [ 1976 −2 over the 16‐year period. This value was slightly mitigated by a small increase in stratospheric aerosols over the last decade. From its start in 1995, the SURFRAD network sampled the surface radiation environment without any major volcanic perturbations other than a minor residual from the explosive Mt. Pinatubo eruption in 1991. Several consecutive minor eruptions in the last decade [Vernier et al., 2011 Solomon et al., 2011 −2 per decade decrease in SW‐down, or −0.16 Wm−2 over the entire period. It is clear that the SW radiative effect at the surface attributed to the steady decrease of AOD (+0.82 Wm−2) is negligible compared with the +11 Wm−2 observed increase of SW‐down over the 16‐year period, by deduction supporting the dominant role of clouds in the observed brightening. Decreased aerosols and/or cloud coverage could cause the observed solar brightening. Although] attributed solar brightening over much of the globe after 1990 to decreasing aerosols,. [] convincingly showed that the increase in SW‐down over the United States from 1995 to 2007 was due to a general decrease in cloud coverage. Furthermore,. [] show that the documented decrease in AOD over the United States during that period [.,] had little to do with the observed brightening. Using the algorithm of.,, we extended the same analysis to 2011 (not shown) and demonstrated that cloud coverage continued to decrease after 2007. Furthermore, from surface measurements, we documented a network‐wide systematic decrease in surface specific humidity of −0.3 g/kg over the 16‐year period, which is consistent with a decrease in clouds or cloud optical depth. With regard to the effect of aerosols, the time series of the network annual average 500‐nm AOD from 1997 through 2009 is shown in Figure 6 . AOD for 1995, 2010, and 2011 are excluded because of incomplete data, but, over the 13 years analyzed, the network mean 500‐nm AOD decreased by −0.016 per decade, which extrapolates to −0.025 over the 16‐year period. The mean impact of that change in AOD on SW‐down at the surface was estimated by computing the mean aerosol SW radiative forcing (SW RF) usingwhereis the extraterrestrial downward SW flux at the top of the atmosphere at 40°N (343 Wm) as described by], α is our computed network mean surface albedo (0.2522), Δτ is the measured change in 500‐nm AOD over 16 years, ωis a representative single scattering albedo of continental aerosols (0.97), andis the typical asymmetry parameter for continental aerosols. The value of(0.87) was obtained from]. From equation 2 , the effect on SW‐down of a 0.025 decrease in 500‐nm AOD is +0.82 Wmover the 16‐year period. This value was slightly mitigated by a small increase in stratospheric aerosols over the last decade. From its start in 1995, the SURFRAD network sampled the surface radiation environment without any major volcanic perturbations other than a minor residual from the explosive Mt. Pinatubo eruption in 1991. Several consecutive minor eruptions in the last decade [.,.,] produced a slight increase in stratospheric AOD of about +0.004 per decade, which, according to equation 2 accounts for a −0.1 Wmper decade decrease in SW‐down, or −0.16 Wmover the entire period. It is clear that the SW radiative effect at the surface attributed to the steady decrease of AOD (+0.82 Wm) is negligible compared with the +11 Wmobserved increase of SW‐down over the 16‐year period, by deduction supporting the dominant role of clouds in the observed brightening. Figure 6 Open in figure viewer PowerPoint Time series of the network annual average 500‐nm aerosol optical depth with linear least squares fit for which the best‐fit equation and coefficient of determination (R2) are given at upper right. Error bars indicate standard error. 4.3 Surface LW Net Radiation [20] Slopes of linear fits to the time series of the net‐LW network mean annual anomaly (Figure 4) and its upward and downward components (Figures 5c and 5d) are smaller than those of SW‐down and total surface net radiation. None of the LW time series passes the MK test for significance at the 95% confidence level (see Table 2). Therefore, we discuss the mean tendencies of the LW annual anomaly time series only in terms of their variability and how they contribute mathematically to the total surface net radiation. [21] Network‐wide surface net LW gradually increases by +2.3 Wm−2 per decade (Figure 4) but has a high degree of variability (σ = 2.14 Wm−2). Its downward component (LW‐down; Figure 5c) also demonstrates large interannual variability (σ = 2.64 Wm−2) and increases by +1.5 Wm−2 per decade. Confidence in the observed LW‐down annual anomalies (Figure 5c) is boosted by the similarity among the individual stations’ time series. Common features include a network‐wide maximum in 1998, a sharp minimum in 2008, and a low‐amplitude, single‐mode sinusoidal variation between 1999 and 2007. Likewise, LW‐up annual anomalies (Figure 5d) show coherence among the various stations, especially in 1998, 2005 − 2007, 2008, and 2009 (discounting Desert Rock). The network mean LW‐up decreases by a token −0.94 Wm−2 per decade, but that tendency is also overwhelmed by its interannual variability (σ = 2.65 Wm−2). Because net‐LW is defined as (LW‐down − LW‐up), the combination of slightly increasing LW‐down (adding more LW to the surface with time) and decreasing LW‐up (removing less LW) explains why the surface net‐LW in Figure 4 increases. [22] Quantities that affect surface LW‐down irradiance are low‐level air temperature, low‐level atmospheric moisture, other greenhouse gases, clouds, and in some cases aerosols. According to NOAA's measurements, the global CO 2 concentration increased by about 30 ppm between 1995 and 2010 (ftp://ftp.cmdl.noaa.gov/ccg/co2/trends/co2_mm_gl.txt). Using the Fu‐Liou radiative transfer model [Fu and Liou, 1992], we determined that such a change in CO 2 concentration would result in a +0.14 Wm−2 increase in LW‐down at the surface, or +0.09 Wm−2 per decade, which is negligible compared with the interannual variability and 16‐year increase of LW‐down that was observed. The time series of the surface air temperature annual anomaly (Figure 7a) also exhibits a high degree of interannual variability, but the linear least squares fit to that series is flat, in agreement with reports that the mean global temperature stagnated during our analysis period [e.g., Foster and Rahmstorf, 2011; Hunt, 2011]. Two features dominate the surface temperature time series, a large coherent network‐wide peak in 1998, which was an El Niño year and the warmest year of the past century (http://lwf.ncdc.noaa.gov/oa/climate/research/1998/ann/global_temperatures.html), and a large amplitude warm − cold couplet from 2004 to 2011. The linear least squares fit to the network mean surface specific humidity annual anomaly (Figure 7b) decreases by about −0.3 g/kg over the 16‐year period but also exhibits a high degree of variability. With cloud coverage and the surface specific humidity decreasing, LW‐down would be expected to decrease, yet Figure 5c shows it increasing. A systematic lowering of boundary layer cloud base heights would also cause LW‐down to increase; however, the measured decrease in surface moisture over the analysis period should have led to higher cloud bases with time and, accordingly, a decrease in LW‐down. Therefore, a systematic change in cloud base height could not explain our measured increasing tendency of LW‐down. Figure 7 Open in figure viewer PowerPoint Annual anomalies of surface air temperature (a) and surface specific humidity (b) for individual stations (colored curves) and the network mean (black solid curve) in time series. The black dashed curve in each panel represents the annual average of the ONI ENSO index lagging by 1 year. The straight dashed line in each panel is the least squares linear fit to the network mean time series. [23] LW‐up responds strongly to the local surface air temperature. A plot of annual anomalies of surface air temperature vs. LW‐up (not shown) is linear and displays little scatter. The least squares linear fit has a correlation coefficient of 0.88, indicating that the surface air temperature explains 77% of the variance of LW‐up on the annual time scale. This indicates that the surface air temperature is governed locally by surface heating from which LW‐up derives the bulk of its signal.

5 Discussion [24] The similarity of features among the individual stations’ LW time series lends confidence that the data presented truly describe the actual variability and that generally all stations sense the same large‐scale phenomena. Given the high degree of interannual variability of LW‐down, its increasing tendency is likely the result of the way in which the parameters that influence its signal, i.e., temperature, clouds, and atmospheric moisture, interact in individual years. Foster and Rahmstorf [2011] have shown that, from 1979 to 2010, the short‐term variability of the global temperature responds to the El Niño/Southern Oscillation (ENSO), large volcanic eruptions, and solar variation; among these, ENSO is the most influential. Accordingly, given that there were no major volcanic eruptions between 1996 and 2011 and that the effect of the steady increase of CO 2 on LW‐down during our analysis period was negligible, ENSO would be expected to have the most significant influence on the variability of LW‐down. The time series of annual averages of the Oceanic Niño Index (ONI), an ENSO indicator, have been plotted along with the temperature and specific humidity time series in Figure 7a,b, but with a 1 year lag because the global response typically lags the ONI by about 4 months [Hansen et al., 2010]. The ONI is defined as the running 3 month sea surface temperature anomaly in the region between 5°N and 5°S and 120°W and 170°W and is used by NOAA to define ENSO events. It is positive for El Niño and negative for La Niña. It is clear that the ENSO index tends to follow the short‐term subdecadal variability of surface temperature and specific humidity. Linear regressions of surface temperature and specific humidity against the 1 year lagging ONI show correlation coefficients of 0.58 and 0.5, respectively. The ONI vs. virtual temperature, which combines temperature and moisture, has a correlation coefficient of 0.61, which specifies that ENSO explains 37% of the variance in the surface virtual temperature. Furthermore, the least squares linear regression of the ONI index vs. LW‐down (Figure 8) has a correlation coefficient of 0.7, indicating that ENSO explains 50% of the variance of LW‐down. Thus, the anomalous circulations associated with ENSO likely are responsible for a large part of the LW variability over the United States that we have documented. Figure 8 Open in figure viewer PowerPoint One‐year lagging ENSO ONI index vs. the network average LW‐down annual anomaly. The linear least squares fit (solid line) and its best‐fit equation and coefficient of determination (R2) are shown. [25] The +13 Wm−2 increase in total surface net radiation over our 16‐year analysis period represents surface absorbed energy that should have stimulated surface energy fluxes. That the surface specific humidity decreases over the same period indicates that the excess SRB energy would have gone mostly to sensible heat. If so, it is puzzling that the tendency of the surface air temperature is flat. Unfortunately, we do not have collocated surface energy flux measurements and thus cannot specify the actual behavior of latent and sensible heat within our analysis period. Perhaps the anomalous convective and advective circulations associated with ENSO events were responsible for the dissipation or dispersion of that excess energy in a way that had a net neutral effect on the air temperature.

6 Summary and Conclusions [26] A high‐quality, accurate, and comprehensive surface net radiation data set from the seven U.S. SURFRAD stations from 1996 to 2011 has been analyzed. The data set was reduced from minute‐scale measurements to monthly averages, from which annual averages and anomalies were generated. Most of the station‐month averages were represented by well over 90% of the possible data. Radiation data from the few station months with less than 70% data coverage were discarded and filled by using physical or empirical relationships. Missing atmospheric state variables were filled from collocated measurements or NCEP Reanalysis data. [27] The most significant result was that the network mean total surface net radiation increased at a statistically significant rate of +8.2 Wm−2 per decade, or nearly +13 Wm−2 over the 16‐year period, and that surface net‐SW accounted for 71% of that increase. The basic net‐SW and net‐LW components of the total surface net radiation were shown to be naturally anticorrelated and dominated by large swings of interannual variability. The trend of the surface net‐SW network mean was statistically significant and dominated by an increase in SW‐down of +6.6 Wm−2 per decade. We relied on a previous study by Long et al. [2009] and our own extension of that analysis to attribute much of that SW brightening to a decrease in cloud coverage. An analysis of 500‐nm AOD over the SURFRAD network from 1997 through 2009 showed AOD decreasing at a rate of ‐0.016 per decade, which extrapolates to ‐0.025 over the full 16 years. That decrease contributed only a negligible +0.82 Wm−2 to the SW brightening signal. A simultaneous increase in stratospheric aerosols from minor volcanic eruptions in the first decade of the present century minimally diminished SW‐down by −0.1 Wm−2 per decade, but that small decrease was included in our +0.82 Wm−2 result. A slight increase in SW‐up was inherently attributed to the increase in SW‐down and also to anomalous snow cover in the latter part of the period. [28] An increase of surface net‐LW from 1996 through 2011 accounted for 29% of the documented increase of total surface net radiation. That increase was attributed to an increase of LW‐down of +1.5 Wm−2 per decade and a decrease of LW‐up of ‐0.9 Wm−2 per decade, neither of which was statistically significant. Network mean annual anomalies of the primary factors that affect LW‐down, i.e., cloud coverage, air temperature, and specific humidity, varied individually in a way that intuitively would cause LW‐down to decrease over time, yet LW‐down increased slightly. The steady systematic increase of CO 2 by 30 ppm in the global atmosphere over the analysis period accounted for an increase in LW‐down of +0.14 Wm−2, which is an order of magnitude smaller than the changes and variability in LW‐down that were observed. Mean annual LW‐up and surface air temperature were strongly correlated, so our observed lack of a trend in the network mean annual air temperature is consistent with a relatively small, insignificant change in LW‐up over the 16‐year period. However, the hiatus in temperature is not consistent with the large, statistically significant increase in total surface net radiation of +8.2 Wm−2 per decade. [29] It is well known that the effect of the long‐term, steady increase of greenhouse gas concentrations will be manifested by changes in the SRB, especially its LW‐down component. However, it is apparent that the large interannual variability of LW‐down during our 16‐year analysis period was more of a response to ENSO that dwarfed the more subtle effects associated with systematically increasing greenhouse gases. Furthermore, it is possible that anomalous circulations associated with ENSO dissipated the excess sensible heat from the significant increase in total surface net radiation that was observed. It is clear that high‐quality SRB measurements must continue over many more decades to allow separation of the anthropogenic climate change signal from the natural variability. In addition, it would be beneficial to add surface heat flux measurements to SRB stations to elucidate the partitioning of the surface net radiation.

Acknowledgments [30] The quality, continuity, and length of the SURFRAD data set are due in large part to the contributions and dedication of Gary Hodges and Christopher Cornwall. We gratefully acknowledge the many positive contributions of Charles Long to the SURFRAD program. This analysis would not have been possible without the foresight and guidance of John DeLuisi, who proposed and directed the establishment of the SURFRAD network. Expert instrument calibrations provided by the National Renewable Energy Laboratory in Golden, Colorado; the WMO Region 4 Regional Solar Calibration Center/NOAA in Boulder, Colorado; and the World Radiation Centre in Davos, Switzerland, are duly credited for the accuracy and precision of the SURFRAD radiation measurements. Finally, expert comments by Martin Wild and the formal reviewers were beneficial and, in the opinion of the authors, improved the quality of this article.