Research article

Earth Syst. Dynam., 12, 1393–1411, 2021
https://doi.org/10.5194/esd-12-1393-2021© Author(s) 2021. This work is distributed under
the Creative Commons Attribution 4.0 License.

Ubiquity of human-induced changes in climate variability
Keith B. Rodgers1,2 , Sun-Seon Lee1,2 , Nan Rosenbloom3 , Axel Timmermann1,2 , Gokhan Danabasoglu3 ,Clara Deser3 , Jim Edwards3 , Ji-Eun Kim1,2 , Isla R. Simpson3 , Karl Stein1,2 , Malte F. Stuecker4 ,Ryohei Yamaguchi1,2 , Tamás Bódai1,2 , Eui-Seok Chung5 , Lei Huang1,2 , Who M. Kim3 ,Jean-François Lamarque3 , Danica L. Lombardozzi3 , William R. Wieder3,6 , and Stephen G. Yeager3
1 Center for Climate Physics, Institute for Basic Science, Busan, South Korea
2 Pusan National University, Busan, South Korea
3 National Center for Atmospheric Research, Boulder, CO, USA
4 Department of Oceanography and International Pacific Research Center, School of Ocean and Earth Scienceand Technology, University of Hawai‘i at Mānoa, Honolulu, HI, USA
5 Korea Polar Research Institute, Incheon, South Korea
6 Institute of Arctic and Alpine Research, University of Colorado, Boulder, CO, USA
Correspondence: Keith B. Rodgers (krodgers@pusan.ac.kr) and Axel Timmermann (axel@ibsclimate.org)
Received: 25 June 2021 – Discussion started: 6 July 2021
Revised: 5 October 2021 – Accepted: 9 October 2021 – Published: 9 December 2021

Abstract. While climate change mitigation targets necessarily concern maximum mean state changes, under-

standing impacts and developing adaptation strategies will be largely contingent on how climate variability
responds to increasing anthropogenic perturbations. Thus far Earth system modeling efforts have primarily focused on projected mean state changes and the sensitivity of specific modes of climate variability, such as the
El Niño–Southern Oscillation. However, our knowledge of forced changes in the overall spectrum of climate
variability and higher-order statistics is relatively limited. Here we present a new 100-member large ensembleof climate change projections conducted with the Community Earth System Model version 2 over 1850–2100 to
examine the sensitivity of internal climate fluctuations to greenhouse warming. Our unprecedented simulations
reveal that changes in variability, considered broadly in terms of probability distribution, amplitude, frequency,
phasing, and patterns, are ubiquitous and span a wide range of physical and ecosystem variables across many
spatial and temporal scales. Greenhouse warming in the model alters variance spectra of Earth system variables
that are characterized by non-Gaussian probability distributions, such as rainfall, primary production, or fire occurrence. Our modeling results have important implications for climate adaptation efforts, resource management,seasonal predictions, and assessing potential stressors for terrestrial and marine ecosystems.

1

Introduction

Faced with the prospect of substantial future climatechange, mitigation and adaptation strategies are increasingly
paramount. While mitigation efforts are concerned chiefly
with limiting mean state changes, successful adaptation will
also require understanding the potentially altered variabilityof the climate system (Sarachik, 2010). However, the way inwhich climate variability will change due to anthropogenic
radiative forcing has not been extensively explored. The
spectrum of observed regional-to-global climate fluctuations

exhibits relatively sharp peaks and a broad noise background(Hasselmann, 1976; Franzke et al., 2020). Spectral peaks
can emerge from a range of mechanisms, including astronomical forcings or internal climate instabilities such as for
the El Niño–Southern Oscillation (ENSO). Moreover, these
distinct features can be further influenced by climate processes acting on different timescales. Examples of non-linear“timescale interactions” are multiplicative (state-dependent)noise (Müller, 1987; Majda et al., 2009; Sardeshmukh andSura, 2009; Sardeshmukh and Penland, 2015; Jin et al., 2007,
2020; Levine and Jin, 2010) and combination mode dynam-

Published by Copernicus Publications on behalf of the European Geosciences Union.

1394

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

ics (Stuecker et al., 2015b). How modes of variability will respond to greenhouse warming has been addressed in a number of previous modeling studies (Timmermann et al., 1999;
Cai et al., 2018), albeit with conflicting results. In contrast,
the sensitivity of the spectral background to human-induced
climate change is less well-known. Identifying and characterizing human-induced changes in this spectral background,using for example Coupled Model Intercomparison Project(CMIP)-type coordinated modeling efforts, has proven difficult due to limited statistics.
The relatively recent advent of large ensemble simulations (henceforth termed large ensembles) conducted with
Earth system models provides a new resource for addressing
how climate and ecosystem statistics may evolve in responseto anthropogenic forcing across a wide range of timescales(Deser et al., 2020; Schlunegger et al., 2020). Such large ensembles with global climate models have existed for morethan 15 years (Zelle et al., 2005; Drijfhout et al., 2008;
Branstator and Selten, 2009), but earlier studies expressed
concern with aspects of process representation and therefore
their results with regard to variability changes were inconclusive. Other studies have employed individual model simulations, small (≤ 10 members) ensembles, or CMIP multimodel ensembles (Rind et al., 1989; Raisanen, 2002; Huntingford et al., 2013; Screen, 2014; Stouffer and Wetherald,2007; Wetherald, 2009) to address whether surface temperature and precipitation variability may change under global
warming. To date large ensemble studies of changes in variance have mainly focused on specific quantities, timescales,or regions (Deser et al., 2020; Pendergrass et al., 2017; Maher et al., 2019, 2021; Haszpra et al., 2020). However, to ourknowledge, the full power of the large ensemble frameworkhas not been harvested to gauge broad-scale forced changes
in climate statistics, including changes in variance, spectrum,
patterns, and phase, for a wide range of quantities, regions,or timescales.
To study the sensitivity of higher-order climate statistics to anthropogenic climate change, we conducted a new
100-member ensemble of climate change simulations using
the Community Earth System Model version 2 (CESM2)(Danabasoglu et al., 2020), which we refer to as CESM2LE (Methods). The initialization and forcing are describedin the Methods section and in Figs. S1–S5 in the Supplement. An ensemble of this size and duration with a CMIP6generation Earth system model at 1◦ spatial resolution is unprecedented. A large number of improvements have occurredsince CESM1-LE (Kay et al., 2015), as documented in the
Methods section. In addition to improved parameterizations
and process representation that increase model skill in representing a number of phenomena, a notable improvement
is also seen in land processes within the Community Land
Model Version 5 (CLM5).
CESM2-LE promises to provide an enhanced frameworkfor documenting and understanding robust forced changes ininternal variability, complementing our knowledge of meanEarth Syst. Dynam., 12, 1393–1411, 2021

state changes (Simpson et al., 2020; Fasullo, 2020). The
simulations were performed for the 1850–2100 period withhistorical (1850–2014) and SSP3-7.0 (2015–2100) forcings.
The choice to use 100 members was motivated by the challenges associated with identifying trends in higher statistical
moments. A previous set of analyses performed with the Max
Planck Institute Grand Ensemble (MPI-GE) (Milinski et al.,2020) explored the relationship between ensemble size andthe accuracy of identifying forced changes in higher-ordermoments. Even taking into account differences in model architecture, and thereby model uncertainty in such estimates,
their analysis with the MPI-GE nevertheless supports our decision to expand well beyond the 40 members chosen for
CESM1-LE (Kay et al., 2015). To facilitate analysis over abroad range of timescales, substantial resources have been
devoted to providing high-frequency output for the atmosphere, land, ocean, and cryosphere. Providing a clearer viewof the patterns of altered climate variability should enable investigation of the mechanistic drivers of such changes andtheir implications for impacts of societal and ecosystem relevance. This study presents initial results on forced changes
in internal variability across a range of fields and timescalesin CESM2-LE, and it will serve as the reference publication
for CESM2-LE.
2
2.1

Methods
Model configuration

The simulations consist of a 100-member large ensemble
suite conducted with CESM2 with the Community Atmosphere Model version 6 (CAM6) (Danabasoglu et al., 2020),
referred to here as CESM2-LE. The simulations cover the period from 1850–2100 and follow the historical and SSP3-7.0
forcing protocols provided by CMIP6 (Eyring et al., 2016),
although with some differences noted below for the representation of biomass burning in 50 of the 100 ensemble members. CESM2 has been demonstrated to fare well when evaluated against skill metrics with other models (Fasullo, 2020).
The choice of the SSP3-7.0 scenario forcing follows CMIP6
recommendations (O’Neill et al., 2016) that emphasize thevalue of this relatively high forcing level precisely for the
purpose of quantifying forced changes in natural variability.
This choice should also provide a useful contribution towards
an eventual CMIP6 large ensemble intercomparison.
The CESM2 components use nominal 1◦ horizontal resolution. Specifically, CAM6 has a resolution of 1.25◦ in longitude and 0.9◦ in latitude, and 32 vertical levels with a top
at 2.26 hPa, or approximately 40 km. The ocean and sea ice
models are the Parallel Ocean Program version 2 (POP2)(Danabasoglu et al., 2020; Smith et al., 2010) and the CICE
Version 5.1.2 (CICE5) (Bailey et al., 2020). The nominal resolution of the ocean is 1◦ horizontally, with uniform spacingof 1.125◦ in the zonal direction and varying significantly inthe meridional direction, with the finest resolution of ∼ 0.25◦
https://doi.org/10.5194/esd-12-1393-2021

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

at the Equator. The ocean model provides 60 vertical levels,with 20 of these layers represented in the upper 200 m of thewater column. CESM2 offers a number of improvements pertinent to our scientific interests relative to what was available
for CESM1-LE (Kay et al., 2015). These improvements include advances in the surface boundary layer representation
for the ocean (Li et al., 2016), as well as for cloud microphysics (Gettelman et al., 2015). The ocean biogeochemistrymodel used with the POP2 model is the Marine Biogeochemistry Library (MARBL), which represents an updated versionof what was previously known as the biogeochemistry elemental cycle (BEC) (Moore et al., 2001, 2004, 2013; Long etal., 2021).
An important advance of great value to large ensemble
investigations is achieved through new developments incorporated into CLM5 (Danabasoglu et al., 2020; Lawrence etal., 2019; Lombardozzi et al., 2020). The model addresses anumber of well-known limitations relative to previous versions of CLM, including major improvements in simulated
cumulative CO2 uptake over the historical period (Bonan etal., 2019) and improved representation of the seasonal cycleof net ecosystem production (NEP) (Lawrence et al., 2019),
which is highlighted in our analysis of projected forced
phenology changes. Other notable features also included in
CLM5 are the explicit representation of agricultural management and improvements in the implementation of the prognostic fire model (Lombardozzi et al., 2020; Li et al., 2013;
Li and Lawrence, 2017). All CLM5 improvements foundbroadly across a range of simulated variables have been documented through evaluation of model simulations against
the International Land Model Benchmarking (ILAMBv2.1)
package and other analyses (Collier et al., 2018; Danabasoglu et al., 2020). We note that land model trajectories are
sensitive to SSP scenarios that determine the spatial distribution and extent of land use and land cover changes (O’Neillet al., 2016).
As a more general complement to the research results considered in this study, we have also made available resultsfrom running the Climate Variability Diagnostics Package
for large ensembles (CVDP-LE) (https://www.cesm.ucar.
edu/working_groups/CVC/cvdp-le/, last access: 19 November 2021) (Phillips et al., 2020) for CESM2-LE, withgraphical output available under https://climatedata.ibs.re.
kr/data/cesm2-lens/lens-diagnostics (last access: 19 November 2021).
2.2

Large ensemble initialization

For the CESM2-LE initialization procedure, the experimental configuration was designed to respond to broad community demand for a mix of macro- and micro-perturbations(where for micro-perturbations members differ only in asmall random perturbation applied at initialization). To satisfy this demand and allow for exploration of the impact
of initialization type, it was decided to initialize members
https://doi.org/10.5194/esd-12-1393-2021

