Accurate assessments of future climate impacts require realistic simulation of interannual–century-scale temperature and precipitation variability. Here, well-constrained paleoclimate data and the latest generation of Earth system model data are used to evaluate the magnitude and spatial consistency of climate variance distributions across interannual to centennial frequencies. It is found that temperature variance generally increases with time scale in patterns that are spatially consistent among models, especially over the mid- and high-latitude oceans. However, precipitation is similar to white noise across much of the globe. When Earth system model variance is compared to variance generated by simple autocorrelation, it is found that tropical temperature variability in Earth system models is difficult to distinguish from variability generated by simple autocorrelation. By contrast, both forced and unforced Earth system models produce variability distinct from a simple autoregressive process over most high-latitude oceans. This new analysis of tropical paleoclimate records suggests that low-frequency variance dominates the temperature spectrum across the tropical Pacific and Indian Oceans, but in many Earth system models, interannual variance dominates the simulated central and eastern tropical Pacific temperature spectrum, regardless of forcing. Tropical Pacific model spectra are compared to spectra from the instrumental record, but the short instrumental record likely cannot provide accurate multidecadal–centennial-scale variance estimates. In the coming decades, both forced and natural patterns of decade–century-scale variability will determine climate-related risks. Underestimating low-frequency temperature and precipitation variability may significantly alter our understanding of the projections of these climate impacts.
Many instrumental and paleoclimate records suggest that the spectrum of climate variations can be described by power-law scaling (e.g., Blender and Fraedrich 2003; Fraedrich and Blender 2003; Ditlevsen et al. 1996; Pelletier 1998; Pelletier and Turcotte 1997; Huybers and Curry 2006; Franke et al. 2013; Laepple and Huybers 2013, 2014a,b; Ault et al. 2013b, 2014). A power-law relationship implies that variations in climate (i.e., temperature or precipitation) exhibit higher (or lower) amplitudes at lower frequencies (Huybers and Curry 2006). However, comparative studies among climate model and observation data have indicated that there may be up to an order of magnitude mismatch between proxy and simulated temperature and precipitation variability at decade–century time scales (e.g., Ault et al. 2013a,b; Franke et al. 2013; Laepple and Huybers 2014b).
Understanding variability at these time scales is important for studies of climate change, including analyses of detection and attribution of climate change impacts. Additionally, understanding variability at these time scales is important for simulation of abrupt change (Valdes 2011; Shuman 2012), prolonged drought (Pelletier and Turcotte 1997; Ault et al. 2013b, 2014), and lower-frequency climate variations (e.g., Thornalley et al. 2009). If climate models underestimate variability on decadal–centennial time scales, they will underestimate the likelihood of abrupt change, megadrought, and other extreme climate events.
Beyond climate prediction, understanding the true shape of local precipitation and temperature spectra is integral to our understanding of climate itself. For instance, where is long-period variability based on coupled dynamics, and where does variability arise solely from autocorrelation and smoothing of weather noise? Manabe and Stouffer (1996) showed that dampened weather noise can explain much of the shorter-period temperature variability in a coupled atmosphere–ocean model. However, ocean circulation may lead to significant longer-period (centennial scale) regional temperature variability (Bjerknes 1964; Manabe and Stouffer 1996).
Here we expand the analysis of the continuum of climate variability to multiple models from phase 5 of the Coupled Model Intercomparison Project (CMIP5) and paleoclimate records. We test whether the latest Earth system models generate local low-frequency temperature and precipitation variability greater than variability expected from a simple autocorrelated process. The tendency for many climate variables to exhibit some degree of “redness” (i.e., increased variance at longer time scales) is well documented (Kutzbach and Bryson 1974; Blender and Fraedrich 2003; Fraedrich and Blender 2003; Ditlevsen et al. 1996; Pelletier 1998; Pelletier and Turcotte 1997; Huybers and Curry 2006; Franke et al. 2013; Laepple and Huybers 2013; Ault et al. 2013b). In its simplest form, this tendency may arise from the large thermal capacity of the oceans and cryosphere, which integrate high-frequency “white noise” forcings, including weather, interannual variations, and even short-lived cooling from volcanic eruptions (Hasselmann 1976; Newman et al. 2003; Newman et al. 2016; Sigl et al. 2015). We evaluate the distribution of temperature and precipitation variance simulated by coupled Earth system models against the null hypothesis that low-frequency variations are generated exclusively from a monthly lag-1 autocorrelated [AR(1)] process (e.g., ocean damping of weather noise). We also test whether pre-anthropogenic external forcings of the last millennium significantly influence modeled variance distributions.
2. Data and methods
a. CMIP5 simulations
We used gridded surface temperature (TS) and precipitation (PR) model output from 10 all-forcing CMIP5 last millennium (CMIP5 LM) and 10 all-forcing NCAR CESM1.1(CAM5) last millennium ensemble (CESM LME) simulations (Taylor et al. 2012; Schmidt et al. 2012; Otto-Bliesner et al. 2015). In addition to the all-forcing CESM LME simulations, we used one control run and 17 single-forcing simulations from the CESM LME, including volcanic (N = 5), solar (N = 3), greenhouse gas (N = 3), orbital (N = 3), and land-use and land-cover changes (N = 3) to examine the influence of external forcing on low-frequency climate variability. We also compared all available NASA GISS-E2-R last millenium (LM) simulations with and without volcanic forcing. We only compared last millennium forcing scenarios in CESM and GISS-E2-R because these were the only simulations available at the time of our analysis. The LM experiment typically spans the period 850–2005 Common Era (CE), but here we limited our analysis to model data covering the time period 850–1850 CE to exclude the influence of anthropogenic warming on long-period variability. We also used CMIP5 preindustrial control (CMIP5 PIC) simulations (Taylor et al. 2012) that span at least 500 yr to estimate unforced interannual–centennial-scale variability from the same modeling groups that conducted the LM experiment (Table 1). Finally, we used CMIP5 historical (CMIP5 HIST) simulations to estimate forced interannual–centennial-scale temperature variability during the time period 1850–2005 CE. Although we calculated β for all models and ensemble members listed in Table 1, we show only three representative ensemble member maps from the CESM LME in the main text to highlight our results. We also show results from the 10 CMIP5 LM members and the same 10 CMIP5 PIC models in the main text to highlight the forcing versus control differences; additional results are included in the supplemental material.
b. Power spectra estimation
We estimated power spectra from model output for each grid box on the model’s native grid structure to avoid the smoothing effects of regridding. We then calculated the power-law scaling parameter β [from power-law relationship for frequency f and variance S(f), where S(f) ∝ f−β], using the methodology of Huybers and Curry (2006) and Ault et al. (2013b). We first calculated the annual mean for each gridbox time series (Figs. 1a,b). We then standardized all temperature and precipitation (using a modified version of the inverse Rosenblatt transform; Van Albada and Robinson 2007) time series to exhibit unit variance. We next estimated the power spectrum for each nondetrended time series using Thomson’s multitaper method and specified 2 as the time–bandwidth product (Slepian sequences) used as data windows (Thomson 1982). We then binned the spectra in evenly spaced frequency bins (spaced at equal 0.2 log intervals) to avoid overemphasizing high-frequency variance in regression calculations. We calculated a least squares fit line for log-transformed spectral power and frequency between interannual- (1/2 yr) and centennial-scale (1/500 yr) frequencies. The slope of this line in log–log space equals the value of β that we mapped in all figures (Figs. 1c–f).
We refer to a time series as “red” if its spectrum has increasing variance at lower frequencies (a negative slope in log–log space), such as surface temperature in CESM in the North Atlantic (Fig. 1c). We refer to a time series as “blue” if its spectrum has decreasing variance at lower frequencies (a positive slope in log–log space), such as surface temperature in CESM in the tropical Pacific (Figs. 1d,f). Although certain regions appear blue on many of the maps of β, the spectral shape in these regions is generally not appropriately described by a linear fit owing to the large interannual variance peak (Figs. 1d,f). Therefore, although power-law behavior can be used to derive expected return times for extreme events (e.g., Roe and Baker 2016; Ault et al. 2013b), the apparent negative scaling in the tropics should not be used to calculate return times for extreme events owing to the lack of negative power-law scaling behavior in this region. We refer to a given time series as “white” if it roughly had equal variance at all frequencies and had a flat slope (between −0.2 and 0.2) (Ault et al. 2013b) in log–log space, such as precipitation over the North Atlantic in CESM (Fig. 1e).
After power spectra estimation and slope β calculation for each native grid box, we then interpolated all data onto a ~2° horizontal grid for plotting and averaging among different models. If multiple runs from the same model were available, we calculated the multirun mean before averaging with other models. We calculated mean β values for a given model and across models for a given experiment to find regions of model agreement (Figs. 2–4). We chose to use β because this parameter allowed us to compare the relative powers of variability across the spectrum, as opposed to individually analyzing standard deviations of various frequencies, which would only allow us to compare distinct portions of the spectrum.
We replicated our analysis with near-surface air temperature (SAT; or reference height temperature) in CESM LME, and there was minimal difference from TS (“skin” temperature) β values. We also calculated temperature β between 1/2- and 1/100-yr frequencies and found slight differences in magnitude, but no difference in the spatial patterns of our results. Additionally, we repeated our analysis for the CESM LME experiments after removing the linear trend from all the time series, and our β value results did not significantly change (Fig. 5).
For the CMIP5 LM and CMIP5 PIC experiments, we excluded all models with a global mean trend (model drift) greater than 1°C per millennium (GFDL CM3, MIROC-ESM, HadGEM2-AO) after finding a relationship between preindustrial control simulations with a large global temperature drift and global mean temperature β values. Once we removed models with a temperature drift greater than 1°C per millennium, the relationship between global mean temperature trend and global mean β in CMIP5 PIC was greatly reduced. There is some degree of regional drift in many of the remaining models that slightly increases our estimates of low-frequency variability, but as in Fredriksen and Rypdal (2016), we find that removing the full time series linear trend in CMIP5 PIC models minimally impacts the distribution of β (Fig. 5).
c. Surface temperature observations
We compared CMIP5 surface temperature spectra to spectra calculated from gridded temperature data in the Niño-3.4 (5°S–5°N, 170°–120°W), North Atlantic (45°–55°N, 50°–20°W), South Atlantic (50°–40°S, 50°–20°W), and North Pacific (20°–40°N, 150°–130°W) regions. We used the full, nondetrended time series from four instrumentally based data products: HadCRUT4 extends back to 1850 CE (Morice et al. 2012), HadISST to 1870 CE (Rayner et al. 2003), Kaplan SST to 1856 CE (Kaplan et al. 1998), and NOAA ERSST.v4 to 1854 CE (Huang et al. 2015). SST observations are increasingly sparse farther back in time, and data coverage is particularly problematic in the South Atlantic and tropical Pacific before the mid-twentieth century (e.g., Deser et al. 2010). The SST data products statistically fill in missing observations with interpolated values based on particular methodologies that are described in the original references (e.g., Smith et al. 1996; Rayner et al. 2003; Huang et al. 2015). We show results using these data products in Figs. 6 and 7.
In Fig. 6, we show the model and observation spectra for the four ocean regions. We did not standardize these time series prior to spectral estimation so we could directly compare the magnitude of variance at a given frequency among the observed and simulated datasets. We calculated the annual mean anomalies for each time series, estimated the power spectrum for each time series, and averaged the individual spectra for each grid box to calculate the regional mean spectrum. We show the average spectrum (and 95% confidence bounds) from the four observation datasets for each region in Fig. 6.
d. Coral data
We compared CESM LME and CMIP5 temperature β to a network of Indo-Pacific coral records. We selected all tropical Indian and Pacific Ocean centennial-scale coral records available from the Past Global Changes 2000 yr (PAGES 2k) network (Ahmed et al. 2013) and the NOAA paleoclimate coral data archive in the spring of 2016 (http://www.ncdc.noaa.gov/paleo/coral/coral_data.html). We excluded coral records with less than 100 yr of continuous data to ensure our analysis would begin to resolve centennial-scale variability in our calculations of β (see Table 2 coral record information). We first binned the monthly coral data into annual (January–December) mean bins, and then we calculated β [(1/2) − (1/record length frequencies)] for each (nondetrended) coral record. We shaded the coral β in a circle using the same color bar as the model data at the corresponding latitude and longitude of the coral sample location (Fig. 8). If multiple coral records at least 100 yr in length were available from the same location, we averaged their β values before plotting. We chose to show the coral β on both model mean temperature and precipitation maps because these records can be sensitive to both sea surface temperatures and sea surface salinity; salinity often reflects changes in evaporation, precipitation, or runoff. We estimated values of coral-record β across multiple frequency ranges (1/2–1/500 and 1/2–1/100 yr) and found that all coral records show some degree of power-law scaling across frequency ranges. Although each coral record’s β changes slightly if a new frequency range is selected, the mean β of all coral sites effectively does not change between the two frequency ranges mentioned above. Additionally, we calculated β for all pre-twentieth-century coral records with at least 100 yr of data in this time period to test whether the low-frequency variability in coral records is solely an artifact of twentieth-century anthropogenic forcing. Although the pre-twentieth-century β value changed slightly for each coral record, we found a minimal difference in the mean β of all the coral records (Table 2).
e. Coral δ18O forward modeling
We estimated spectra from observed proxy coral δ18O and model-derived (NCAR CCSM4) “pseudocoral” δ18O at the PAGES 2k network coral sites (PAGES 2k δ18O records highlighted with bold font in Table 2). We used the Thompson et al. (2011) forward model, which applies a linear regression to combine temperature and salinity variability, to calculate forward-modeled pseudocoral δ18O. We used the temperature and salinity data from the CCSM4 last millennium and historical simulations at the PAGES 2k network δ18O coral sites. For each coral site, we used the model data from the ocean grid box closest to the coral proxy site and the model data that overlapped temporally with the coral proxy data.
f. Monte Carlo simulations of climate variability
We used a Monte Carlo procedure to test if model β values were significantly different than β values generated by a simple autoregressive process. Here we are testing the null hypothesis grid point by grid point; in the climate system, there is a spatially coherent structure to variability, and we are not applying a multivariate approach to our analysis that takes into account the spatial autocorrelation in the climate system (Newman et al. 2016). We first generated 10 000 time series with AR(1) set to the autocorrelation for each model grid box monthly time series (after removing the annual cycle). We then calculated the β values for each of these time series (as with the actual data) to provide upper and lower confidence limits for AR(1) β values (Fig. 9b). A simple autoregressive process can generate a wide range of β values, but the AR(1) spectrum flattens at decadal–multidecadal scales, and a power-law fit is only a rough approximation of AR(1) scaling (Fig. 9) (Ault et al. 2013b). An AR(1) process with monthly lag-1 autocorrelation less than 0.6 typically generates a nearly flat spectrum, with β values rarely greater than 0.2 (Fig. 9). The mean lag-1 autocorrelation of monthly temperature in CMIP5 is 0.56 (±0.27), whereas precipitation is essentially white noise with a lag-1 autocorrelation mean of 0.08 (±0.1).
g. Global impacts of the temperature scaling in the Niño-3.4 region
We compared simulated values of β for Niño-3.4 SST to decadal–multidecadal drought frequency in El Niño–Southern Oscillation (ENSO) teleconnected regions in CMIP5 LM. We calculated drought frequency using a “running mean method”: we defined a drought as a unique instance when the 11- (decadal drought) or 35-yr (multidecadal drought) running mean of the normalized precipitation time series passed the −0.5 standard deviation threshold below the full time series mean (Ault et al. 2014). We then calculated the mean drought count (per century) over land for the following regions: western North America (WNA; 26°–50°N, 125°–100°W), Australia (AUS; 40°–10°S, 112°–154°E), northern Amazonia (AMAZ; 12°S–10°N, 70°–35°W), and the northern Middle East and central Asia (MIDEAST; 30°–45°N, 35°–70°E). We picked these regions based on strong ENSO teleconnections in CMIP5 models and in observations.
We also compared tropical Pacific temperature variability to global temperature variability in CMIP5 LM. We calculated the mean temperature β from all nontropical Pacific grid boxes and compared this average to the β of the Niño-3.4 index. Additionally, we compared the mean simulated Niño-3.4 β to the value of β from the first mode of surface temperature variability in each CMIP5 LM model. We calculated the modes of variability after removing the linear trend from all grid boxes. We also tested the degree to which local temperature β is related to the magnitude of the local ENSO teleconnection in the CESM LME control simulation; we calculated the local r value of the ENSO temperature teleconnection and regressed the magnitude of this relationship with ENSO against the local temperature β value.
We find that CMIP5 temperature spectral shapes (values of β) vary regionally (consistent with Henriksson et al. 2015; Fredriksen and Rypdal 2016). We observe the lowest temperature β values in the tropics and the highest β values over mid- and high-latitude oceans (Figs. 2–4). CMIP5 temperature spectra in the tropical Pacific tend to have the most power at interannual frequencies, whereas over higher-latitude oceans the temperature spectra exhibit more power at lower frequencies (Fig. 6) (Fredriksen and Rypdal 2016). Although the tropical Pacific and Indian Oceans appear to be blue, the spectral shape in these regions is not appropriately described by a linear fit due to the large interannual variance peak (Figs. 1d and 2; consistent with Henriksson et al. 2015; Fredriksen and Rypdal 2016).
Last millennium and control simulations show similar patterns. In CESM LME, temperature values of β are typically positive across the globe, except in the tropical Pacific and Indian Oceans and in some parts of the Southern Ocean (Fig. 2). CMIP5 LM temperature values of β are also positive across the globe except in the tropical Pacific and some regions in the Southern Ocean in certain models (Fig. 3). CMIP5 PIC temperature values of β are slightly lower than CMIP5 LM values, but the spatial patterns remain consistent between the LM and PIC forcing, with the highest values at higher latitudes over the oceans and lower values near the tropics (Figs. 4 and 10).
Simulated temperature values of β exhibit both higher amplitudes (positive and negative) and greater consistency in the sign of β than precipitation β across models and experiments (Figs. 11–13). Simulated precipitation is nearly white at most locations across the globe (Fig. 1e), and the sign of precipitation β is much more variable than the sign of temperature β (Fig. 1f), even among different ensemble members of the same model with identical forcing (Fig. 11). However, the sign of CMIP5 precipitation β tends to loosely follow the sign of temperature β; the tropical Pacific tends to have low values of β for both variables (Figs. 1d,f), and many regions over high-latitude oceans are consistently slightly red for both variables (Fig. 10).
In CESM LME, there is agreement among different ensemble members that precipitation spectra over the Indian Ocean, the central and eastern tropical Pacific, and land in the tropics are dominated by interannual variability. By contrast, precipitation spectra over the far-western tropical Pacific, isolated patches of the South and North Pacific gyres, North Atlantic, and Southern Ocean regions are all slightly red in CESM LME (Fig. 11). CMIP5 LM models also show that precipitation β values are low over the tropical Pacific and slightly positive over much of the Atlantic Ocean and polar oceans (Fig. 12). Similarly, many CMIP5 PIC models agree that precipitation is dominated by interannual variability in the tropical Pacific and by low-frequency variability over the high-latitude oceans (Fig. 13). The largest disagreement among different models, in terms of the magnitude of precipitation β, is in the North Atlantic (Figs. 14d,f).
c. Testing against an AR(1) null hypothesis
Some degree of redness in the climate system may arise from ocean damping of weather noise (Hasselmann 1976; Newman et al. 2003, 2016). Here we evaluate the magnitude of CMIP5 β values against the null hypothesis that increasing variability with time scale is generated exclusively from smoothing of high-frequency forcing. We find that CMIP5 temperature β values are distinct from β values generated by a simple autocorrelated [AR(1)] process over most of the high-latitude oceans and continents, where temperature β values exceed ~0.15 (Figs. 2–4). However, over much of the central and eastern tropical Pacific and Southern Oceans, CMIP5 LM models do not generate temperature β values distinct from AR(1) β values. The low-latitude area with low β values expands in unforced CMIP5 models; most CMIP5 PIC models do not generate β values distinct from β values generated by an AR(1) process across most of the tropics. In fact, much of the midlatitude temperature variability in many CMIP5 PIC models is difficult to distinguish from variability generated by an AR(1) process (Fig. 4).
Our results indicate that most CMIP5 models generate temperature variability distinct from monthly autocorrelation over high-latitude oceans (Fig. 6) (Manabe and Stouffer 1996). These model-based results support the assertion that Hasselmann-type processes can explain much of CMIP5 models’ unforced temperature behavior over low-to-midlatitude oceans (Hasselmann 1976; Manabe and Stouffer 1996; Ault et al. 2013a). However, although an AR(1) process can generate β values that are similar in magnitude to the tropical Pacific CMIP5 temperature β values, the shape of the spectrum in the cold tongue region differs from an AR(1) spectrum in many CMIP5 models (Fig. 6a). Although there is some degree of surface ocean memory process (lag-1 monthly autocorrelation of ~0.88) in the tropical Pacific, the mechanisms that generate interannual variability include well-described coupled ocean–atmosphere interactions (e.g., Zebiak and Cane 1987), and an AR(2) model more accurately describes this variability (e.g., Neumaier and Schneider 2001).
Most CMIP5 models generate precipitation variability that falls outside the variability expected from a simple autocorrelated process in isolated midlatitude and high-latitude regions (Figs. 11–13). However, not all CMIP5 models show this behavior in the same geographic locations, even among different ensemble members of the same model (Fig. 11). We find that ~15% of all CMIP5 LM and PIC precipitation β values fall outside the 95% confidence interval (greater than the statistically expected 5% and far less than the ~75% of grid boxes for simulated TS in CMIP5 LM). Based on these results, we suggest simulated precipitation variability is effectively white in most geographic locations, but local ocean–atmosphere coupling or other climate variability can influence precipitation β. For example, in the tropical Pacific where sea surface temperature and precipitation tend to be coupled, interannual variance dominates the spectrum in both temperature and precipitation (Figs. 1d,f). However, in other regions with low-frequency temperature variability, precipitation variability remains white (Figs. 1c,e).
d. Coral comparison
CMIP5 models generally show that the variance spectra of central and eastern tropical Pacific temperature and precipitation are white or dominated by interannual variability (Figs. 1, 2–4, and 11–13). However, when we estimate β from coral paleoclimate records across the tropical Indo-Pacific Oceans, we find that none of these records reproduce this spectral pattern (Fig. 8). All of the coral records analyzed here show red spectra, and there is no interannual variance peak dominating the spectrum or disrupting the continuum of variability. In fact, all coral records we analyze suggest power-law scaling behavior in the tropical Pacific and Indian Oceans (mean coral β = 0.85 ± 0.23). This finding is consistent regardless of our use of temperature-, precipitation-, or salinity-sensitive coral records or of our use of skeletal δ18O, Sr/Ca, or luminescence records. In other words, these high values of tropical Pacific and Indian Ocean coral β are robust regardless of coral record type or location, suggesting that either CMIP5 models tend to underestimate the magnitude of β in this region or that corals are incorporating a nonclimatic red noise signal.
e. External forcing in CMIP5 LM
External forcings used in the CESM LME experiments include orbital, well-mixed greenhouse gases, land use and land cover, volcanic forcing of stratospheric aerosols, and solar variations (Schmidt et al. 2012). We compare the impact of preindustrial era forcing on β in the CESM LME and find greenhouse gas, land-use and land-cover, orbital, and solar variability have minimal impacts on LM variance distributions (Fig. 15). The volcanic-only forcing in the CESM LME is the only external forcing that appreciably increases β values in simulated temperature (Fig. 15b). This volcanically forced increase in temperature β is particularly strong in the tropics in CESM (land only) and in GISS-E2-R (land and ocean) (Figs. 10 and 16). We observe a similar signal when we compare the CMIP5 LM (all forcing, including volcanic forcing) and CMIP5 PIC simulations: a large increase in temperature β, particularly over the tropics (Figs. 8c, 8e, and 10c). Although volcanic forcing impacts the temperature spectrum in CESM and GISS-E2-R, external, preindustrial forcing has little impact on the modeled precipitation spectrum (Figs. 8d, 8f, and 10d). Precipitation β values in volcanically forced CMIP5 simulations are only slightly higher than β values in unforced CMIP5 simulations (Figs. 10, 12, and 13).
We provide an overview of temperature and precipitation scaling in multicentury simulations of the latest (CMIP5/IPCC AR5) Earth system models, under both unforced and pre-industrially forced conditions. Similar to Fredriksen and Rypdal (2016), we find that temperature scaling in these models is highly variable regionally, particularly weak over the tropical Indo-Pacific, and stronger over the high-latitude oceans. Precipitation scaling is similarly weak over much of the tropical Indo-Pacific and is dominated largely by interannual variability in the central and eastern tropical Pacific. We focus on the tropical Indo-Pacific for data–model comparison, where corals provide an unusually widespread and unsmoothed paleoclimate record that indicates consistently strong power-law (red) behavior. Preindustrial external forcing is insufficient to bring model and paleodata results into agreement for decadal and longer variability estimates. Here we address specific issues raised by these discrepancies.
a. Intermodel variance disagreements in CMIP5
Among different CMIP5 models, we see large differences in the magnitudes of temperature β. These disagreements are regionally focused, with the largest differences in the Southern, tropical Pacific, and North Atlantic Oceans (Fig. 14). Although we removed CMIP5 models with large global temperature drifts from our analysis, local climate drift over the Southern Ocean could contribute to temperature β disagreement in this region (Sen Gupta et al. 2013). Additionally, this disagreement could be due to internal variability in climate behavior; the length of the model runs (500+ yr) may not capture the full spectrum of variance (e.g., Wittenberg 2009). Disagreement among models could also result from differences in how internal variability is modeled in key parts of the global climate system [e.g., Atlantic overturning circulation (MacMartin et al. 2013), Southern Ocean circulation (Russell et al. 2006), and Pacific dynamics (Bellenger et al. 2014)].
The largest magnitude of precipitation β disagreement among different models is in the North Atlantic (Fig. 14). This disagreement could be due to the differences in the simulation of Atlantic overturning circulation among CMIP5 models (MacMartin et al. 2013) and the impact of differing local low-frequency temperature variability on the precipitation spectrum. Models that tend to have high temperature β in the North Atlantic also tend to have high precipitation β (r = 0.74).
b. Possible causes of scaling behavior over high-latitude oceans
Our results indicate that most CMIP5 models generate temperature variability distinct from monthly autocorrelation over high-latitude oceans (Fig. 6) (Manabe and Stouffer 1996). Although ocean circulation and coupled dynamics are prime candidates for generating low-frequency variability in CMIP5 models, scaling behavior can be present in high-latitude oceans even without fully coupled, complex deep ocean circulation. For example, Barsugli and Battisti (1998) find that coupling between the atmosphere and ocean can explain the increase in low-frequency variability observed in Manabe and Stouffer (1996), even in the absence of oceanic circulation. Studies have also shown that adding deeper ocean circulation and various mixed layer depths to the simple Hasselmann (1976) model can increase low-frequency variability (e.g., Fraedrich et al. 2004; Barsugli and Battisti 1998; Frankignoul and Hasselmann 1977; Lemke 1977). Therefore, although CMIP5 temperature β over the high-latitude oceans may not be consistent with a simple Hasselmann-type model, this apparent increase in low-frequency variability with time scale may be consistent with a two-layer Hasselmann-type model with atmospheric forcing (e.g., Fraedrich et al. 2004).
c. Discrepancies among variance estimates in models, proxies, and instrumental observations
The coral records analyzed here indicate power-law scaling behavior across the tropical Pacific and Indian Oceans, but the instrumental record β does not always agree with historical-era coral β (Figs. 7a,b). Coral and instrumental records tend to agree that low-frequency variability defines the temperature spectrum across much of this region. However, instrumental records show no such behavior in the cold tongue region of the eastern equatorial Pacific (Fig. 7a), and CMIP5 spectra agree quite well with the short observation record (Fig. 6a). The two coral records in or near this region (Galapagos and Palmyra) show stronger scaling than the instrumental record (Fig. 7a). Yet the observation record of tropical Pacific temperatures is extremely sparse as recently as 1955 CE, and the record loses consistent coverage before ~1975 CE (e.g., Fairbanks et al. 1997; Deser et al. 2010). Furthermore, these statistically in-filled instrumental data depend on stationary relationships between poorly and better-observed ocean regions and therefore may not accurately resolve multidecadal–centennial-scale spectral estimates. Additionally, different sea surface temperature (SST) datasets disagree on the magnitude of long-term trend in the Pacific and thus may disagree on low-frequency variability estimates (Deser et al. 2010). The observation record therefore cannot rule out low-frequency climate variability in the tropical Pacific.
Paleoclimate records show that the magnitude of ENSO variance has changed throughout the Holocene (e.g., Cobb et al. 2013; McGregor et al. 2013), and multicentury coral records indicate that decadal-scale tropical Pacific variability was stronger before the mid-twentieth century (Damassa et al. 2006; Ault et al. 2009). Similarly, long climate model simulations show intervals of stronger and weaker interannual ENSO variability (e.g., Wittenberg 2009). However, lower-frequency variability in both temperature and precipitation appears considerably different in CMIP5 models (weak) as compared to paleoclimate data (strong) (Fig. 8).
We further examine the coral–CMIP5 tropical temperature spectral discrepancies by forward modeling coral δ18O at the PAGES 2k network δ18O coral sites, using SST and salinity from the historical and last millennium CCSM4 simulations to create pseudocoral records (Thompson et al. 2011). We find a general agreement in interannual–decadal pseudocoral and observed coral δ18O variability. However, observed coral δ18O shows more variance at multidecadal–century time scales than pseudocoral δ18O (Fig. 17). Similarly, the interannual instrumental temperature variance at these sites agrees well with CCSM4 SST variance, but the model and mean instrumental spectra begin to deviate at multidecadal time scales. However, the least variable of the observation-based data products (HadISST and Kaplan SST) agrees with multidecadal temperature variability in CCSM4 (Fig. 17).
Outside the tropics, CMIP5 models and instrumental observations agree that temperature variance increases with time scale (Fig. 6). However, even if the overall shape of the instrumental and the simulated spectra are similar, many models are considerably more variable than the instrumental record across the spectrum in regions such as the North Atlantic (Fig. 6b) (Laepple and Huybers 2014a). In other regions, such as the South Atlantic, many models overestimate interannual variability and underestimate multidecadal variability (Fig. 6c) (Laepple and Huybers 2014a).
d. Possible explanations for proxy–model mismatch
Possible reasons for the proxy–model mismatch include biases in paleoclimate proxies, shortcomings in Earth system models, and uncertainties in the last millennium forcing reconstructions (Schmidt et al. 2012). Potential model shortcomings include insufficient energy cascades in models (Ferrari and Wunsch 2009), modeled ocean–atmosphere interactions that are too weak (e.g., Cane 1998), or lack of modeled slow climate feedbacks related to carbon cycling or ice sheets (e.g., Lovejoy et al. 2013; Bakker et al. 2017). In terms of coral biases, environmental stress may alter biological fractionation of stable isotopes in some corals (Damassa et al. 2006), and few observations exist to calibrate coral records over century-plus time scales. Moreover, proxies may track multiple climate variables. For example, coral isotope records reflect a mixed signal of temperature and salinity, which have distinct spectra of variability. Combining SST and sea surface salinity (SSS) in forward models shows that trends (one possible source of apparent low-frequency variability) in coral δ18O are consistent with modeled twentieth-century trends (Thompson et al. 2011).
The shape of the CMIP5-derived pseudocoral spectra closely resembles the shape of CMIP5 SST spectra, and local SSS is even more white than SST at the same geographic locations (Fig. 17). Therefore, including modeled SSS in our calculations does not resolve the coral–CMIP5 spectral mismatch. Although forward modeling of proxies using climate model data decreases much of the interannual paleo-model frequency disagreement, neither proxy forward modeling (e.g., Dee et al. 2017) nor inverse modeling (Laepple and Huybers 2014b) can remove the low-frequency (centennial–millennial scale) proxy–model mismatch. However, the proxy forward models used in these studies may exclude processes that influence low-frequency climate variability.
The mismatch between coral and modeled variance spectra may result from biases at either end of the frequency spectrum. Our coral analysis suggests that CMIP5 models may overestimate relative high- versus low-frequency variability in the tropical Pacific, consistent with previous work (e.g., Guilyardi et al. 2009; Ault et al. 2013a; Bellenger et al. 2014; Henriksson et al. 2015). However, individual proxy records of surface temperature indicate that CMIP5 models may underestimate local low-frequency temperature variability (Laepple and Huybers 2014b). Despite the presence of low-frequency internal processes in climate models, such as centennial-scale variations in Pacific SST, deep ocean circulation, and multicentury variability in interhemispheric ocean transport (Manabe and Stouffer 1996; Delworth and Zeng 2012; Karnauskas et al. 2012), paleodata show local low-frequency temperature variability that is up to 100 times stronger (Laepple and Huybers 2014b). This mismatch suggests that unless we are misinterpreting paleoclimate records, the latest generation of climate models substantially underrepresents low-frequency processes.
Additionally, missing or inaccurate external forcing estimates used in the CMIP5 LM experiment could explain part of the paleo-model mismatch. Although there is a short satellite record of solar activity, prior to the satellite era, proxy and geophysical measurements must be used to reconstruct solar activity, and the magnitude of this presatellite era solar activity is not well known (Schmidt et al. 2012). Furthermore, volcanic eruptions drive stratospheric aerosol concentrations in the last millennium experiment; these aerosols block solar radiation from reaching Earth’s surface and help drive short-lived cooling in CMIP5 models (Masson-Delmotte et al. 2013). Various methods have been used to infer the pre-instrumental concentration of aerosols in the stratosphere; these methods often disagree in the magnitude of this forcing and each has particular errors and uncertainties (Gao et al. 2008; Crowley et al. 2008).
e. Role of external forcing in temperature variance distributions
External forcing, such as solar and orbital variability, is well expressed in multimillennial paleoclimate reconstructions (Masson-Delmotte et al. 2013). However, in a comparison of individual forcing runs over the past millennium in the CESM LME, solar and orbital forcing do not appreciably change the relative strength of low- versus high-frequency variance (Figs. 15e,f). An underestimation of the magnitude of solar forcing in the last millennium experiment, or the climate system’s sensitivity to such forcing, could contribute to this paleo-model mismatch.
Various CMIP5 models also indicate that volcanic forcing increases long-period temperature variability in the climate system over both the last millennium and the instrumental era (Zhong et al. 2011; Laepple and Huybers 2014b; Henriksson et al. 2015; Vyushin et al. 2004; Vyushin and Kushner 2009). Indeed, the volcanic-only forcing in the CESM LME is the only external forcing that seems to appreciably increase simulated temperature β values (Fig. 15).
Although external (volcanic) forcing appears to increase low-frequency temperature variance in the tropics in CMIP5, this forcing does not increase model temperature β above AR(1) β values in the Pacific cold tongue region (Figs. 6a, 2, and 3). We suggest that simulated interannual ENSO variance dominates the spectrum in “strong ENSO” models (such as BCC_CSM1.1, CCSM4, CESM, and FGOALS) and that external natural forcing is insufficient to overpower the simulated high-frequency tropical internal variability (Ault et al. 2013a). By contrast, ENSO variance plays a smaller role in the overall spectrum in “weak ENSO” models (such as GISS-E2-R). In these models, integration of short-lived volcanic forcing leads to a larger impact on temperature β in the tropics (Figs. 6, 10, and 16). Outside the tropics, volcanic eruptions have a weaker impact on temperature values of β in most CMIP5 models (Fig. 6), likely because 1) the temperature impact of eruptions is small compared to the large internal variability of temperature at high latitudes, and 2) low-frequency internal variability already dominates the spectrum at high latitudes in CMIP5 PIC (Fig. 6), so volcanic forcing does not change the relationship of low- versus high-frequency variance.
The nonvolcanic, preindustrial forcings do not appreciably impact the β values of temperature (Fig. 15). If the forcing estimates and the modeled climate sensitivity to the forcing are correct, then this result suggests that volcanic forcing is the primary external forcing of simulated temperature of the preindustrial era in the last millennium (850–1850 CE). By contrast, paleoclimate records suggest that solar (and internal) variability influence decadal–centennial-scale temperature variability (Masson-Delmotte et al. 2013) and that temperature impacts of volcanic forcing are limited to interannual–decadal scales over the past 2000 yr (Masson-Delmotte et al. 2013; Sigl et al. 2015).
f. Role of external forcing in precipitation variance distributions
Although volcanic forcing impacts the temperature spectrum in CESM and GISS-E2-R, external, preindustrial forcing has less impact on the modeled precipitation spectrum (Figs. 10d and 8d,f). Precipitation β values in volcanically forced CMIP5 simulations are slightly higher than β values in unforced simulations in the tropics and at higher latitudes in the Northern Hemisphere (Fig. 10). However, this result should be interpreted cautiously. Previous work suggests that in models, ENSO and the monsoon respond to volcanic forcing in the opposite direction of that indicated by paleoclimate records (Ault et al. 2013a; Anchukaitis et al. 2010). This ambiguity contrasts with the impacts of volcanic eruptions on temperature; paleoclimate data agree with models that volcanic aerosols cause at least short-lived cooling.
g. Global impact of tropical Pacific variance
Paleoclimate records from the tropical Pacific suggest that low-frequency temperature variability should dominate across the tropics, whereas CMIP5 models and the relatively short observation record suggest the temperature spectrum is either flat or dominated by interannual variance in the cold tongue region (Fig. 8). Tropical Pacific temperature variability has a near-global climate imprint (e.g., Trenberth et al. 1998); power-law behavior in the tropical Pacific therefore has the potential to impact variability in many regions. To examine this possibility, we compare models that simulate weak and strong low-frequency variability in the tropical Pacific. The models that simulate a tropical Pacific dominated by high-frequency variability tend to experience less decadal–multidecadal drought in ENSO-teleconnected regions (Figs. 18a–d), whereas prolonged drought is more common in the models that simulate stronger low-frequency Pacific variability. We also compare tropical Pacific temperature variability to global temperature variability in CMIP5 LM and find a strong correlation between β values in the Niño-3.4 region and global temperatures (Figs. 18e and 18f, respectively). Recent work related to the slowdown of recent atmospheric temperature increase supports our findings that Pacific climate variability influences global atmospheric temperature variability (e.g., Trenberth and Fasullo 2013; Dai et al. 2015; Thompson et al. 2015).
Although we observe a strong relationship among tropical Pacific β and global temperature (and precipitation) β, correlation does not imply causation; tropical Pacific temperature β may not directly impact global temperature and precipitation variance distributions. These different parts of the climate system may all experience low-frequency variability due to a combination of external forcing and internal dynamics that increases low-frequency variability globally and does not necessarily originate in the tropical Pacific. Nonetheless, maps of CESM temperature β are strikingly similar to ENSO-teleconnection maps in this model. In fact, when we compare the absolute magnitude of the local ENSO temperature teleconnection to local β from the same CESM simulation, we find a strong relationship (r = −0.60). This finding suggests that interannual tropical Pacific variability can play an important role in shaping temperature variance distributions in ENSO-teleconnected regions, and ENSO biases in models may impact regions far from the tropical Pacific.
5. Summary and conclusions
Accurate estimates of the spectrum of climate variability are essential to understanding the background dynamics of climate itself. Low-frequency climate variability likely arises from a combination of physical modes of variability and smoothing of weather noise by oceans (Manabe and Stouffer 1996). One outcome of our work is that we can identify where simulated climate variability could arise solely from local autocorrelation and where other factors also help determine the shape of the spectrum.
Our analysis of CMIP5 precipitation and temperature variability suggests that simulated precipitation is effectively atmospheric noise with little internal memory or low-frequency variability, regardless of forcing, except over isolated patches of mid- and high-latitude ocean regions (Figs. 12 and 13). By contrast, simulated temperature is slightly red on interannual–century time scales over most of the globe outside the tropics (Figs. 3 and 4) (Fredriksen and Rypdal 2016). Unforced CMIP5 models produce internally generated, significant low-frequency temperature variability over most of the high-latitude oceans (Manabe and Stouffer 1996; Fredriksen and Rypdal 2016). However, unforced CMIP5 temperature scaling is difficult to distinguish from a simple AR(1) process across most of the tropics and midlatitudes; outside the central and eastern tropical Pacific much of this variability could arise from an integration of weather noise by thermal inertia in the oceans. In forced CMIP5 simulations of the last millennium, low-frequency temperature variability, especially in the tropics, appears to arise from the persistence of short-lived cooling from volcanic eruptions.
Many CMIP5 models suggest that the variance spectrum of tropical temperature and precipitation is dominated by interannual variability. By contrast, many paleoclimate records, including multicentennial coral records that should be robust at decade–century time scales, suggest that low-frequency variability dominates the last millennium across the tropical oceans. Although the tropical Pacific temperature spectrum in CMIP5 matches well with the spectrum from the instrumental record over interannual–decadal scales, the instrumental record is too short to accurately resolve multidecadal–centennial-scale frequencies, so we cannot independently test whether paleoclimate records or climate models are more accurately representing the low-frequency behavior of the climate system (Fig. 6a). The proxy–model mismatch could be attributed either to an unknown aspect of the proxy system or to model variance that is too strong on interannual scales and/or too weak on multidecadal–millennial scales (Laepple and Huybers 2014b), perhaps due to model and/or forcing issues.
Clearly, CMIP5 models and paleoclimate data paint drastically different pictures of the natural variability of tropical climate. Especially in the near term, the natural patterns of decade–century-scale variability will determine future climate and related risks. Understanding and anticipating these risks requires a deeper knowledge of the internal dynamics that lead to the apparent red continuum of climate variability.
The National Science Foundation EaSM2 Grant (AGS-1243125), the Directorate for Geosciences (Grant 3008610), and the Graduate Research Fellowship under Grant DGE-1143953 supported this work. We also thank the Kartchner Caverns scholarship fund and the Department of Geosciences at the University of Arizona for their funding support. J. Emile-Geay provided a script for normalizing precipitation time series. B. Otto-Bliesner at NCAR provided the CESM last millennium ensemble data. B. Otto-Bliesner, J. Fasullo, and S. Stevenson at NCAR also provided valuable feedback. We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table 1) for producing and making available their model output. For CMIP the U.S. Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-16-0863.s1.