Changes in the phase of the annual cycle of surface temperature

The annual cycle in the Earth’s surface temperature is extremely large—comparable in magnitude to the glacial–interglacial cycles over most of the planet. Trends in the phase and the amplitude of the annual cycle have been observed, but the causes and significance of these changes remain poorly understood—in part because we lack an understanding of the natural variability. Here we show that the phase of the annual cycle of surface temperature over extratropical land shifted towards earlier seasons by 1.7 days between 1954 and 2007; this change is highly anomalous with respect to earlier variations, which we interpret as being indicative of the natural range. Significant changes in the amplitude of the annual cycle are also observed between 1954 and 2007. These shifts in the annual cycles appear to be related, in part, to changes in the northern annular mode of climate variability, although the land phase shift is significantly larger than that predicted by trends in the northern annular mode alone. Few of the climate models presented by the Intergovernmental Panel on Climate Change reproduce the observed decrease in amplitude and none reproduce the shift towards earlier seasons.


ARTICLES
The annual cycle in the Earth's surface temperature is extremely large-comparable in magnitude to the glacial-interglacial cycles over most of the planet. Trends in the phase and the amplitude of the annual cycle have been observed, but the causes and significance of these changes remain poorly understood-in part because we lack an understanding of the natural variability. Here we show that the phase of the annual cycle of surface temperature over extratropical land shifted towards earlier seasons by 1.7 days between 1954 and 2007; this change is highly anomalous with respect to earlier variations, which we interpret as being indicative of the natural range. Significant changes in the amplitude of the annual cycle are also observed between 1954 and 2007. These shifts in the annual cycles appear to be related, in part, to changes in the northern annular mode of climate variability, although the land phase shift is significantly larger than that predicted by trends in the northern annular mode alone. Few of the climate models presented by the Intergovernmental Panel on Climate Change reproduce the observed decrease in amplitude and none reproduce the shift towards earlier seasons.
Climate change is often described by trends in annual mean temperature, but large seasonal temperature changes exist independent of changes in the annual mean. A small literature exists concerning the variability in the phase of the annual cycle. Thomson 1 examined the Central England Temperature time series , and identified a trend in the phase of the annual cycle towards later seasons, beginning around 1950, that is anomalously large in the context of the preceding several-hundred-year record. He argued that this excursion is associated with increases in atmospheric CO 2 concentration. He also presented evidence of trends in the phase of the annual cycle over larger spatial scales and an increase in the spatial variance of the trends. Mann and Park 2 and Wallace and Osborn 3 demonstrated that the hemispheric averaged observations contain trends in amplitude and phase. The amplitude trend is negative and is related to the observation that winter is, on average, warming more quickly than summer [4][5][6] . The hemispheric phase trend, however, is towards earlier seasons, opposite in direction to that found by Thomson 1 for central England.
The importance of these observed amplitude and phase trends is hard to judge because we lack a good model for natural variability. Wallace and Osborn 3 used two criteria to evaluate whether the observed trends are unusual: (1) a statistical test for the presence of a trend and (2) a comparison of trends with natural variability as represented in a general circulation model. Neither of these approaches is altogether satisfactory. We expect low-frequency variability always to be present, so the presence of a trend in and of itself is not surprising 7 . Furthermore, general circulation models may not give us an accurate picture of low-frequency variability, particularly in phase, because of two shortcomings. First, the models that have been used to evaluate phase and amplitude variability have used seasonal heat and freshwater flux adjustments to match the mean annual cycle, which may artificially stabilize the modelled annual cycle. Second, and more troubling, Northern Hemisphere phase trends predicted by models forced with twentieth-century anthropogenic forcing are in the opposite direction to the observed trend 2,3 . Modelled Northern Hemisphere amplitude trends also disagree with observations when compared using a temporally fixed network 3 .
Models that are unable to replicate observed trends are clearly not ideal for constraining the range of natural variability. Instead, we appeal to the early observational record to estimate the natural spatial and temporal variability of the seasonal cycle and ask if the trends observed in the recent record are anomalous in nature.
The basic state of the annual cycle Two distinct temperature-based methods for discussing the timing of the seasons have been used in the literature. The more common is a threshold-based model wherein seasonal transitions are defined as the times of year when the temperature rises above or drops below some specific value. In this framework, the 'spring' threshold will be reached earlier if temperature increases uniformly through the year. This type of change is of first-order importance for explaining changes in seasonality observed both in biological systems (for example flowering dates 8,9 , bird migration timing 8,9 and terrestrial surface carbon uptake 10 ) and in components of climate that exhibit threshold responses (for example the freezing and melting of ice 11 ). However, threshold-based definitions conflate changes in the phase of the annual cycle with changes in the annual mean (see Supplementary  Information). We instead describe the seasonal cycle by the amplitude and phase of the yearly-period sinusoidal component in surface temperature, a measure of seasonality that is distinct from changes in the annual mean [12][13][14][15][16] , and reference it to the yearly-period sinusoidal component in local solar insolation. The difference between the temperature and local insolation phases (l 5 w T 2 w sun ) is the lag 17 , and the ratio of the amplitudes (G 5 A T /A sun ) is the gain (see Methods Summary). We examine gridded 5u 3 5u temperature records from the University of East Anglia's Climate Research Unit 18,19 and analyse long-term mean, detrended variability and trend fields for land (l land , G land ) and ocean (l ocean , G ocean ).
The spatial patterns of l and G (Fig. 1a, b) are dominated by the contrast between land and ocean. The larger ocean thermal mass causes it to respond more sluggishly to oscillatory forcing than land, which results in a smaller and later oceanic response. Ocean points have a mean gain of 28 uC (kW m 22 ) 21 (standard deviation of point-wise means, s m 5 15 uC (kW m 22 ) 21 ) and a mean lag of 56 days (s m 5 11 days); in comparison, the more rapidly adjusting land has a mean gain of 74 uC (kW m 22 ) 21 (s m 5 23 uC (kW m 22 ) 21 ) and a mean lag of 29 days (s m 5 6 days).
Superimposed on the dominant land-ocean contrast is an eastwest gradient in G and l. As we move from west to east (following the prevailing winds) across the mid-latitude continents, there is a tendency towards a more rapid response (large G land , small l land ). This cross-continent gradient in G land is quite strong, whereas the gradient in l land is relatively weak (l land adjusts rapidly to interior values along the western continental margin 1,2 ). Conversely, as we move from west to east across the mid-latitude ocean basins, there is a tendency towards a more sluggish response (small G ocean , large l ocean ), and the relative strengths of the G ocean and l ocean gradients is reversed relative to that of the land (G ocean adjusts rapidly to interior values along the western margin of ocean basins).
The role of land-sea contrast in setting the climatological distribution of the annual cycle is not a new observation 17,[20][21][22] , but its dominance is particularly obvious when considering the relationship between G and l. Pairs of G and l fall along an arc (Fig. 2a). We define a 'seasonal response index' to represent a point's position in this laggain space as and find that 75% of the variance in this index is explained by the distance between a grid point and the coast to its west, where distance is taken as positive for land and negative for ocean (see Supplementary Information for more discussion on the structure of variability). The relationship between SRI and distance from the coast holds best in Eurasia, where southern mountains constrain the transport to be zonal and isolate the interior from tropical moisture. Deviations from this general east-west pattern are found in regions where transport is less zonal, such as the southeastern North American monsoonal region, where there is strong poleward moisture transport onto land, and in the western United States, where the north-south alignment of the Rocky Mountains effectively blocks oceanic influence from the Pacific Ocean.
The observed arc in the relationship between G and l is a ubiquitous feature of seasonally driven models that contain interacting land and ocean regions, and can be understood as the natural consequence of interacting sinusoids. Consider two sine waves with different phases and amplitudes, S 1 5 Asin(w) and S 2 5 (A/r)sin(w 1 Dw). A weighted average of the two sine waves, wS 1 1 (1 2 w)S 2 , with 0 , w , 1, yields a sine wave with amplitude A w~A r ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi w 2 r 2 z2rw(1{w) cos Dwz(1{w) 2 q and phase A mixing line using this equation (Fig. 2a) is consistent with the general form of the observed arc. Apparently, the spatial structure associated with the seasonal cycle can be understood, to first order, as the result of variable mixing between continental and marine influence. Variability in G (Fig. 1d) tends to be largest where the climatological G is large (coefficient of correlation, R 5 0.83), and is about twice as large over land (mean of point-wise standard deviations, s s 5 5.2 uC (kW m 22 ) 21 ) as it is over the ocean ( s s 5 2.5 uC (kW m 22 ) 21 ). Conversely, temporal variability in l (Fig. 1c) is correlated with G 21 (R 5 0.62) and is larger over the ocean ( s s 5 5.0 days) than it is over land ( s s 5 4.0 days). The inverse relationship and larger l ocean variability arises because finite perturbations more readily alter the phase of a smaller amplitude sinusoid (see the Supplementary Information discussion on natural variability). We thus expect that it will be more difficult to detect the presence of any true phase trend over the ocean.

Trends in the phase and gain of the annual cycle
The 1954-2007 l land trends ( Fig. 1e) are predominantly towards earlier seasons, with a mean decrease of 1.7 days (that is, 6%) over the past 54 years. The l ocean trends are large but regionally disparate. For example, the interior of the North Pacific, and the Atlantic north of 50u N, exhibit trends towards later seasons, whereas along the eastern edge of the North Pacific, and in the North Atlantic south of 50u N, trends are primarily towards earlier seasons. The mean l ocean shift is towards later seasons by 1.0 days over the past 54 years.
A comparison of trend maps (Fig. 1e, f) and variability maps (Fig. 1c, d) reveals that the trends are large where the detrended variability is large. This suggests the obvious null hypothesis that the trends are merely a manifestation of natural variability. One test for whether the trends observed in the recent record are consistent with natural variability is to compare them with trends observed in earlier periods. We consider the distribution of point-wise trends ( Fig. 2b) for the 1900-1953 and 1954-2007 intervals using those records which have good temporal coverage during both intervals (see Methods; this is the default distribution of records that we use below, unless specifically stated otherwise). Land and ocean are considered separately because the characters of their annual cycles are so different. We plot long-term-mean value (a, b), temporal standard deviation of the detrended time series (c, d); and trend in days per 54 years (e) and uC (kW m 22 ) 21 per 54 years (f). Both variability and trend maps are plotted on the 'dense network' , without land and ocean masks applied. Results have been excluded in the tropics, where data availability is poor, and where less than 85% of the variance in an average year is explained by the yearly component.
We adopt a null hypothesis that the mean of each distribution of trends is zero. Testing this null hypothesis requires a knowledge of the effective spatial degrees of freedom 23,24 , and we use estimates obtained from the moment-matching method of ref. 25 (see Methods). Of the four distributions of l trends that we consider, only those over land during the 1954-2007 interval have a mean that differs significantly from zero (Table 1), and here the significance is marked (21.9 6 1.0 days per 54 years, P , 0.001). Repeating the tests for 1954-2007, using the larger spatial network that is available for this interval (the dense network; see Methods), supports the significance of the l land trend (21.7 6 0.8 days per 54 years, P , 0.001). We also detect a significant 1954-2007 l ocean trend towards later seasons (1.0 6 0.9 days per 54 years, P 5 0.02) that is only detectable in the more extensive dense network.
The dominant signal in the 1954-2007 G trend (Fig. 1f) is a decrease in the amplitude of the annual cycle over land, averaging 22.5 uC (kW m 22 ) 21 over the past 54 years on the dense network, a 3% drop. This is the well-known amplification of winter warming 4,5,26 , which is strongest in the interior of Eurasia and in the boreal forests of western Canada. Note, however, that large regions exist where the amplitude has increased. In western Europe and the Middle East, the observed increase in G is associated with greater warming in summer than in winter. In the central North Pacific and the southeastern United States, the increase in G results from winter cooling. In Quebec, the summer warming and winter cooling trends are of comparable magnitudes, leaving little trend in mean temperature but a measurable increase in G. Ocean G trends are almost everywhere small and show a mean increase of 0.4 uC (kW m 22 ) 21 over the past 54 years.
We apply tests to the G trends similar to those made on the l trends (Fig. 2c). Of the four distributions of G trends, only those over land during 1954-2007 have a mean that differs significantly from zero (22.6 6 2.4 uC (kW m 22 ) 21 per 54 years, P , 0.03; Table 1). If the dense network is used instead for this interval, the G land trends remain significant (P 5 0.05) and the G ocean trends emerge as being marginally significant (P 5 0.07). It is noteworthy that the well-reported changes Westward distance from coast (10 3 km) Relative frequency a b d c Land, 1900-1953Land, 1954Land, 1954, dense Ocean, 1900-1953Ocean, 1954Ocean, 1954 a, Observed relationship between local gain, G, and lag, l, for Northern Hemisphere extratropical locations. Colour represents the distance between a grid point and the coast to its west (positive for land, negative for ocean). Outliers with l , 20 days are from the Indian subcontinent and presumably reflect monsoon dynamics. The black line shows the nonlinear relationship between amplitude and phase for weighted averages of two end-member sine waves (see Supplementary Information)   in the amplitude of the annual cycle 19,27 are less significant than the less-reported land-phase trend. The low significance of the amplitude results is related to the large natural variability in wintertime temperature. Winter warming is considerably stronger than summer warming over land during 1954-2007, but the variance in winter land temperatures is about four times that in summer land temperatures, making the winter trend less significant and making detection of changes in amplitude difficult 28 . In fact, we are unable to detect a significant difference between summer and winter warming when temperature trends are analysed as the difference between threemonth seasonal averages (Table 1). We focus on the l land trends because their significance is markedly higher than that of any other observed trend. Furthermore, G trends are more easily discussed in the seasonal-average-temperature framework than are l trends and have received more attention elsewhere 5,6,27 .
There are two steps in establishing the presence of an anomalous trend. The first is establishing that a trend is statistically distinguishable from zero, which we demonstrated for the 1954-2007 l land observations. The second is establishing that this trend is different in character from what we would expect in the naturally varying system, which is more difficult given the finite length of the instrumental temperature records. We are particularly concerned about low-frequency variability being misinterpreted as an anomalous trend 7 . The absence of a significant l land trend for the 1900-1953 test period indicates that the trend during 1954-2007 is anomalous. By restricting ourselves to a smaller set of locations, we can also extend our analysis back to 1850. We construct an average l land time series by averaging the phase time series from all of the land grid boxes with perfect temporal coverage between 1850 and 2007, and adopt the null hypothesis that the 1954-2007 trends result from natural lowfrequency variability as represented in the 1850-1953 record. We build a distribution for this null hypothesis by calculating the trends of many synthetic time series having the same spectral amplitude structure as the 1850-1953 record, but with randomized phases 29 (see Supplementary Methods), and are able to reject it with very high confidence (P 5 0.006). The phase trend over the past 54 years is not consistent with the structure of natural variability found in the earlier record. Furthermore, there is no 54-year period in the 1850-1953 control period that would allow rejection of this null hypothesis. (Note that we are unable to meaningfully compare the 1954-2007 trends in l ocean with the 1850-1953 period because instrumental coverage over the ocean during these early times was poor.) Finally, it is reasonable to ask whether the observed variability is really best thought of as a shift in the yearly sinusoidal component, or if it would be better described by changes in individual months. A change in a single month's temperature will map into a shift in the annual cycle, although the yearly frequency component provides a poor description of such an anomaly. We calculate the mean annual cycles for four 27-year periods, using land grid boxes with good temporal coverage between 1900 and 2007 (see Supplementary Methods), and consider their anomalies from the 108-year mean annual cycle (Fig. 2d). (Consideration of the means in these four periods gives insight into the origins of the trends in the 1900-1953 and 1954-2007 intervals.) The most recent anomaly time series (1981-2007) exhibits the largest departures from the long-term mean, and 80% of its variance is explained by the yearly component. The most recent period has more variability at the annual period than the total variability during all preceding periods, highlighting both that these shifts are well described by a yearly sinusoidal component and the anomalous nature of the recent changes (see Supplementary  Table 2).

Origins of the changes in the annual cycle
To explore the origins of the shifts in G land and l land , we first turn to the global climate model results summarized in the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) 30 . In particular, we analyse the 72 simulations of twentieth-century climate that use observed forcings conducted for the World Climate Research Programme's (WCRP) Coupled Model Intercomparison Project Phase 3 (ref. 31). The distributions of mean trends in G land and l land found in the models easily encompass the observed trends during the 1900-1953 interval. However, during the 1954-2000 interval (the simulations stop at the end of the century), the observed decrease in G land has a larger magnitude than all but six of the model simulations, and no model reproduces the observed shift towards earlier seasons (Fig. 3). The mean of the model l land trends for 1954-2000 is towards later, not earlier, seasons. Furthermore, 25 atmospheric model runs forced with observed sea surface temperatures for the 1978-2000 period (the Atmospheric Model Intercomparison Project models 32 ) do no better at replicating the observed l land trends (see Supplementary Methods and Discussion).
The IPCC model results do not appear to give us an explanation of the observed trends, except to suggest that the answer involves something that the models do not capture. We thus retreat to a simple, conceptual model to explore how local processes may cause variability in l and G consistent with the observations. Our goal here is to explore some obvious candidates and roughly estimate the size of the perturbation needed to explain the observation. We use a one-box energy Many of the mechanisms invoked to explain variability in annual mean temperature are unlikely to be directly responsible for the observed shift in phase. Doubling the atmospheric long-wave optical depth to simulate the radiative effect of a very large increase in greenhouse gases has essentially no effect on seasonal timing (Dl 5 0.1 days, DG 5 20.6 uC (kW m 22 ) 21 ). Increasing solar luminosity by a fixed percentage increases the amplitude of the temperature response by the same percentage and has a negligible effect on phase 1,35 . Decreases in sea ice (not represented in our model) present the atmosphere with a larger thermal mass, implying a delayed seasonal response (although threshold responses at the time of spring melt may induce changes in the opposite direction). Consistent with this intuition, the (incorrect) phase delays found in the model results of ref. 2 are attributed to decreases in sea-ice cover.
However, there exist numerous mechanisms that may shift the seasonal cycle in the observed direction. A decrease in thermal mass on land of 8 6 4% is sufficient to produce the observed offset in phase of 1.7 6 0.8 days. Thermal mass on land is largely modulated by soil moisture. Compare, for example, the effective thermal mass of a dry desert sand (1.9 J m 23 uC 21 ) with that of a saturated loam soil (3.2 J m 23 uC 21 ). For a typical soil, the observed phase shift would require a 13 6 7% decrease in soil moisture. IPCC Fourth Assessment Report model runs disagree with each other on the sign of recent soil moisture trends and show little skill at explaining the (sparse) observations 36 . The few available long-term measurements suggest increased soil moisture over the latter part of the twentieth century 37,38 , which is inconsistent with the thermal mass hypothesis, although we observe that drought reconstructions 39 indicate these observations may not be representative of continental-scale variations. The paucity of records with more than 40 years of data prohibits a more detailed comparison. We consider large-scale decreases in soil moisture to be a viable candidate for inducing the observed shift towards earlier seasons.
Perturbations to atmospheric short-wave optical properties are also effective at modifying the annual cycle, and it appears that short-wave absorptivity has been changing, perhaps because of aerosols [40][41][42] . The Earth's short-wave optical properties are not constant throughout the year, and their potential range of variability is not captured by this simple model. Nonetheless, the model indicates that variability in atmospheric annual mean reflectivity, absorptivity or transmissivity on the order of 10% will change l land by the observed amount. Note that Wallace and Osborn 3 were unable to replicate the observed hemispheric phase shifts using a general circulation model, but that the inclusion of aerosol forcing did decrease the modelled (incorrect) shift towards later seasons. We see no indication of shifts in mean l land after any of the major volcanic eruptions of the past century, although the effects of stratospheric and tropospheric aerosols on phase are likely to be quite different.
Thomson 1 makes the intriguing, although difficult-to-evaluate, proposal that decreases in phase are due to an increased local sensitivity to anomalistic year forcing (associated with the annual cycle in Earth-Sun distance) relative to tropical year forcing (due to the annual cycle in the orientation of the Earth's rotation axis relative to the Sun).
The above-mentioned changes in albedo, soil moisture and shortwave forcing have all been implicated in changing modes of atmospheric circulation [43][44][45] . This raises the further possibility that shifts in atmospheric circulation participate in the modulation of the annual cycle. We focus on the northern annular mode (NAM) and the Pacific North American pattern, as these have been shown to represent the bulk of the variability in standard atmospheric climate indices 46 (but see Supplementary Information for a more complete analysis). The NAM shows significant cross-correlations with time series of 1950-2007 spatially averaged l land (R 5 20.5, P , 0.001) and G land (R 5 0.42, P 5 0.007), whereas the Pacific North American pattern has significant correlation with G ocean (R 5 0.3, P 5 0.04). Apparently, atmospheric dynamical processes respond to similar forcing mechanisms as l land or themselves participate in altering l land through the advection of heat and moisture or other indirect processes. Northern Hemisphere snow cover, for example, is known to interact with the NAM 43,47 , and winddriven changes in mixed layer depth affect the thermal mass that the ocean presents to the atmosphere 48 .
The Atmospheric Model Intercomparison Project simulations have been found to capture the spatial pattern, but not the temporal pattern, of NAM variability 49 , just as we find that the models fail to capture the long-term trends in phase. However, the recent phase excursion appears to be only partly explained by the late-twentiethcentury excursion in the NAM. A regression of the spatial average l land time series against the NAM index from 1950-2007 removes 25% of the variance and 40% of the trend in the l land time series, but still leaves a significant trend (P , 0.02) of 21.0 days per 57 years (see Supplementary Information).
The statistics of the distribution of l land trends are well described as natural variability from 1900-1953, but the distribution shifts in 1954-2007, the period in which anthropogenic interference with mean temperature becomes apparent. If we extend our natural control period back to 1850, the recent trends appear yet more anomalous. Numerous climate factors can influence the phase of the annual cycle, and it appears that some portion of the trend in the annual cycle is associated with changes in the NAM during the late twentieth century. We expect that a complete explanation for long-term shifts in atmospheric circulation will also encompass an explanation of the variability in the phase of the annual cycle. Although the mechanism is still uncertain, the tests we apply to the 1954-2007 trends in land phase indicate that they are inconsistent with natural variability, and thus appear to be due to anthropogenic influence.

METHODS SUMMARY
For each year of data, we calculate the yearly (one cycle per year) sinusoidal component using the Fourier transform, as where x(t 1 t 0 ), t 5 0.5, …, 11.5, are 12 monthly values of either the de-meaned monthly temperature or de-meaned monthly insolation and t is time in months. The factor of two is to account for both positive and negative frequencies. Phase is given by w x 5 tan 21 (Im(Y x )/Re(Y x )). To discuss both hemispheres in a common framework, we reference the temperature phase, w T , to the local solar insolation phase, w sun . The difference between these two phases is the lag 17 , l 5 w T 2 w sun . Amplitude is given by A x 5 jY x j. For the purpose of understanding the response of the Earth's temperature to forcing, we examine the gain, which is defined as the ratio of the amplitudes of the annual cycles in temperature and insolation, G 5 A T /A sun . Unlike insolation or temperature amplitudes alone, G has very little latitude dependence.
If any of the 12 monthly temperature values is missing in the data set at a given location, then no estimate of the annual cycle is made at that location for that year. Analyses using longer record pieces and more advanced filter techniques do not change our conclusions regarding the significance of phase and amplitude changes.
Full Methods and any associated references are available in the online version of the paper at www.nature.com/nature.

METHODS
Data sets. When plotting full spatial fields (Figs 1, 2a), we use the HadCRUT3 blended land-and-ocean 5u 3 5u gridded surface temperature anomalies 18 plus gridded climatology 19 from the Climate Research Unit at the University of East Anglia. For comparison with IPCC model archive output (Fig. 3), we use the HadCRUT3 data set (because models do not calculate separate land and ocean surface temperatures), but restrict ourselves to grid points that are more than 50% land. When land and ocean are considered separately (all other calculations), we instead use the CRUTEM3 (for land) and HadSST2 (for ocean) data sets so that the land and ocean signals are more cleanly separated. Dominantly land and ocean boxes are identified using the Clark US Navy Fleet Numerical Oceanographic Center Land/Ocean Mask 51 . Because a grid box with land cover of only a few per cent is not as representative of land as a continental interior box, only CRU grid points that are more than 50% land in the Navy mask are used in land calculations. Those that are less than 50% land are considered ocean. Data mask. We desire a high ratio of signal (annual amplitude) to noise (errors in observations and high frequency temperature variability), so that we can isolate variability that is associated with changes in the annual cycle. In fact, the yearly sinusoidal component dominates the extratropical records; on average, it explains 96% and 90% of the within-year variance in monthly temperatures for land and ocean grid boxes, respectively. With the time series of yearly G and l at each grid box, we estimate the long-term means, the long-term trends and the standard deviations of the departures from the long-term trend. To do so, we exclude from analysis (1) those extratropical grid boxes where less than 85% of the average within-year variance is explained by the yearly sinusoidal component (primarily the Southern Ocean) and (2) all tropical grid boxes (23.5u S to 23.5u N), because the two-cycles-per-year harmonic in forcing and response is strong in this region. For calculating long-term-mean l and G, we exclude grid boxes with fewer than ten yearly estimates over the entire record. For calculating 54-year trends and detrended standard deviation in G and l, we exclude grid boxes with fewer than 40 yearly estimates. Trends are calculated using a least-squares fit. For comparing 1900-1953 and 1954-2007 trends, we use a 'comparison network' of grid boxes that meet these data-density criteria for both periods (180 land points, 120 ocean grid points). We also use a 'dense network' of all of the grid boxes that meet the data inclusion criteria for 1954-2007, to obtain a best estimate for the most recent period (299 land points, 345 ocean points). The dense network has good spatial coverage between 25u N and 60u N (with some missing values in the interior of Eurasia and at higher latitudes) and more sporadic coverage between 25u S and 40u S. For the comparison network, all of the Southern Ocean, most of the Pacific Ocean and much of Asia are excluded. Trend distribution testing. Tests for the deviation of distribution means from zero are done using the t-test (two-tailed) and confidence intervals are t-intervals. The standard deviation for the t-test is calculated from the observed distribution and the degrees of freedom are estimated as the effective spatial degrees of freedom (ESDOF) of the time-varying field using the moment-matching method of ref. 25, which they describe as appropriate when testing for the difference of a realization from the mean. This method estimates 21 (l land ), 19 (l ocean ), 12 (G land ) and 20 (G ocean ) ESDOF for l and G variability. For the late dense network, the estimates are 29 (l land ), 58 (l ocean ), 12 (G land ) and 47 (G ocean ) ESDOF. For summer temperature field variability we use 15 (land) and 9 (ocean) ESDOF for the comparison network and 17 (land) and 30 (ocean) ESDOF for the dense network. For winter temperature field variability we use 8 (land) and 9 (ocean) ESDOF for the comparison network and 6 (land) and 31 (ocean) ESDOF for the dense network. For seasonal-difference (summer temperature minus winter temperature) hypothesis testing we use ESDOF values calculated from fields of annual mean temperature and we use values of 10 (land) and 9 (ocean) ESDOF for the comparison network and 12 (land) and 22 (ocean) ESDOF for the dense network. Recovered ESDOF estimates are comparable to the observation-based estimates of ref. 23 and are notably smaller than the model-based estimate of ref. 24. For testing the average summer and average winter trend distributions, summer is defined as June, July and August in the Northern Hemisphere and as December, January and February in the Southern Hemisphere. Winter is defined as December, January and February in the Northern Hemisphere and June, July and August in the Southern Hemisphere. Energy balance model. The one-box conceptual model is a one-atmosphericlayer energy balance model, with a black-body surface and a black-body atmosphere, forced with sinusoidally varying short-wave radiation (S 5 S 0 cos(2pt)) characteristic of the annual cycle in radiation at 40u N. We add two complications: (1) to consider sensitivity to atmospheric optical properties, we specify atmospheric short-wave absorptivity (A 5 0.15), transmissivity (T 5 0.6), and reflectivity (R 5 0.25), and calculate the effects of multiple reflections following ref. 34; (2) to consider the effects of increasing long-wave optical depth (t), we calculate a different atmospheric upward-radiating temperature (T up 5 T a 2 CH(ln(3t/ 2) 2 1)) and downward-radiating temperature (T down 5 T a 1 CH), which are related to the interior atmospheric temperature (T a ) by the atmospheric height (H) and lapse rate (C), following ref. 33 (ch. 5). Surface temperature (T s ) tendency is a function of the sum of energy fluxes divided by the thermal mass of the surface (c s ). On land, the depth of soil that contributes to the thermal inertia is estimated as the square root of the soil diffusivity times the timescale in question, and for annual timescales we use a depth of 4.7 m in calculating c s . where a s is the surface albedo. The model is run for parameter values typical for land, and we then perturb these values to estimate sensitivity. Thermal mass changes are equated with soil moisture changes assuming a soil consisting of 10% inorganic matter, 45% organic matter, 5% unfilled airspace and with a soil water content of 40%. A 13% drop in soil moisture then implies that the soil water content drops to 35%.