1395

from various years between 1001 and 1301 of a preindustrial simulation conducted with CESM2 (Danabasoglu et al.,2020). This was as far as the CESM2(CAM6) preindustrial
simulation had reached at the time when the CESM2-LE
project began, and by this point the top-of-the-atmosphere(TOA) global energy imbalance was relatively small and stable with a correspondingly small model drift (Danabasoglu etal., 2020). The years from the preindustrial control run used
for initialization are highlighted in Figs. S1 and S3.
Micro-initializations start from four different years: 1231,
1251, 1281, and 1301. A total of 20 members were run foreach start year, with ensemble spread introduced by a random perturbation to the atmospheric temperature field at initialization (through a CAM6 namelist variable referred to as“pertlim”), as was used for all members of CESM1-LE (Kayet al., 2015). Macro-initialization (one run for each initialization date) used initialization years {1001, 1011, 1021, . . .,1191} by using 20 independent restart files at 10-year intervals over 1001–1191. It warrants mention for the case of
the macro-perturbations that no explicit perturbation was required from the preindustrial control simulation. Taken together, if one includes one member from each of the microperturbation runs, then a total of 24 macro-perturbation runs
are available.
Importantly, as can be seen in Fig. S1b, for the initialization points of years 1231, 1251, 1281, and 1301 were specifically chosen for the micro-initializations to correspond toyears of maximum, decreasing, minimum, and increasing
Atlantic Meridional Overturning Circulation (AMOC) transport, respectively, relative to the preindustrial control simulation. It is important to note that when using the large
ensemble output, the initialization procedure should not be
considered to produce members that are independent, or tohave randomized modes of climate variability, for the years
immediately subsequent to 1850. Considering the AMOC
strength at 26.5◦ N as an example (Fig. S2), the ensemble
mean AMOC strength for each of the micro-perturbation
clusters initialized for years 1231, 1251, 1281, and 1301 of
the preindustrial control run (averaged across 20 members
for each case) converge only after several decades, indicativeof the timescale over which the initial condition memory persists for AMOC. For this reason, our analysis with internalvariability focuses on the period after 1960, more than a full
century after initialization. Further quantitative explorationof the specific duration over which initial condition memory
is retained is the subject of a separate ongoing study.
A generalized schematic for the initialization procedure isshown in Fig. S3, illustrating the organization of the simulations. The schematic also includes mention of the biomassburning emissions differences between two groups of 50
simulations, as described more fully in the next section. The
macro-perturbation runs initialized at {1011, 1031, 1051,. . ., 1191} have greatly enhanced output at high frequency
to meet the needs of broader community interests for large
ensemble output. The temporally high-resolution output inEarth Syst. Dynam., 12, 1393–1411, 2021

1396

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

cludes 6-hourly snapshots of three-dimensional temperature,
winds, and specific humidity for the Coordinated Regional
Climate Downscaling Experiment (https://cordex.org, lastaccess: 19 November 2021) simulations, as well as output
appropriate for the Cloud Feedback Model Intercomparison Project (CFMIP) Observation Simulation Package(COSP)(https://climatedataguide.ucar.edu/climate-data/cosp-cloud-feedback-model-intercomparison-project-cfmipobservation-simulator-package, last access: 19 November 2021).
2.3

Large ensemble forcing

A choice was made to use two different sets of forcing fields
to represent the effects of variability in biomass burningemissions for CESM2-LE (see Figs. S4–S5). The biomassburning aerosol fluxes in CESM2 are imposed at the surface. As such, they are not prognostic, meaning that they
are not generated by the model’s internal prognostic firemodel. The first 50 members of our large ensemble follow
CMIP6 protocols (Van Marle et al., 2017), with biomassburning following the description in the CESM2 overview
paper (Danabasoglu et al., 2020), and this forcing is referred to as BB_CMIP6. For the second set of 50 members, which we refer to as BB_CMIP6_SM (for smoothedbiomass burning fluxes), the BB_CMIP6 biomass burningemissions of all relevant species for CAM6 were smoothedin time with an 11-year running mean filter. The averaging
impacted variability in biomass burning fluxes over 1990–
2020. Due to the inclusion of observations, the variability
in biomass burning emissions during 1990–2020 is considerably stronger for BB_CMIP6 than the preceding and following periods. The smoothed forcing with BB_CMIP6_SM wasdesigned to nearly conserve total emissions, while reducing
the strong changes in interannual variability. The temporalsmoothing of the forcing is applied to the biomass burningemissions at each grid point subsequent to being regridded to
the CESM2 grid. The high 1990–2020 biomass burning variability case (ensemble members 1–50, or BB_CMIP6) relative to the smoothed forcing (ensemble members 51–100,
or BB_CMIP6_SM) has a discernible impact on large-scale
climate, as documented by the accelerated loss of September Arctic sea ice and northern hemispheric and tropical Pacific warming (Fig. S5a and c). Outside of the period 1990–
2010, the impact of BB_CMIP6_SM relative to BB_CMIP6
for biomass burning emissions is not pronounced for simulated surface temperature, sea ice, or precipitation. It is forthis reason that we selected the time intervals 1960–1989 and
2070–2099 for our analysis of variance changes in Figs. 2
and 4, for which the 100 ensemble members can realisticallybe considered to be part of the same population.

Earth Syst. Dynam., 12, 1393–1411, 2021

2.4

Minor corrections relative to previous versions

The code base for the BB_CMIP6_SM simulations (the second set of 50 members) incorporates corrections for two
sets of errors that were present in the first set of 50 ensemble members (BB_CMIP6). The first pertains to the SO2 ,SO4 , and gas-phase semi-volatile secondary organic aerosol(SOAG) emission datasets. For SO2 and SO4 , the spatial
patterns of the “shipping” and “agriculture+solvents+waste”
components of forcing were inadvertently switched during
the historical-to-projection transition, or more specifically at
the start of 2015. The incorrect partitioning of SO2 does
not impact the results considered here, given that its components are summed before use. In contrast, the issue with
SO4 datasets can impact the model state evolution as the particle sizes and numbers differ for the SO4 components. The
SOAG emissions are calculated from several hydrocarbons,and they were not recalculated after an earlier bug correction in covering units of the lumped species for the biomassburning emissions. This issue was corrected, and diagnostics indicate that there are not any pronounced changes in themodel solutions from these particular aerosol corrections.
The second correction introduced for the 50
BB_CMIP6_SM simulations concerns the presence ofa sporadic large CO2 uptake over land that was identified
for the BB_CMIP6 runs. This large uptake is associatedwith a negative flux of carbon occurring at crop harvest
time over a single time step. Although these large negative
carbon flux component terms in autotrophic respiration arenecessary for maintaining carbon balance, such CO2 spikes
are not realistic. To avoid these spikes, the associated CO2
fluxes that occur over a single time step are distributed tothe atmosphere over a timescale of approximately 6 months
for the BB_CMIP6_SM simulations. Analysis indicatesthat these modifications for carbon between the BB_CMIP6
and BB_CMIP6_SM simulations did not result in any
climate-changing impacts.
3
3.1

Results
Mean state changes

During the historical period the evolution of key simulated annual-mean climate indicators in CESM2-LE (Figs. 1
and S6) agrees well with observations. The range acrossthe ensemble members, which results from internal variability and its forced changes, spans the observed climate statemuch of the time, with a notable exception being Southern Ocean sea ice (Fig. 1e). The results here and the general model behavior are qualitatively consistent with those ofsimilarly forced CMIP6-generation models (Fasullo, 2020;
Kwiatkowski et al., 2020; Arora et al., 2020), although projected temperature changes (Fig. 1c) are in the upper range of
the CMIP6 models owing to the relatively high climate sensitivity of CESM2 (Gettelman et al., 2019). The progressive
https://doi.org/10.5194/esd-12-1393-2021

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

1397

Figure 1. Annual mean evolution of global fields over 1850–2100 for 100 ensemble members. For model fields, bold lines represent ensemblemeans, and dark and light shading represent 1 standard deviation (SD) and 2 SD variability. Observational data are shown in red whenappropriate. Portions of the figure with light-blue background shading indicate the historical period (1850–2014) while portions with lightred background shading indicate the projection period (2015–2100). (a) Top-of-atmosphere radiative imbalance (W m−2 ) along with the
CERES-EBAF product (Loeb et al., 2009, 2018); (b) anomalies of the global mean precipitation (mm d−1 ) increasing 5.4 % between the
1850s and the 2090s, compared with the Global Precipitation Climatology Project (GPCP) (Adler et al., 2003, 2012); (c) anomalies ofglobal mean surface temperature, increasing by 4.4 ◦ C between the 1850s and 2090s, along with HadCRUT4 (Morice et al., 2012) anomalies
over 1950–2019; (d) anomalies of ocean heat content integrated over the upper 2000 m, along with an observation-based product (Ishii etal., 2017); (e) anomalies of sea ice extent for the Arctic (black) and Southern Ocean (blue), with observed sea ice extent over 1979–2020(Fetterer et al., 2017), and with the vertical scales of the anomaly plots offset to facilitate comparison; (f) Atlantic Meridional Overturning
Circulation (AMOC) transport anomalies at 26.5◦ N, with RAPID array observations (Frajka-Williams et al., 2019); (g) globally integrated
net primary productivity (NPP) over the ocean (blue; increase of 2.7 % between the 1850s and 2090s), and over land (green); and (h) globally
integrated net CO2 fluxes over the ocean (solid blue) and integrated net CO2 flux (net biome production, or NBP, including fire and land-use
change) over land (green) with all quantities in (g) and (h) in units of Pg C yr−1 . For each case, where observational products are included,
anomalies are calculated with respect to the period spanned by the observations. For anomaly fields, printed numbers represent the absolutemean of the ensemble mean of CESM2-LE (black or blue numbers) and the observational product (red numbers).

https://doi.org/10.5194/esd-12-1393-2021

Earth Syst. Dynam., 12, 1393–1411, 2021

1398

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

weakening of the AMOC at 26.5◦ N in CESM2 over the 21st
century (Fig. 1f) is largely consistent with other CMIP6 models (Weijer et al., 2020). We also find a substantial increasein land primary productivity (Fig. 1g), which contributes tothe uptake of carbon in the terrestrial biosphere. Marine net
primary productivity (NPP) (Fig. 1h) remains relatively constant throughout the simulation, and the overall uptake of carbon by the ocean reflects the re-emergence of anthropogenic
carbon into the mixed layer (Toyama et al., 2017; Rodgers
et al., 2020) and changes in the CO2 buffering capacity ofseawater (Revelle and Suess, 1957). For the analysis that ispresented in Fig. 1 for sea ice, daily-mean output fields areused for both the model and the data product. In representing
sea ice extent a threshold of 15 % was used, whereby a gridcell is identified as being ice covered if it has a concentrationof sea ice above 15 %. For the net land fluxes of CO2 , we use
the variable net biome production which includes the effects
of not only photosynthesis and respiration, but also fire andland-use change.
The pattern of mean state surface temperature change,
shown as the difference between the periods 2070–2099 and
1960–1989 (Fig. 2, central panel; 2 m reference temperatureshown in Fig. S6), exhibits preferential warming of the eastern relative to the western equatorial Pacific, Arctic amplification, and a pronounced warming hole over the subpolar North Atlantic. These features are associated with the
known mechanisms of the enhanced equatorial warming pattern (Xie et al., 2010), more positive polar feedbacks (Goosse
et al., 2018) including the Arctic heat capacitor (Chung etal., 2021), and the slowdown of the AMOC (Rahmstorf etal., 2015; Menary and Wood, 2018), respectively. For precipitation (Fig. 2, central panel; Fig. S6e), changes include
marked precipitation increases along the equatorial Pacific,
within the Arctic Ocean and decreases over the subtropicalregions (Stocker et al., 2013).
3.2

Forced changes in amplitude, frequency, and phase

Figure 2 illustrates the ensemble- and grid-point-aggregatedFourier amplitude spectra and probability density functions(PDFs) for five key climate and ecosystem quantities (complementary quantities are shown in Fig. S7). The choice ofvariables reflects an interest in both climate and ecosystem
dynamics, as well as societal relevance in terms of adaptation
and resource management. The decision to represent Fourier
amplitude spectra was motivated by our desire to enrich ourunderstanding of the amplitude of perturbations across different timescales. For the spectral analysis in Fig. 2, each fast
Fourier transform (FFT) spectrum is calculated for the timeseries of raw data over a given variable for the full 30-year
interval. This includes all timescales shorter than 30 years
and longer than 2 days (months) for daily (monthly) timeresolution data. The spectrum is calculated first at each horizontal grid point and for each ensemble member and thenaveraged over the designated region and over the 100 enEarth Syst. Dynam., 12, 1393–1411, 2021

semble members. Due to the relatively large degree of aggregation for each field, it was not necessary to apply windowing to avoid spectral leakage. The surface chlorophyll
concentration fields analyzed here represent total chlorophyll
concentrations taken as a sum of diatom, diazatroph, andsmall phytoplankton chlorophyll concentrations. The AMOC
in Fig. S7 is defined as a maximum transport at 26.5◦ N.
For the spectrum of internal variability of the AMOC, theensemble-mean is subtracted from the raw data to remove
the forced response
For a wide range of Earth system variables, we find substantial changes of the projected 21st century probability distributions, impacting mean state, variance, and higher-order
statistical moments (Fig. 2). Human-induced alterations ofclimate spectrum and probability distribution could translate into changes in the average return time of climate andextreme events. Averaging the spectra over 100 ensemble
members and individual grid boxes within each region creates enough data to reveal spectral characteristics that mightotherwise be obscured. The dominant feature for most quantities examined is the spectral peak at the annual frequency,
along with higher-order harmonics that result from deviations of the seasonal cycle from a pure sinusoid. Futurechanges of the annual cycle overtone spectrum can be causedby forced non-sinusoidal distortions of the annual cycle, generated for example by shifts in phenology, as discussed below. For nearly all variables under consideration, the seasonal cycle amplitude responds to the external forcing. Nearannual combination modes (C modes) of ENSO and the seasonal cycle (Stuecker et al., 2015a) and its overtones can beclearly identified in some spectra, particularly for precipitation over the equatorial Pacific. In addition to representing the C modes as deterministic components of the system, CESM2-LE also exhibits shifts in the frequency of the
C modes due to future reductions in ENSO’s dominant frequency (Fig. 3a). The C-mode peaks also strengthen in the
future, reflecting that the amplitude of precipitation and thecorresponding C-mode-generating non-linearity increase atboth ENSO and annual frequencies.
For most of the variables shown in Fig. 2 (and Fig. S7)
there are changes in the amplitude of the spectrum across theentire range of frequencies from synoptic to intra-seasonal
to interannual to decadal, revealing the ubiquity of variancechanges. Importantly, frequency-independent shifts in variance can be seen in the three variables shown here, which exhibit a strong non-Gaussian skewed PDF, namely the spectraof California wildfire occurrence, surface chlorophyll concentrations over the subpolar North Atlantic (40◦ –60◦ N,
60◦ –15◦ W), and precipitation over the Niño 3.4 region(5◦ S–5◦ N, 170◦ –120◦ W). For these positive variables with
their highly skewed probability distributions, forced changes
in the mean state are accompanied by a stretching (squeezing) of the associated PDFs, thereby causing enhancement(or reduction) of variance and extremes. Changes of this typehave previously been considered for more specialized cases
https://doi.org/10.5194/esd-12-1393-2021

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

1399

Figure 2. Changes in the Fourier amplitude spectrum of historical (1960–1989) to future (2070–2099) climate variability in CESM2-LE.

The center map shows historical-to-future changes in surface temperature (shaded, ◦ C) and precipitation (solid blue/cyan dots, mm d−1 ).
Each pairing of panels shows historical (cyan) and future (red) spectra and PDFs for five different variables over four different regions.
The spectra are considered over the respective periods, 1960–1989 (historical) and 2070–2099 (future), thereby including the trend, and
PDFs are considered for all days over 1980–1989 and 2090–2099 to minimize the impact of the trend. From upper-left clockwise, eachpair of panels shows fire occurrences in California (32◦ –41◦ N, 125◦ –118◦ W, land only), surface chlorophyll concentrations in the North
Atlantic subpolar gyre (40◦ –60◦ N, 60◦ –15◦ W), net ecosystem production (NEP) in the Amazon (10◦ S–10◦ N, 80◦ –50◦ W, land only),
precipitation over the Niño 3.4 regions (5◦ S–5◦ N, 170◦ –120◦ W), and sea surface temperature (SST) over the Niño 3.4 region. The spectra
are calculated for daily time series at individual grid points including both forced responses and internal variability and using 30-year
intervals. Subsequently the spectra are averaged over the grid points in each region. Sharp spectral peaks are associated with the annual
cycle and its non-sinusoidal components, which generate high-order harmonics. Shaded areas for spectra of precipitation and temperature in
the Niño 3.4 region correspond to the timescales of the El Niño–Southern Oscillation (ENSO) and ENSO-annual-cycle combination modes(Stuecker et al., 2013) (C modes). Spectra are shown as amplitude, with the units being the same as the x axes for the PDFs. PDFs of positive
variables (California fire counts, North Atlantic surface chlorophyll, and Niño 3.4 precipitation) are shown with logarithmic y axes. The
fields in the center panel are presented in more detail in Fig. S6, except that there 2 m reference temperature is used rather than surface
temperature. A suite of complementary spectral and PDF analyses to those shown here are presented in Fig. S7.

using the Wasserstein distance (Ghil, 2015; Robin et al.,2017; Vissio et al., 2020). For white noise processes, the associated variance changes manifest as timescale-independent
variance changes, thereby accounting for the shown spectral
background shifts. For California fire counts and Niño 3.4
precipitation, mean state increases are therefore also accompanied by increases in variance occurring over a wide rangeof timescales. For North Atlantic chlorophyll concentrations, the mean state decrease is associated with a timescale-https://doi.org/10.5194/esd-12-1393-2021

independent decrease in variance, with expected impacts forhigher trophic levels in the ocean, leading to potential disruptions to ecosystems.
For variables that are less skewed, a diversity of responses
is found. Forced changes in sea surface temperature (SST)
variability in the Niño 3.4 region are confined to interannual
timescales in association with a decrease in ENSO amplitude and a slight shift toward higher frequencies. On the other
hand, for NEP over the Amazon, reflecting natural CO2 ex-

Earth Syst. Dynam., 12, 1393–1411, 2021

1400

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

Figure 3. Changes in the dominant frequencies and seasonal variance of sea surface temperature (SST, left) and precipitation (right), in the

Niño 3.4 region (5◦ S–5◦ N, 170◦ –120◦ W). The top row shows the wavelet power spectra of Niño 3.4 (a) SST and (b) precipitation using a
Morlet wavelet, normalized by σ −2 , where σ is the ensemble mean standard deviation of the respective Niño 3.4 time series (Torrence andCompo, 1998). The y axis shows the equivalent Fourier period in years. The hatching indicates regions where the wavelet spectrum is nottrustworthy due to edge effects. Prior to calculating the wavelet spectra, the time series were detrended by subtracting the ensemble-meanannual means, which were linearly interpolated to a monthly time step. The middle row shows the ensemble-mean of Niño 3.4 (c) SST
and (d) precipitation indicated for each day (ordinate) and year (abscissa) using daily output. The red (blue) dots indicate the maximum(minimum) daily values of each year. The black line to the right in panels (c–f) indicates the linear trend over 1960–2100. The bottom row
shows the same as for (c) and (d), but for the across-ensemble SDs of (e) SST and (f) precipitation.

change between the land and the atmosphere, there is an increase in variance over all timescales, accompanied by a shift
in the broad interannual peak towards higher frequencies.
To demonstrate the significance of the forced change forthe spectra considered in Figs. 2 and S7, an example is given

Earth Syst. Dynam., 12, 1393–1411, 2021

in Fig. S8 for precipitation over the Niño 3.4 region forthe same 1960–1989 (blue) and 2070–2099 (red) time intervals. The grey shading indicates the 95 % confidence interval(1.96 × standard error). For each ensemble member, we first
spatially averaged the spectra at individual grid points over

https://doi.org/10.5194/esd-12-1393-2021

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

the Niño 3.4 region and then calculated the ensemble-wise(or across-ensemble) error using the 100 spectra for the full
100 ensemble members. This approach appropriately handles
spatially correlated data in the calculation of the confidence
intervals and is known as the block bootstrap method, where
the block means are the spatial means. The estimated confidence intervals in Fig. S8 indicate that the spectra for 1960–
1989 and 2070–2099 are statistically different. It is worth
noting here that if samples at different grid points are treatedas being independent samples, the confidence intervals become much narrower, in which case the two spectra in Fig. S8
are even more statistically significant in their difference.
We next turn our attention to an expanded view of the temporal evolution of both frequency and amplitude modulationsof SST and precipitation over the Niño 3.4 region over theperiod 1960–2100. For the wavelet analysis (Torrence andCompo, 1998) in Fig. 3, we apply a Morlet wavelet normalized by σ −2 , where σ is the ensemble mean standard deviation of the respective time series. Ensemble wavelet analysis of SST (Fig. 3a) and precipitation (Fig. 3b) within the
Niño 3.4 region has been conducted after first removing theensemble-mean trend over the full period from each ensemble member, while retaining the seasonal cycle. The wavelet
analysis is conducted for each ensemble member and then averaged. Our motivation for retaining the seasonal cycle stems
from an interest in illustrating timescale interactions between
ENSO and the seasonal cycle with the full power of large ensemble statistics. The annual cycle and ENSO interact witheach other in a complex way, with the annual cycle itself
being a forced coupled air–sea mode (Xie, 1994). This interaction gives rise to combination modes (Stuecker et al.,2015b), frequency entrainment (Timmermann et al., 2007),and ENSO’s phase-locking and seasonal variance modulations (Stein et al., 2010, 2014). Not only does the annual
cycle in the equatorial Pacific influence the amplitude andphase of ENSO, but ENSO also impacts the seasonal cycle.
We consider the normalized variance to highlight the amplification above the white noise level, and in contrast to
Fig. 2 represent variance with a linear scale to emphasize
temporal modulation of the amplitude of the maxima. For
SST a clear separation is seen between the maxima for interannual variability and the annual cycle (Fig. 3a). At interannual timescales, there are two notable features. The firstis a shift in the ENSO peak period from 3.5 to 2.5 years between the end of the 20th century and the end of the 21st century. The second feature with interannual variability is thatvariance does not change monotonically but rather exhibits amaximum midway through the 21st century, similar to what
has been reported elsewhere (Kim et al., 2014). This stands
in contrast to precipitation over the same region (Fig. 3b),for which there is a monotonic increase in variance, following a similar shift in the period of the peak that was found for
SST. For precipitation, the amplitude of the seasonal cycle
increases over 1960–2100, consistent with the notion of vari-https://doi.org/10.5194/esd-12-1393-2021

1401

ability enhancement over the tropics due to thermodynamic
and dynamic processes (Yun et al., 2021).
The forced changes over 1960–2100 in the structure of theseasonal cycle for the ensemble mean of SST (Fig. 3c) and
precipitation (Fig. 3d), as well as the across-ensemble standard deviation of SST (Fig. 3e) and of precipitation (Fig. 3f)
are also considered for the Niño 3.4 region using daily-mean
model output. The maximum (red dots) of ensemble-mean
SST occurs in May and the minimum (blue dots) in October in the late 20th century (Fig. 3c), with both showing
monotonic increases over 1960–2100. The maximum shifts
to 2 weeks later and the minimum shifts to 2 weeks earlier
by the end of the 21st century, with this modest perturbation to the phase of the seasonal cycle being accompanied
by a modulation of seasonal amplitude. The ensemble-meanseasonal amplitude in precipitation (Fig. 3d) occurs approximately 1 month before the ensemble-mean maximum in SST(Fig. 3c), and a second maximum in precipitation in late
January becomes evident during the second half of the 21st
century. On the other hand, the ensemble-mean minimum in
precipitation occurs approximately 2 weeks after the local
minimum in temperature. The increase in the amplitude of
the seasonal cycle is thereby accompanied by changes in thephasing of the seasonal cycle for both SST and precipitation.
The mechanisms responsible for the phasing of maximum precipitation leading maximum temperature over the
Niño 3.4 region over seasonal timescales (red dots in Fig. 3cand d) have been considered previously in published literature (Xie, 1996; Xie et al., 2010; Williams and Patricola,
2018; Stuecker et al., 2020). Current understanding maintains that seasonal precipitation phasing is largely driven by
meridional SST gradients and is thereby not directly tied tothe phasing of seasonal SST variations in the Niño 3.4 region.
In other words, the phase relationship between precipitationand SST is not surprising, as moisture convergence is in part
determined by non-local SST conditions.
The seasonally stratified maximum across-ensemble SD in
SST (Fig. 3e), associated with peak ENSO variability, exhibits a trend towards an earlier occurrence by approximately
1 month over 1960–2070. This is accompanied by a modest
decrease in amplitude (line plot). The across-ensemble SD
minimum for SST occurs in July for the 20th century, witha secondary minimum in the across-ensemble SD developing over the first half of the 21st century in May. Subsequently the across-ensemble SD minimum in May becomes
more pronounced and becomes the dominant minimum inthe across-ensemble SD of SST by the end of the 21st century. For the across-ensemble SD of precipitation (Fig. 3f),
there is a monotonic strengthening of the seasonal maximum
in late January, corresponding roughly to the time of peak
ENSO variability, and a weakening of the seasonal minimum in October, over the interval 1960–2100. Whereas theseasonal minimum in the across-ensemble SD of precipitation (Fig. 3f) occurs nearly in phase with the seasonal minimum of ensemble-mean SST (Fig. 3c), the seasonal maxEarth Syst. Dynam., 12, 1393–1411, 2021

1402

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

imum for the across-ensemble SD of precipitation does notcoincide with the seasonal maximum of ensemble-mean SST.
Rather, it coincides with the secondary seasonal maximum inensemble-mean precipitation in late January (Fig. 3d).
3.3

Changes in variance and co-variance patterns

Along with modulations in the frequency domain, the spatial patterns of variance are altered in response to changing
climate conditions. The analysis of patterns of variance andco-variance in Fig. 4 uses across-ensemble calculations of
seasonal mean SDs (top) and Pearson correlation coefficients(bottom), respectively. These calculations entail first calculating the SD or correlation coefficients across all ensemble
members for the same time record. Subsequently averaging
is done across time. This sequence was chosen to avoid spurious amplification of variability due to the non-trivial forced
variations in precipitation and surface temperature driven by
volcanic aerosols over the historical period. To examine peak
ENSO variance, surface temperature (Fig. 4a and c) and precipitation (Fig. 4b and d) were averaged over December–
January–February (DJF), and then the across-ensemble SDs
and correlation coefficients were calculated separately over
all years spanning 1960–1989 and 2070–2099, and then averaged over the two respective periods. The intention with thecalculation of both across-ensemble SDs and correlations is
to harness the full power of the large ensemble and is analogous to the empirical orthogonal function (EOF) EOF-E
snapshot method (Maher et al., 2018).
To calculate the statistical significance of the changes, thegeneral approach is to first calculate the equivalent sample
size n̂, to account for potential serial correlations of the timeseries. This is then used to calculate the degrees of freedom
for the Welch’s t test, which is an adjusted version of theStudent’s t test that allows for the two samples to have unequal variance (i.e., heteroskedasticity). First, the decorrelation timescale Te was calculated at each grid point, and foreach period, based on the e-folding timescale of the autocorrelation function r(τ ), defined as the smallest lag τ for whichr(τ ) < e−1 . Then the equivalent sample size n̂ was defined asn̂ = N/Te , where N = 30 is the total sample size in our case.
The equivalent sample size was then used to calculate thedegrees of freedom of the standard Welch’s t test. Note that
this test may still be liberal if the equivalent sample sizes are
small, i.e., in areas of high serial correlation.
We begin by considering interannual variance changes inboreal winter (DJF) by evaluating relative changes in theacross-ensemble SD of surface temperature and precipitation
for the same periods as with the spectra in Fig. 2 (1960–1989
and 2070–2099). The background across-ensemble SD averaged over 1960–1989 is shown in shading (Fig. 4a and b),
and their respective future changes are shown in circles.
Surface temperature (Fig. 4a) reveals modest decreases invariability across the equatorial Pacific and Indian oceans,consistent with Fig. 2. Variability decreases over much ofEarth Syst. Dynam., 12, 1393–1411, 2021

the higher latitudes of the Northern Hemisphere (Screen,
2014; Screen et al., 2015; Holmes et al., 2016; Sun et al.,2015; Schneider et al., 2015), with exceptions over the Arctic
and the North Atlantic, and with exceptions in the Southern
Hemisphere found over southern Africa and parts of Antarctica (Fig. 4a). For precipitation (Fig. 4b) a relative increasein SD is seen over most regions with particularly pronounced
enhancements occurring in the eastern equatorial Pacific, theIndo-Pacific warm pool including the South Pacific Convergence Zone, the western Arabian Sea, the poles, and mostland areas. The equatorial Pacific changes represent an eastward broadening in the centers of convection in response to
the enhanced equatorial Pacific warming and the reductionof the overall zonal SST gradient (Fig. 2, center). In contrast,
there is a decrease in the northern equatorial Atlantic Ocean
as well as in some trade wind regions of the eastern Pacific.
Another important question to address is whether greenhouse warming can also impact the co-variability of different
climate components and the global teleconnections of major
modes of climate variability. This is illustrated here by examining the projected changes in the local correlation coefficients between the Niño 3.4 SST index and surface temperature from 1960–1989 and 2070–2099 (Fig. 4c), with thebackground correlation coefficients shown in shading andtheir respective future changes shown in circles. Our analysis reveals a systematic strengthening of ENSO’s remotetemperature correlation over the Amazon basin and in theequatorial Atlantic, the Philippines and Japan in the western
Pacific, throughout Africa, in northern India and across eastern Canada and the southern US. Co-variance decreases overwestern Canada and Alaska, and zonally across the equatorial Indian Ocean.
The future changes in the correlation between the Niño 3.4
index and precipitation (Fig. 4d) indicate a pattern of enhanced co-variance over the western Pacific region surrounding the Philippines, much of Africa and South America, andwestern China, as documented by the background correlation
coefficients and their future changes having the same sign. In
other words, in these regions we see stronger ENSO teleconnections under future global warming, which in turn could
translate to increased predictability of climate in the regions
on seasonal to interannual timescales, but also stronger impacts. In contrast, decreased precipitation co-variance with
ENSO is found for North America over the Pacific Northwest as well as much of the southern US and Mexico,as well as over Columbia and Venezuela, Bangladesh andMyanmar, parts of eastern Australia, and parts of eastern
Siberia. Taken together, the global pattern of ENSO precipitation co-variance changes (Fig. 4d) is due to a combinationof simulated weakening of ENSO SST variability (Fig. 4a)and eastward expansion of the region of maximum convective activity in the equatorial Pacific (Fig. 4b) (analysis for
the June–July–August (JJA) season in shown in Fig. S9),and likely other projected changes of the background atmospheric circulation. There are a number of outstanding chalhttps://doi.org/10.5194/esd-12-1393-2021

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

1403

Figure 4. Changes in the patterns of interannual variability and Niño 3.4 correlation coefficients of December–January–February (DJF)

surface temperature and precipitation. In the top row the color shading shows the time-averaged across-ensemble SD of the DJF seasonalmean surface temperature (a) and precipitation (b) for the period 1960–1989. Circles show the relative changes in the SD between 2070–2099
and 1960–1989, where insignificant change (p ≥ 0.05) has been removed. Statistical significance of the changes (circles) was determinedbased on the p values of the two-sample Welch’s t tests for the equality of temporal means of the SDs, with the equivalent sample sizes
adjusted to account for serial correlations (Sect. 3.3). In the bottom row the color shading shows ensemble-wise correlations of the Niño 3.4
index with surface temperature (c) and precipitation (d) anomalies for DJF, averaged over the period 1960–1989. Circles show the absolutechange in correlations between 2070–2099 and 1960–1989, where statistically insignificant changes (p ≥ 0.05) have been removed. The
Niño 3.4 index for ENSO is the spatial average of sea surface temperature within 5◦ S–5◦ N, 170◦ –120◦ W. Statistical significance of thechanges (circles) was determined based on the p values of two-sample Student’s t test of the Fisher z-transformed correlation coefficients(Timmermann et al., 2014). Note that the t test treats the ensemble standard deviations and correlations as stationary and serially uncorrelated
within either of the two periods. For all four panels, the circles represent subsampled fields at 10◦ intervals over the global domain. The
corresponding analysis for June–July–August (JJA) is presented in Fig. S9.

lenges in interpreting mechanistically how ENSO teleconnections change in response to anthropogenic forcing, including the relative role of local diabatic forcing and modulations of ENSO (Taschetto et al., 2020). We anticipate thatthe large ensemble analyses here will complement efforts directed at understanding mechanistic controls.

https://doi.org/10.5194/esd-12-1393-2021

3.4

Forced changes in phenology of net ecosystemproduction

Finally in this overview, we illustrate how anthropogenic
forcing impacts the phase of the seasonal cycle by focusing
on the phenology of NEP in the Northern Hemisphere middle to high latitudes (over 50◦ –80◦ N). NEP as a flux quantity
represents the difference between gross primary productionand ecosystem respiration, and thereby the net exchange ofcarbon with the atmosphere when fire and human land usechanges are ignored. Our interest in NEP is motivated by ecoEarth Syst. Dynam., 12, 1393–1411, 2021

1404

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

logical concerns that a shift to an earlier spring bloom, in particular over the land regions adjacent to the Arctic, can drivea phenological mismatch in ecological interactions between
plants and animals (Renner and Zohner, 2018). For the seasonality and phenology analysis in the upper panel of Fig. 5,
an area integral of daily-mean NEP, is performed for each ensemble member separately. A total of 90 ensemble members
are used, as daily-mean CLM5 output was not saved for the
first 10 members, namely for members {1001.001, 1021.002,
1041.003, . . ., 1181.010}.
Ensemble-mean NEP is integrated over the region in 5year intervals, with aggregation performed for individualyears and with a binning interval of 1 day (colors in Fig. 5,
upper panel). We find an evolving amplitude of the seasonal
cycle and of the growing season length (the interval during
which NEP is positive, indicating net land uptake of carbon).
This representation of forced changes in the non-sinusoidal
seasonal cycle reveals that the growing season length is projected to increase by almost 4 weeks, with the onset shifting 3 weeks earlier and termination shifting 1 week later.
The forced changes in growing season length are mostly attributable to changes in the mean temperature (Lawrence etal., 2019; Lombardozzi et al., 2020). The analysis also reveals a more than doubling of the amplitude of the seasonal
cycle in NEP as a forced response. This represents an increase in the “breathing” of the terrestrial high-latitude biosphere. Information from individual ensemble members in
20-year intervals regarding the timing of (i) first zero crossing, (ii) maximum NEP, (iii) second zero crossing, and (iv)maximum negative NEP (Fig. 5, lower panel) reveals that interannual variability (identified using one SD) is in general
smaller than the forced trend evident in the ensemble meanin spring. Our analysis indicates that for the aggregated NEP
signal, the phenological shift as a decadal trend already becomes emergent relative to the natural variability within thefirst decades of the 21st century. The trend itself is broadly
consistent with observations (Zhu et al., 2016; Myers-Smithet al., 2020). Internal variability in the date of onset of thegrowing season decreases by 35 % over the course of thesimulations and the date of the end of the growing season
decreases by 18 % (Fig. 5, lower panel). In deriving thesepercentages, the transitions (zero crossings) were first calculated individually for each ensemble member for each timeinterval (across 90 members).

4

Summary and discussion

This study introduced a new, publicly available large ensemble of climate change simulations conducted with the global
fully coupled CESM2 model. This large ensemble (CESM2LE) is unprecedented in terms of its combination of size (100members), duration (1850–2100), and spatial resolution in
the atmosphere and ocean (nominally 1◦ horizontally). As
such, it offers a unique opportunity to study not only forced
Earth Syst. Dynam., 12, 1393–1411, 2021

changes in the mean state, but also forced changes in internalvariability, including higher-order statistical moments. Here
we showcase aspects of the remarkable diversity of forced
responses in amplitude, frequency, patterns, co-variance, andseasonal characteristics of internal variability in CESM2-LE
across a broad suite of key physical and ecosystem quantities,
spanning the atmosphere, land, cryosphere, and ocean. Importantly, and contrary to conventional wisdom, the changes
are not solely centered on the frequency of specific climate
modes such as ENSO and the Madden–Julian oscillation
but are instead broadly distributed over nearly all timescales(Fig. 2), in particular for non-Gaussian distributed variables.
The mechanistic underpinnings of the changes in variability go beyond amplification or damping of major climate
modes, and possibly include state dependence of linear stabilities, non-linearities, rectification, and changes in damping
timescales and noise characteristics, many of which will be
investigated in forthcoming studies analyzing the breadth of
the CESM2-LE output fields.
If the ubiquitous changes in variance across temporal andspatial scales described here are realized in the real world,
they will have several important implications for informing
adaptation strategies and assessing potential impacts. This
holds for water resource management and agriculture, fisheries, and occurrence of wildfires. Forced changes in phenology and phasing of the seasonal cycle for ecosystem productivity pose risks of mismatches with trophic-level interactions and energy transfers. The ubiquity of such changes
in variability also points to the importance of moving beyond the assumption of stationary variability in detection
and attribution studies of climate change (Hegerl et al.,2007) and underscores the necessity of recalibrating climateeconomy models (Diaz and Moore, 2017) to account foran entirely different probability distribution for variability(Figs. 2 and S7) than what is currently used when projecting
future climate change scenarios. The non-stationary natureof climate noise under anthropogenic forcing (Fig. 2) and the
evolving teleconnections patterns (Fig. 4) also have implications for seasonal to multi-year climate predictability.
Although our analysis of CESM2-LE has revealed a broad
range of forced changes in variance across physical scales
and Earth system variables, it nevertheless should be emphasized that model uncertainty has not been consideredhere. There is already evidence for the narrower case of interannual variability in surface temperature and precipitation that model uncertainty in forced changes exhibits pronounced differences between models (Maher et al., 2021)(their Figs. 7 and 8 in the Supplement). Thus, it is our hopethat our work will motivate further investigations of forcedchange in Earth system variance across a broad range oftimescales under existing archives of large ensemble simulations (Deser et al., 2020; Schlunegger et al., 2020).
Taken together, we have provided support with new examples and new global emphasis that the Earth system is sensitive in its statistical characteristics to anthropogenic forcing,https://doi.org/10.5194/esd-12-1393-2021

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

1405

Figure 5. Expansion of growing season length, or equivalently the carbon uptake period, over 50◦ –80◦ N (shown here for all 90 members

for which daily-mean land output was saved). (a) Evolution of ensemble-mean seasonal cycle (one line for every five years, color-coded)of integrated net ecosystem productivity (NEP), with positive values indicating net terrestrial carbon uptake and negative values indicating
loss of carbon from the aggregated land region. The first zero crossing arks the start of the growing seasons, and the second zero crossing
marks the end of the growing seasons. (b) Histograms of first occurrence of zero crossing, peak, second zero crossing, and minimum as afunction of the day of the year. The horizonal axis for both panels is a climatological calendar day of the year, and aggregation is done across90 members. The histograms represent model output sampled at 20-year intervals. The inlay map (upper right) shows the ensemble mean
amplitude of the seasonal cycle of NEP averaged over 1960–1989 (g C m−2 d−1 ).

thereby building upon and complementing previous studies
that have focused on mechanistic analyses for specific phenomena (Swain et al., 2018; Tamarin-Brodsky et al., 2020;
Taschetto et al., 2020; Burger et al., 2020). Although only asmall fraction of such forced changes could be documented
in this study, we expect that the diagnostic ensemble analysis
tools applied here, along with the open access to our datasets,
will inspire further investigations into the non-stationarityof Earth system processes in the presence of anthropogenic
forcing.

The Python wavelet software used for Fig. 3 was provided by Evgeniya Predybaylo (Torrence and Compo, 1998)
and is available at http://atoc.colorado.edu/research/wavelets/(FIG_CODE_WAVELET, 2021).

Data availability. The CESM2-LE model output is available

through https://www.cesm.ucar.edu/projects/community-projects/
LENS2/data-sets.html (CESM2_LE_OUTPUT, 2021).

Supplement. The supplement related to this article is available
Code availability. The code used to generate the figures inonline at: https://doi.org/10.5194/esd-12-1393-2021-supplement.

this study is available here: https://github.com/kj-stein/CESM2-LE(FIG_CODE_GENERAL, 2021).

https://doi.org/10.5194/esd-12-1393-2021

Earth Syst. Dynam., 12, 1393–1411, 2021

1406

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

Author contributions. The CESM2-LE project was initiated

by KBR, AT, GD, and CD. The scientific framing of this paper wasdeveloped by KBR, AT, JEK, RY, KS, SSL, and MFS. Analyses andscientific post-processing were performed by RY, JEK, KS, LH, TB,
and WK. The CESM2-LE model runs were set up, performed, andextracted through a joint effort by the team of SSL, NR, and JE.
The initialization procedure for the model was developed throughthe joint efforts of CD, GD, IS, WK, SGY, and NR. All authors
discussed the results and contributed to the writing of the paper.

Competing interests. The contact author has declared that nei-

ther they nor their co-authors have any competing interests.

Disclaimer. Publisher’s note: Copernicus Publications remains
neutral with regard to jurisdictional claims in published maps andinstitutional affiliations.

Acknowledgements. The authors would like to thankAlexis Tantet and an anonymous reviewer for their thoughtful and constructive comments.
The CESM2 large ensemble (CESM2-LE) simulations presented
here for the first time have been conducted through a partnership
between the IBS Center for Climate Physics (ICCP) in South Korea and the Community Earth System Model (CESM) group at the
National Center for Atmospheric Research (NCAR) in the US, representing a broad collaborative effort between scientists from both
centers.
The authors would like to thank Woncheol Roh at the ICCP andJohn Fasullo, Keith Lindsay, Adam S. Phillips, and Gary Strand at
NCAR for their input and support. We would also like to thank all of
the other scientists, software engineers, and administrators at both
NCAR and the ICCP that contributed to this project. The framework for the macro- and micro-perturbation initialization strategy
employed here also benefited from the US CLIVAR Workshop on
large ensembles held in July 2019 in Boulder CO, USA, for which
we wish to acknowledge the support of Mike Patterson, Jennie Zhu,
and Jeff Becker at US CLIVAR.
The CESM project is supported primarily by the US National
Science Foundation (NSF). This material is based upon work supported by the NCAR, which is a major facility sponsored by the
US NSF under cooperative agreement 1852977. The CESM2 preindustrial control run was performed on the Cheyenne supercomputer (https://doi.org/10.5065/D6RX99HX) operated by the Computational and Information Systems Laboratory (CISL) at NCAR.
The simulations presented here with CESM2-LE were conducted on the IBS/ICCP supercomputer “Aleph”, a 1.43 petaflop
high-performance Cray XC50-LC Skylake computing system with
18 720 processor cores, 9.59 PB of disk storage, and 43 PB of tapearchive storage. The CESM2-LE project duration was 15 months,
and it generated 5.3 PB of data and used approximately 80 million
CPU hours of computing time.

Financial support. The work of Keith B. Rodgers, Sun-Seon Lee,
Axel Timmermann, Ryohei Yamaguchi, Ji-Eun Kim, Karl Stein,and Lei Huang was supported by the Institute for Basic Sci-

Earth Syst. Dynam., 12, 1393–1411, 2021

ences (IBS), Republic of Korea, under IBS-R028-D1. Tamas Bodai
was supported by the Institute for Basic Sciences (IBS), Republic of Korea, under IBS-R028-Y1. Malte F. Stuecker was supportedby the NOAA Climate Program Office Modeling Analysis, Predictions, and Projections (MAPP) program, grant NA20OAR4310445,
and participates in the MAPP Marine Ecosystem Task Force.
This is IPRC publication 1540 and SOEST contribution 11412.
William R. Wieder and Danica L. Lombardozzi were supportedby the National Institute of Food and Agriculture, US Department of Agriculture (2015-67003-23485). William R. Wieder was
also supported by NASA Interdisciplinary Science Program award
NNX17AK19G. The work of Nan Rosenbloom was supported by
the Regional and Global Model Analysis (RGMA) component of
the Earth and Environmental System Modeling Program of the
US Department of Energy’s Office of Biological and Environmental Research (BER) via the US NSF IA 1947282. Eui-Seok Chung
was supported by the project PE21010 of the Korea Polar ResearchInstitute. NCAR is a major facility sponsored by the US NSF undercooperative agreement 1852977.

Review statement. This paper was edited by Valerio Lucarini andreviewed by Alexis Tantet and one anonymous referee.

References
Adler, R. F., Huffman, G. J., Chang, A., Ferraro, R., Xie,P. P., Janowiak, J., Rudolf, B., Schneider, U., Curtis, S.,
Bolvin, D., Gruber, A., Susskind, J., Arkin, P., and Nelkin,
E.: The version-2 global precipitation climatology project(GPCP) monthly precipitation analysis (1979–present), J.
Hydrometeorol., 4, 1147–1167, https://doi.org/10.1175/15257541(2003)004<1147:tvgpcp>2.0.co;2, 2003.
Adler, R. F., Gu, G. J., and Huffman, G. J.: Estimating Climatological Bias Errors for the Global Precipitation Climatology Project (GPCP), J. Appl. Meteorol. Clim., 51, 84–99,
https://doi.org/10.1175/jamc-d-11-052.1, 2012.
Arora, V. K., Katavouta, A., Williams, R. G., Jones, C. D., Brovkin,V., Friedlingstein, P., Schwinger, J., Bopp, L., Boucher, O., Cadule, P., Chamberlain, M. A., Christian, J. R., Delire, C., Fisher,
R. A., Hajima, T., Ilyina, T., Joetzjer, E., Kawamiya, M., Koven,
C. D., Krasting, J. P., Law, R. M., Lawrence, D. M., Lenton,
A., Lindsay, K., Pongratz, J., Raddatz, T., Séférian, R., Tachiiri,
K., Tjiputra, J. F., Wiltshire, A., Wu, T., and Ziehn, T.: Carbon–
concentration and carbon–climate feedbacks in CMIP6 models
and their comparison to CMIP5 models, Biogeosciences, 17,4173–4222, https://doi.org/10.5194/bg-17-4173-2020, 2020.
Bailey, D. A., Holland, M. M., DuVivier, A. K., Hunke, E. C., andTurner, A. K.: Impact of a New Sea Ice Thermodynamic Formulation in the CESM2 Sea Ice Component, J. Adv. Model. Earth
Sy., 12, 15, https://doi.org/10.1029/2020ms002154, 2020.
Bonan, G. B., Lombardozzi, D. L., Wieder, W. R., Oleson, K. W.,
Lawrence, D. M., Hoffman, F. M., and Collier, N.: Model Structure and Climate Data Uncertainty in Historical Simulations of
the Terrestrial Carbon Cycle (1850–2014), Global Biogeochem.
Cy., 33, 1310–1326, https://doi.org/10.1029/2019gb006175,
2019.

https://doi.org/10.5194/esd-12-1393-2021

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability
Branstator, G. and Selten, F.: “Modes of Variability”
and Climate Change, J. Climate, 22, 2639–2658,
https://doi.org/10.1175/2008jcli2517.1, 2009.
Burger, F. A., John, J. G., and Frölicher, T. L.: Increasein ocean acidity variability and extremes under increasing atmospheric CO2 , Biogeosciences, 17, 4633–4662,
https://doi.org/10.5194/bg-17-4633-2020, 2020.
Cai, W. J., Wang, G. J., Dewitte, B., Wu, L. X., Santoso, A., Takahashi, K., Yang, Y., Carreric, A., and McPhaden, M. J.: Increased
variability of eastern Pacific El Niño under greenhouse warming, Nature, 564, 201–206, https://doi.org/10.1038/s41586-0180776-9, 2018.
CESM2_LE_OUTPUT: Model output for CESM2-LE presentedin Rodgers et al., 2021, available at: https://www.cesm.ucar.
edu/projects/community-projects/LENS2/data-sets.html, UCAR
[data set], last access: 19 November 2021.
Chung, E.-S., Ha, K.-J., Timmermann, A., Stuecker, M. F., Bodai,T., and Lee, S.-K.: Coldseason Arctic amplification driven by
Arctic Ocean-mediated seasonal energy transfer, Earth’s Future,9, e2020EF001898, https://doi.org/10.1029/2020EF001898,
2021.
Collier, N., Hoffman, F. M., Lawrence, D. M., Keppel-Aleks, G.,
Koven, C. D., Riley, W. J., Mu, M. Q., and Randerson, J. T.:
The International Land Model Benchmarking (ILAMB) System:
Design, Theory, and Implementation, J. Adv. Model. Earth Sy.,10, 2731–2754, https://doi.org/10.1029/2018ms001354, 2018.
Danabasoglu, G., Lamarque, J. F., Bacmeister, J., Bailey, D. A., DuVivier, A. K., Edwards, J., Emmons, L. K., Fasullo, J., Garcia, R.,
Gettelman, A., Hannay, C., Holland, M. M., Large, W. G., Lauritzen, P. H., Lawrence, D. M., Lenaerts, J. T. M., Lindsay, K.,
Lipscomb, W. H., Mills, M. J., Neale, R., Oleson, K. W., OttoBliesner, B., Phillips, A. S., Sacks, W., Tilmes, S., van Kampenhout, L., Vertenstein, M., Bertini, A., Dennis, J., Deser, C., Fischer, C., Fox-Kemper, B., Kay, J. E., Kinnison, D., Kushner, P.
J., Larson, V. E., Long, M. C., Mickelson, S., Moore, J. K., Nienhouse, E., Polvani, L., Rasch, P. J., and Strand, W. G.: The Community Earth System Model Version 2 (CESM2), J. Adv. Model.
Earth Sy., 12, 35, https://doi.org/10.1029/2019ms001916, 2020.
Deser, C., Lehner, F., Rodgers, K. B., Ault, T., Delworth, T. L.,
DiNezio, P. N., Fiore, A., Frankignoul, C., Fyfe, J. C., Horton, D. E., Kay, J. E., Knutti, R., Lovenduski, N. S., Marotzke,J., McKinnon, K. A., Minobe, S., Randerson, J., Screen, J. A.,
Simpson, I. R., and Ting, M.: Insights from Earth system model
initial-condition large ensembles and future prospects, Nat. Clim.
Change, 10, 277–286, https://doi.org/10.1038/s41558-020-07312, 2020.
Diaz, D. and Moore, F.: Quantifying the economic risksof climate change, Nat. Clim. Change, 7, 774–782,
https://doi.org/10.1038/nclimate3411, 2017.
Drijfhout, S., Hazeleger, W., Selten, F., and Haarsma, R.: Future changes in internal variability of the Atlantic Meridional Overturning Circulation, Clim. Dynam., 30, 407–419,
https://doi.org/10.1007/s00382-007-0297-y, 2008.
Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B.,
Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled
Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958,
https://doi.org/10.5194/gmd-9-1937-2016, 2016.

https://doi.org/10.5194/esd-12-1393-2021

1407

Fasullo, J. T.: Evaluating simulated climate patterns from the CMIP
archives using satellite and reanalysis datasets using the Climate
Model Assessment Tool (CMATv1), Geosci. Model Dev., 13,3627–3642, https://doi.org/10.5194/gmd-13-3627-2020, 2020.
Fetterer, F., Knowles, K., Meier, W., Savoie, M., and Windnagel,
A.: Sea Ice Index, version 3, National Snow and Ice Data Center,
2017.
FIG_CODE_GENERAL: Figure code for Rodgers et al., 2021,
available at: https://github.com/kjstein/, CESM2-LE [code], lastaccess: 19 November 2021.
FIG_CODE_WAVELET: Figure code Fig. 3 wavelet analysis inRodgers et al., 2021, available at: https://atoc.colorado.edu/
research/wavelets/, Department of Atmospheric and Oceanic
Sciences at the University of Colorado Boulder [code], last access: 19 November 2021.
Frajka-Williams, E., Ansorge, I. J., Baehr, J., Bryden, H. L.,
Chidichimo, M. P., Cunningham, S. A., Danabasoglu, G., Dong,
S. F., Donohue, K. A., Elipot, S., Heimbach, P., Holliday, N.
P., Hummels, R., Jackson, L. C., Karstensen, J., Lankhorst, M.,
Le Bras, I. A., Lozier, M. S., McDonagh, E. L., Meinen, C. S.,
Mercier, H., Moat, B. I., Perez, R. C., Piecuch, C. G., Rhein,
M., Srokosz, M. A., Trenberth, K. E., Bacon, S., Forget, G.,
Goni, G., Kieke, D., Koelling, J., Lamont, T., McCarthy, G. D.,
Mertens, C., Send, U., Smeed, D. A., Speich, S., van den Berg,
M., Volkov, D., and Wilson, C.: Atlantic Meridional Overturning Circulation: Observed Transport and Variability, Front. Mar.
Sci., 6, https://doi.org/10.3389/fmars.2019.00260, 2019.
Franzke, C. L. E., Barbosa, S., Blender, R., Fredriksen, H. B., Laepple, T., Lambert, F., Nilsen, T., Rypdal, K., Rypdal, M., Scotto,
M. G., Vannitsem, S., Watkins, N. W., Yang, L. C., and Yuan,
N. M.: The Structure of Climate Variability Across Scales, Rev.
Geophys., 58, 44, https://doi.org/10.1029/2019rg000657, 2020.
Gettelman, A., Morrison, H., Santos, S., Bogenschutz, P.,
and Caldwell, P. M.: Advanced Two-Moment Bulk Microphysics for Global Models. Part II: Global Model Solutions
and Aerosol-Cloud Interactions, J. Climate, 28, 1288–1307,
https://doi.org/10.1175/jcli-d-14-00103.1, 2015.
Gettelman, A., Hannay, C., Bacmeister, J. T., Neale, R. B., Pendergrass, A. G., Danabasoglu, G., Lamarque, J. F., Fasullo, J.
T., Bailey, D. A., Lawrence, D. M., and Mills, M. J.: HighClimate Sensitivity in the Community Earth System Model
Version 2 (CESM2), Geophys. Res. Lett., 46, 8329–8337,
https://doi.org/10.1029/2019gl083978, 2019.
Ghil, M.: A mathematical theory of climate sensitivity or, How todeal with both anthropogenic forcing and natural variability?, in:Climate change: Multidecadal and beyond, edited by: Chang, C.
P.,Ghil, M., Latif, M., and Wallace, J. M., World Scientific Publ.
Co./Imperial College Press, London, 31–51, 2015.
Goosse, H., Kay, J. E., Armour, K. C., Bodas-Salcedo, A., Chepfer, H., Docquier, D., Jonko, A., Kushner, P. J., Lecomte, O.,
Massonnet, F., Park, H. S., Pithan, F., Svensson, G., and Vancoppenolle, M.: Quantifying climate feedbacks in polar regions, Nat.
Commun., 9, 1919, https://doi.org/10.1038/s41467-018-041730, 2018.
Hasselmann, K.: Stochastic Climate Models. 1. Theory, Tellus, 28,
473–485, https://doi.org/10.1111/j.2153-3490.1976.tb00696.x,
1976.
Haszpra, T., Herein, M., and Bódai, T.: Investigating ENSO
and its teleconnections under climate change in an ensemble

Earth Syst. Dynam., 12, 1393–1411, 2021

1408

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

view – a new perspective, Earth Syst. Dynam., 11, 267–280,
https://doi.org/10.5194/esd-11-267-2020, 2020.
Hegerl, G. C., Zwiers, F. W., Braconnot, P., Gillett, N. P., Luo, Y.
M., Marengo Orsini, J. A., Nicholls, N., Penner, J. E., and Stott,
P. A.: Understanding and Attributing Climate Change, in: Climate Change 2007: The Physical Science Basis. Contribution ofWorking Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Solomon, S.,
Qin, D., Manning, M., Chen, Z., Marquis, M., Averyt, K. B., Tignor, M., and Miller, H. L., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 2007.
Holmes, C. R., Woollings, T., Hawkins, E., and de Vries, H.: Robust Future Changes in Temperature Variability under Greenhouse Gas Forcing and the Relationship with Thermal Advection, J. Climate, 29, 2221–2236, https://doi.org/10.1175/jcli-d14-00735.1, 2016.
Huntingford, C., Jones, P. D., Livina, V. N., Lenton, T. M.,
and Cox, P. M.: No increase in global temperature variability despite changing regional patterns, Nature, 500, 327–330,
https://doi.org/10.1038/nature12310, 2013.
Ishii, M., Fukuda, Y., Hirahara, S., Yasui, S., Suzuki, T., and Sato,K.: Accuracy of Global Upper Ocean Heat Content Estimation
Expected from Present Observational Data Sets, Sola, 13, 163–
167, https://doi.org/10.2151/sola.2017-030, 2017.
Jin, F. F., Lin, L., Timmermann, A., and Zhao, J.: Ensemblemean dynamics of the ENSO recharge oscillator under statedependent stochastic forcing, Geophys. Res. Lett., 34, L03807,
https://doi.org/10.1029/2006gl027372, 2007.
Jin, F. F., Chen, H. C., Zhao, S., Hayashi, M., Karamperidou, C.,
Stuecker, M. F., Xie, R., and Geng, L.: Simple ENSO models, in:
El Niño Southern Oscillation in a Changing Climate, edited by:
McPhaden, M. J., Santoso, A., and Cai, W., Geophysical Monograph Series, AGU, 119–151, 2020.
Kay, J. E., Deser, C., Phillips, A., Mai, A., Hannay, C., Strand,G., Arblaster, J. M., Bates, S. C., Danabasoglu, G., Edwards,
J., Holland, M., Kushner, P., Lamarque, J. F., Lawrence, D.,
Lindsay, K., Middleton, A., Munoz, E., Neale, R., Oleson, K.,
Polvani, L., and Vertenstein, M.: The Community Earth System
Model (CESM) Large Ensemble Project: A Community Resource for Studying Climate Change in the Presence of Internal Climate Variability, B. Am. Meteorol. Soc., 96, 1333–1349,
https://doi.org/10.1175/bams-d-13-00255.1, 2015.
Kim, S. T., Cai, W. J., Jin, F. F., Santoso, A., Wu, L. X., Guilyardi,
E., and An, S. I.: Response of El Niño sea surface temperature
variability to greenhouse warming, Nat. Clim. Change, 4, 786–
790, https://doi.org/10.1038/nclimate2326, 2014.
Kwiatkowski, L., Torres, O., Bopp, L., Aumont, O., Chamberlain,M., Christian, J. R., Dunne, J. P., Gehlen, M., Ilyina, T., John, J.
G., Lenton, A., Li, H., Lovenduski, N. S., Orr, J. C., Palmieri,
J., Santana-Falcón, Y., Schwinger, J., Séférian, R., Stock, C.
A., Tagliabue, A., Takano, Y., Tjiputra, J., Toyama, K., Tsujino, H., Watanabe, M., Yamamoto, A., Yool, A., and Ziehn, T.:
Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline
from CMIP6 model projections, Biogeosciences, 17, 3439–3470,
https://doi.org/10.5194/bg-17-3439-2020, 2020.
Lawrence, D. M., Fisher, R. A., Koven, C. D., Oleson, K. W., Swenson, S. C., Bonan, G., Collier, N., Ghimire, B., van Kampenhout, L., Kennedy, D., Kluzek, E., Lawrence, P. J., Li, F., Li,

Earth Syst. Dynam., 12, 1393–1411, 2021

H. Y., Lombardozzi, D., Riley, W. J., Sacks, W. J., Shi, M. J.,
Vertenstein, M., Wieder, W. R., Xu, C. G., Ali, A. A., Badger,
A. M., Bisht, G., van den Broeke, M., Brunke, M. A., Burns,
S. P., Buzan, J., Clark, M., Craig, A., Dahlin, K., Drewniak, B.,
Fisher, J. B., Flanner, M., Fox, A. M., Gentine, P., Hoffman,
F., Keppel-Aleks, G., Knox, R., Kumar, S., Lenaerts, J., Leung,
L. R., Lipscomb, W. H., Lu, Y. Q., Pandey, A., Pelletier, J. D.,
Perket, J., Randerson, J. T., Ricciuto, D. M., Sanderson, B. M.,
Slater, A., Subin, Z. M., Tang, J. Y., Thomas, R. Q., Martin,
M. V., and Zeng, X. B.: The Community Land Model Version
5: Description of New Features, Benchmarking, and Impact ofForcing Uncertainty, J. Adv. Model. Earth Sy., 11, 4245–4287,
https://doi.org/10.1029/2018ms001583, 2019.
Levine, A. F. Z. and Jin, F. F.: Noise-Induced Instability inthe ENSO Recharge Oscillator, J. Atmos. Sci., 67, 529–542,
https://doi.org/10.1175/2009jas3213.1, 2010.
Li, F. and Lawrence, D. M.: Role of Fire in the Global Land Water
Budget during the Twentieth Century due to Changing Ecosystems, J. Climate, 30, 1893–1908, https://doi.org/10.1175/jcli-d16-0460.1, 2017.
Li, F., Levis, S., and Ward, D. S.: Quantifying the role of fire in
the Earth system – Part 1: Improved global fire modeling in theCommunity Earth System Model (CESM1), Biogeosciences, 10,2293–2314, https://doi.org/10.5194/bg-10-2293-2013, 2013.
Li, Q., Webb, A., Fox-Kemper, B., Craig, A., Danabasoglu, G.,
Large, W. G., and Vertenstein, M.: Langmuir mixing effects on
global climate: WAVEWATCH III in CESM, Ocean Model., 103,
145–160, https://doi.org/10.1016/j.ocemod.2015.07.020, 2016.
Loeb, N. G., Wielicki, B. A., Doelling, D. R., Smith, G. L., Keyes,
D. F., Kato, S., Manalo-Smith, N., and Wong, T.: Toward Optimal
Closure of the Earth’s Top-of-Atmosphere Radiation Budget,
J. Climate, 22, 748–766, https://doi.org/10.1175/2008jcli2637.1,
2009.
Loeb, N. G., Doelling, D. R., Wang, H. L., Su, W. Y., Nguyen,
C., Corbett, J. G., Liang, L. S., Mitrescu, C., Rose, F. G.,
and Kato, S.: Clouds and the Earth’s Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) Top-ofAtmosphere (TOA) Edition-4.0 Data Product, J. Climate, 31,
895–918, https://doi.org/10.1175/jcli-d-17-0208.1, 2018.
Lombardozzi, D. L., Lu, Y. Q., Lawrence, P. J., Lawrence,
D. M., Swenson, S., Oleson, K. W., Wieder, W. R., andAinsworth, E. A.: Simulating Agriculture in the Community
Land Model Version 5, J. Geophys. Res.-Biogeo., 125, 19,https://doi.org/10.1029/2019jg005529, 2020.
Long, M. C., Moore, J. K., Lindsay, K., Levy, M. N., Doney, S.
C., Luo, J. Y., Krumhardt, K. M., Letscher, R. T., Grover, M.,
and Sylvester, Z. T.: Simulations with the Marine Biogeochemistry Library (MARBL), J. Adv. Model. Earth Syst., accepted,
https://doi.org/10.1029/2021MS002647, 2021.
Maher, N., Matei, D., Milinski, S., and Marotzke, J.: ENSO
Change in Climate Projections: Forced Response or Internal Variability?, Geophys. Res. Lett., 45, 11390–11398,
https://doi.org/10.1029/2018gl079764, 2018.
Maher, N., Milinski, S., Suarez-Gutierrez, L., Botzet, M., Dobrynin,M., Kornblueh, L., Kroger, J., Takano, Y., Ghosh, R., Hedemann, C., Li, C., Li, H. M., Manzini, E., Notz, D., Putrasahan, D., Boysen, L., Claussen, M., Ilyina, T., Olonscheck, D.,
Raddatz, T., Stevens, B., and Marotzke, J.: The Max Planck Institute Grand Ensemble: Enabling the Exploration of Climate

https://doi.org/10.5194/esd-12-1393-2021

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability
System Variability, J. Adv. Model. Earth Sy., 11, 2050–2069,
https://doi.org/10.1029/2019ms001639, 2019.
Maher, N., Power, S. B., and Marotzke, J.: More accurate quantification of model-to-model agreement in externally forced climatic responses over the coming century, Nat. Commun., 12, 13,
https://doi.org/10.1038/s41467-020-20635-w, 2021.
Majda, A. J., Franzke, C., and Crommelin, D.: Normal forms for
reduced stochastic climate models, P. Natl. Acad. Sci. USA, 106,3649–3653, https://doi.org/10.1073/pnas.0900173106, 2009.
Menary, M. B. and Wood, R. A.: An anatomy of the projected North
Atlantic warming hole in CMIP5 models, Clim. Dynam., 50,
3063–3080, https://doi.org/10.1007/s00382-017-3793-8, 2018.
Milinski, S., Maher, N., and Olonscheck, D.: How large does alarge ensemble need to be?, Earth Syst. Dynam., 11, 885–901,
https://doi.org/10.5194/esd-11-885-2020, 2020.
Moore, J. K., Doney, S. C., Kleypas, J. A., Glover, D. M., andFung, I. Y.: An intermediate complexity marine ecosystem modelfor the global domain, Deep-Sea Res. Pt. II, 49, 403–462,
https://doi.org/10.1016/s0967-0645(01)00108-4, 2001.
Moore, J. K., Doney, S. C., and Lindsay, K.: Upper ocean
ecosystem dynamics and iron cycling in a global threedimensional model, Global Biogeochem. Cy., 18, 21,
https://doi.org/10.1029/2004gb002220, 2004.
Moore, J. K., Lindsay, K., Doney, S. C., Long, M. C.,
and Misumi, K.: Marine Ecosystem Dynamics and Biogeochemical Cycling in the Community Earth System Model
CESM1(BGC): Comparison of the 1990s with the 2090s under
the RCP4.5 and RCP8.5 Scenarios, J. Climate, 26, 9291–9312,
https://doi.org/10.1175/jcli-d-12-00566.1, 2013.
Morice, C. P., Kennedy, J. J., Rayner, N. A., and Jones, P.
D.: Quantifying uncertainties in global and regional temperature change using an ensemble of observational estimates:
The HadCRUT4 data set, J. Geophys. Res.-Atmos., 117, 22,https://doi.org/10.1029/2011jd017187, 2012.
Müller, D.: Bispectra of Sea-Surface Temperature Anomalies,
J. Phys. Oceanogr., 17, 26–36, https://doi.org/10.1175/15200485(1987)017<0026:bossta>2.0.co;2, 1987.
Myers-Smith, I. H., Kerby, J. T., Phoenix, G. K., Bjerke, J. W., Epstein, H. E., Assmann, J. J., John, C., Andreu-Hayles, L., AngersBlondin, S., and Beck, P. S.: Complexity revealed in the greeningof the Arctic, Nat. Clim. Change, 10, 106–117, 2020.
O’Neill, B. C., Tebaldi, C., van Vuuren, D. P., Eyring, V., Friedlingstein, P., Hurtt, G., Knutti, R., Kriegler, E., Lamarque, J.-F.,
Lowe, J., Meehl, G. A., Moss, R., Riahi, K., and Sanderson, B. M.: The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6, Geosci. Model Dev., 9, 3461–3482,
https://doi.org/10.5194/gmd-9-3461-2016, 2016.
Pendergrass, A. G., Knutti, R., Lehner, F., Deser, C., and Sanderson, B. M.: Precipitation variability increases in a warmer climate, Sci. Rep., 7, 9, https://doi.org/10.1038/s41598-017-17966y, 2017.
Phillips, A. S., Deser, C., Fasullo, J., Schneider, D. P., andSimpson, I. R.: Assessing Climate Variability and Change inModel Large Ensembles: A User’s Guide to the “Climate Variability Diagnostics Package for Large Ensembles”, version 1,
https://doi.org/10.5065/h7c7-f961, 2020.
Rahmstorf, S., Box, J. E., Feulner, G., Mann, M. E., Robinson, A., Rutherford, S., and Schaffernicht, E. J.: Exceptional twentieth-century slowdown in Atlantic Ocean

https://doi.org/10.5194/esd-12-1393-2021

1409

overturning circulation, Nat. Clim. Change, 5, 475–480,
https://doi.org/10.1038/nclimate2554, 2015.
Raisanen, J.: CO2 -induced changes in interannual temperature and precipitation variability in 19 CMIP2 experiments, J. Climate, 15, 2395–2411, https://doi.org/10.1175/15200442(2002)015<2395:ciciit>2.0.co;2, 2002.
Renner, S. S. and Zohner, C. M.: Climate Change and Phenological Mismatch in Trophic Interactions Among Plants, Insects, and Vertebrates, Annu. Rev. Ecol. Evol. S., 49, 165–182,
https://doi.org/10.1146/annurev-ecolsys-110617-062535, 2018.
Revelle, R. and Suess, H. E.: Carbon Dioxide Exchange Between
Atmosphere and Ocean and the Question of an Increase of Atmospheric CO2 During the Past Two Decades, Tellus, 9, 18–27,
1957.
Rind, D., Goldberg, R., and Ruedy, R.: Change in Climate
Variability in the 21st-Century, Climatic Change, 14, 5–37,
https://doi.org/10.1007/bf00140173, 1989.
Robin, Y., Yiou, P., and Naveau, P.: Detecting changes inforced climate attractors with Wasserstein distance, Nonlin. Processes Geophys., 24, 393–405, https://doi.org/10.5194/npg-24393-2017, 2017.
Rodgers, K. B., Schlunegger, S., Slater, R. D., Ishii, M.,
Frolicher, T. L., Toyama, K., Plancherel, Y., Aumont, O.,
and Fassbender, A. J.: Reemergence of Anthropogenic Carbon Into the Ocean’s Mixed Layer Strongly Amplifies
Transient Climate Sensitivity, Geophys. Res. Lett., 47, 9,
https://doi.org/10.1029/2020gl089275, 2020.
Sarachik, E. S.: The tools of adaptation, Second International Conference on Climate, Sustainability, and Development in SemiArid Regions, Fortaleza, Brazil, 2010.
Sardeshmukh, P. D. and Penland, C.: Understanding the distinctively skewed and heavy tailed character of atmosphericand oceanic probability distributions, Chaos, 25, 036410,
https://doi.org/10.1063/1.4914169, 2015.
Sardeshmukh, P. D. and Sura, P.: Reconciling Non-Gaussian Climate Statistics with Linear Dynamics, J. Climate, 22, 1193–
1207, https://doi.org/10.1175/2008jcli2358.1, 2009.
Schlunegger, S., Rodgers, K. B., Sarmiento, J. L., Ilyina,
T., Dunne, J. P., Takano, Y., Christian, J. R., Long, M.
C., Frolicher, T. L., Slater, R., and Lehner, F.: Time of
Emergence and Large Ensemble Intercomparison for Ocean
Biogeochemical Trends, Global Biogeochem. Cy., 34, 18,
https://doi.org/10.1029/2019gb006453, 2020.
Schneider, T., Bischoff, T., and Plotka, H.: Physics of Changes
in Synoptic Midlatitude Temperature Variability, J. Climate, 28,
2312–2331, https://doi.org/10.1175/jcli-d-14-00632.1, 2015.
Screen, J. A.: Arctic amplification decreases temperature variancein northern mid- to high-latitudes, Nat. Clim. Change, 4, 577–
582, https://doi.org/10.1038/nclimate2268, 2014.
Screen, J. A., Deser, C., and Sun, L.: Reduced risk of North American cold extremes due to continued Arctic sea ice loss, B. Am.
Meteorol. Soc., 96, 1489–1503, https://doi.org/10.1175/bams-d14-00185.1, 2015.
Simpson, I. R., Bacmeister, J., Neale, R. B., Hannay, C., Gettelman,
A., Garcia, R. R., Lauritzen, P. H., Marsh, D. R., Mills, M. J.,
Medeiros, B., and Richter, J. H.: An Evaluation of the LargeScale Atmospheric Circulation and Its Variability in CESM2and Other CMIP Models, J. Geophys. Res.-Atmos., 125, 42,
https://doi.org/10.1029/2020jd032835, 2020.

Earth Syst. Dynam., 12, 1393–1411, 2021

1410

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability

Smith, R., Jones, P., Briegleb, B., Bryan, F., Danabasoglu, G.,
Dennis, J., Dukowicz, J., Eden, C., Fox-Kemper, B., Gent, P.,
Hecht, M., Jayne, S., Jochum, M., Large, W., Lindsay, K., Maltrud, M., Norton, N., Peacock, S., Vertenstein, M., and Yeager, S.: The Parallel Ocean Program (POP) Reference Manual,
Ocean Component of the Community Climate System Model(CCSM), available at: https://opensky.ucar.edu/islandora/object/manuscripts:825 (last access: 19 November 2021), 2010.
Stein, K., Schneider, N., Timmermann, A., and Jin, F.
F.: Seasonal Synchronization of ENSO Events in a
Linear Stochastic Model, J. Climate, 23, 5629–5643,
https://doi.org/10.1175/2010jcli3292.1, 2010.
Stein, K., Timmermann, A., Schneider, N., Jin, F. F., and Stuecker,
M. F.: ENSO Seasonal Synchronization Theory, J. Climate, 27,
5285–5310, https://doi.org/10.1175/jcli-d-13-00525.1, 2014.
Stocker, T., Intergovernmental Panel on Climate Change, W. G. I.,
WMO, and UNEP: Climate Change 2013: The Physical Science
Basis: Working Group I Contribution to the 5th Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, England, 2013.
Stouffer, R. J. and Wetherald, R. T.: Changes of variability in response to increasing greenhouse gases. Part I: Temperature, J.
Climate, 20, 5455–5467, https://doi.org/10.1175/2007jcli1384.1,
2007.
Stuecker, M. F., Timmermann, A., Jin, F. F., McGregor, S., andRen, H. L.: A combination mode of the annual cycle andthe El Niño/Southern Oscillation, Nat. Geosci., 6, 540–544,
https://doi.org/10.1038/ngeo1826, 2013.
Stuecker, M. F., Jin, F. F., and Timmermann, A.: El Niño-Southern
Oscillation frequency cascade, P. Natl. Acad. Sci. USA, 112,13490–13495, https://doi.org/10.1073/pnas.1508622112, 2015a.
Stuecker, M. F., Jin, F. F., Timmermann, A., and McGregor, S.: Combination Mode Dynamics of the Anomalous
Northwest Pacific Anticyclone, J. Climate, 28, 1093–1111,
https://doi.org/10.1175/jcli-d-14-00225.1, 2015b.
Stuecker, M. F., Timmermann, A., Jin, F. F., Proistosescu, C., Kang,
S. M., Kim, D., Yun, K. S., Chung, E. S., Chu, J. E., Bitz, C. M.,
Armour, K. C., and Hayashi, M.: Strong remote control of future
equatorial warming by off-equatorial forcing, Nat. Clim. Change,
10, 124–129, https://doi.org/10.1038/s41558-019-0667-6, 2020.
Sun, L. T., Deser, C., and Tomas, R. A.: Mechanisms ofStratospheric and Tropospheric Circulation Response to Projected Arctic Sea Ice Loss, J. Climate, 28, 7824–7845,
https://doi.org/10.1175/jcli-d-15-0169.1, 2015.
Swain, D. L., Langenbrunner, B., Neelin, J. D., and Hall,
A.: Increasing precipitation volatility in twenty-firstcentury California, Nat. Clim. Change, 8, 427–433,
https://doi.org/10.1038/s41558-018-0140-y, 2018.
Tamarin-Brodsky, T., Hodges, K., Hoskins, B. J., and Shepherd,
T. G.: Changes in Northern Hemisphere temperature variability
shaped by regional warming patterns, Nat. Geosci., 13, 414–421,
https://doi.org/10.1038/s41561-020-0576-3, 2020.
Taschetto, A. S., Ummenhofer, C. C., Stuecker, M. F., Dommenget,
D., Ashok, K., Rodrigues, R. R., and Yeh, S. W.: ENSO atmospheric teleconnections, in: El Niño Southern Oscillation in a
Changing Climate, edited by: McPhaden, M. J., Santoso, A., andCai, W., Geophysical Monograph Series, AGU, 309–335, 2020.
Timmermann, A., Oberhuber, J., Bacher, A., Esch, M., Latif, M.,
and Roeckner, E.: Increased El Niño frequency in a climate

Earth Syst. Dynam., 12, 1393–1411, 2021

model forced by future greenhouse warming, Nature, 398, 694–
697, 1999.
Timmermann, A., Lorenz, S. J., An, S. I., Clement, A., and Xie,
S. P.: The effect of orbital forcing on the mean climate andvariability of the tropical Pacific, J. Climate, 20, 4147–4159,
https://doi.org/10.1175/jcli4240.1, 2007.
Timmermann, A., Sachs, J., and Timm, O. E.: Assessing divergent
SST behavior during the last 21 ka derived from alkenones andG. ruber-Mg/Ca in the equatorial Pacific, Paleoceanography, 29,
680–696, https://doi.org/10.1002/2013pa002598, 2014.
Torrence,C.
andCompo,G.
P.:
A
practical
guidetowaveletanalysis,
B.
Am.
Meteorol.
Soc.,
79,
61–78,
https://doi.org/10.1175/15200477(1998)079<0061:apgtwa>2.0.co;2, 1998.
Toyama, K., Rodgers, K. B., Blanke, B., Iudicone, D., Ishii, M., Aumont, O., and Sarmiento, J. L.: Large Reemergence of Anthropogenic Carbon into the Ocean’s Surface Mixed Layer Sustainedby the Ocean’s Overturning Circulation, J. Climate, 30, 8615–
8631, https://doi.org/10.1175/jcli-d-16-0725.1, 2017.
van Marle, M. J. E., Kloster, S., Magi, B. I., Marlon, J. R., Daniau, A.-L., Field, R. D., Arneth, A., Forrest, M., Hantson,
S., Kehrwald, N. M., Knorr, W., Lasslop, G., Li, F., Mangeon, S., Yue, C., Kaiser, J. W., and van der Werf, G. R.: Historic global biomass burning emissions for CMIP6 (BB4CMIP)
based on merging satellite observations with proxies and firemodels (1750–2015), Geosci. Model Dev., 10, 3329–3357,
https://doi.org/10.5194/gmd-10-3329-2017, 2017.
Vissio, G., Lembo, V., Lucarini, V., and Ghil, M.: Evaluating the performance of climate models based on Wasserstein distance, Geophys. Res. Lett., 47, e2020GL089385,
https://doi.org/10.1029/2020GL089385, 2020.
Weijer, W., Cheng, W., Garuba, O. A., Hu, A., and Nadiga, B. T.:
CMIP6 models predict significant 21st century decline of the Atlantic Meridional Overturning Circulation, Geophys. Res. Lett.,
47, e2019GL086075, https://doi.org/10.1029/2019GL086075,
2020.
Wetherald, R. T.: Changes of Variability in Response to Increasing Greenhouse Gases. Part II: Hydrology, J. Climate, 22, 6089–
6103, https://doi.org/10.1175/2009jcli2834.1, 2009.
Williams, I. N. and Patricola, C. M.: Diversity of ENSO Events
Unified by Convective Threshold Sea Surface Temperature: A
Nonlinear ENSO Index, Geophys. Res. Lett., 45, 9236–9244,
https://doi.org/10.1029/2018gl079203, 2018.
Xie, S. P.: On the Genesis of the Equatorial Annual Cycle, J. Climate, 7, 2008–2013, https://doi.org/10.1175/15200442(1994)007<2008:otgote>2.0.co;2, 1994.
Xie, S. P.: Effects of seasonal solar forcing on
latitudinal
asymmetry
of
the
ITCZ,
J.
Climate,9,2945–2950,
https://doi.org/10.1175/15200442(1996)009<2945:eossfo>2.0.co;2, 1996.
Xie, S. P., Deser, C., Vecchi, G. A., Ma, J., Teng, H. Y., andWittenberg, A. T.: Global Warming Pattern Formation: Sea
Surface Temperature and Rainfall, J. Climate, 23, 966–986,
https://doi.org/10.1175/2009jcli3329.1, 2010.
Yun, K. S., Lee, J. Y., Timmermann, A., Stein, K., Stuecker, J. F.,
Fyfe, J. C., and Chung, E. S.: Increasing ENSO–rainfall variability due to changes in future tropical temperature–rainfall relationship, Commun. Earth Environ., 2, 1–7, 2021.

https://doi.org/10.5194/esd-12-1393-2021

K. B. Rodgers et al.: Ubiquity of human-induced changes in climate variability
Zelle, H., Van Oldenborgh, G. J., Burgers, G., and Dijkstra, H.:
El Niño and Greenhouse warming: Results from ensemble simulations with the NCAR CCSM, J. Climate, 18, 4669–4683,
https://doi.org/10.1175/jcli3574.1, 2005.

https://doi.org/10.5194/esd-12-1393-2021

1411

Zhu, Z. C., Piao, S. L., Myneni, R. B., Huang, M. T., Zeng, Z. Z.,
Canadell, J. G., Ciais, P., Sitch, S., Friedlingstein, P., Arneth, A.,
Cao, C. X., Cheng, L., Kato, E., Koven, C., Li, Y., Lian, X., Liu,
Y. W., Liu, R. G., Mao, J. F., Pan, Y. Z., Peng, S. S., Penuelas,J., Poulter, B., Pugh, T. A. M., Stocker, B. D., Viovy, N., Wang,
X. H., Wang, Y. P., Xiao, Z. Q., Yang, H., Zaehle, S., and Zeng,
N.: Greening of the Earth and its drivers, Nat. Clim. Change, 6,
791–795, https://doi.org/10.1038/nclimate3004, 2016.

Earth Syst. Dynam., 12, 1393–1411, 2021