Interannual and Seasonal Drivers of Carbon Cycle Variability Represented by the Community Earth System Model (CESM2).

Earth system models are intended to make long-term projections, but they can be evaluated at interannual and seasonal time scales. Although the Community Earth System Model (CESM2) showed improvements in a number of terrestrial carbon cycle benchmarks, relative to its predecessor, our analysis suggests that the interannual variability (IAV) in net terrestrial carbon fluxes did not show similar improvements. The model simulated low IAV of net ecosystem production (NEP), resulting in a weaker than observed sensitivity of the carbon cycle to climate variability. Low IAV in net fluxes likely resulted from low variability in gross primary productivity (GPP)-especially in the tropics-and a high covariation between GPP and ecosystem respiration. Although lower than observed, the IAV of NEP had significant climate sensitivities, with positive NEP anomalies associated with warmer and drier conditions in high latitudes, and with wetter and cooler conditions in mid and low latitudes. We identified two dominant modes of seasonal variability in carbon cycle flux anomalies in our fully coupled CESM2 simulations that are characterized by seasonal amplification and redistribution of ecosystem fluxes. Seasonal amplification of net and gross carbon fluxes showed climate sensitivities mirroring those of annual fluxes. Seasonal redistribution of carbon fluxes is initiated by springtime temperature anomalies, but subsequently negative feedbacks in soil moisture during the summer and fall result in net annual carbon losses from land. These modes of variability are also seen in satellite proxies of GPP, suggesting that CESM2 appropriately represents regional sensitivities of photosynthesis to climate variability on seasonal time scales.

Introduction

Terrestrial ecosystems continue to provide a sink for about a quarter of anthropogenic carbon dioxide (CO2) emissions (Ballantyne et al., 2012; Friedlingstein et al., 2019), but the long‐term strength and locations of this sink remain uncertain (Gaubert et al., 2019; Schimel, Stephens & Fisher; 2015; Tagesson et al., 2020). The net terrestrial flux of CO2 to or from the atmosphere depends on the balance between much larger carbon fluxes that are driven by plant productivity, ecosystem respiration, and disturbances like fire (Keppel‐Aleks et al., 2014). Observational constraints on these gross fluxes are difficult to make globally, which results in persistently high uncertainty in terrestrial carbon cycle projections (Anav et al., 2013; Arora et al., 2020; Friedlingstein et al., 2014; Liddicoat et al., 2021). Thus, capturing the appropriate carbon cycle sensitivities to climate driven variability at interannual and seasonal time scales over the observational record may be important to improving longer term projections of terrestrial carbon balance.

At decadal‐to century‐time scales, the net exchange of CO2 between the land and atmosphere remains one of the more robust benchmarks by which to evaluate the representation of terrestrial biogeochemistry in land models (Collier et al., 2018; Hoffman et al., 2014; Keppel‐Aleks et al., 2013). Separately, the trends in the seasonal amplitude of terrestrial CO2 fluxes provides insights to how global change drivers are altering the net exchange of terrestrial CO2 with the atmosphere (Graven et al., 2013; Ito et al., 2016; Piao et al., 2018). Indeed, successive generations of the Community Earth System Model (CESM) and its terrestrial component, the Community Land Model (CLM) show improvements in the globally integrated net terrestrial carbon flux over the historical period (Bonan et al., 2019; Danabasoglu et al., 2020; Lawrence et al., 2019). The model also shows improved seasonal amplitude of terrestrial CO2 fluxes and better represents their increasing trend over the observational record (Lawrence et al., 2019; Lombardozzi et al., 2020). This suggests that the model adequately represents dominant features influencing longer term trends in the terrestrial carbon cycle dynamics, namely land‐use and land cover change as well as potential CO2 fertilization effects (Wieder et al., 2019). These longer term benchmarks, however, offer little insight into the environmental sensitivities of terrestrial CO2 fluxes at annual to seasonal time scales, which are important for understanding carbon cycle sensitivities to climate variability and potential carbon cycle responses to future climate change.

Measurements of the atmospheric CO2 growth rate provide an integrated estimate of the interannual variability (IAV) in global carbon cycle (Keeling et al., 1995; Zeng et al., 2005). Since much of the observed IAV is driven by terrestrial processes (Bousquet et al., 2000), variability in the atmospheric CO2 record provides a top‐town constraint on the climate sensitivity of land‐atmospheric CO2 exchange (reviewed by Piao et al., 2020). Specifically, natural climate variability in temperature and precipitation over land are the primary drivers of terrestrial carbon cycle variability that can be inferred from IAV in the atmospheric CO2 growth rate. This includes the temperature sensitivity of gross primary productivity (GPP) and ecosystem respiration (Reco) on net ecosystem exchange (NEE, calculated as GPP ‐ Reco), especially in the tropics (Anderegg et al., 2015; Ballantyne et al., 2017; Cox et al., 2013; Liu et al., 2017; Rödenbeck et al., 2018b; Wang et al., 2013). Meanwhile, other studies emphasize the importance of soil moisture variability on GPP and NEE, especially in arid and semi‐arid ecosystems (Anderegg et al., 2015; Humphrey et al., 2018; Poulter et al., 2014). Jung et al. (2017) suggest that compensating moisture‐driven variation in gross, local‐scale fluxes (GPP and Reco), as well as spatial compensation in moisture anomalies among regions leaves a dominant temperature‐driven signal in the IAV of land‐atmosphere CO2 exchange. By using a mechanistic model, we aim to provide simultaneous insights into the drivers of carbon cycle variability across a range of spatial scales, including finer spatial and temporal resolutions.

Disentangling the contributions of gross carbon fluxes and their impact on the IAV of NEE remains a challenge. While GPP and Reco anomalies are strongly correlated, mounting evidence suggests that variance in NEE is more strongly correlated with variance in GPP than Reco (Baldocchi et al., 2018; Schwalm et al., 2010). Globally gridded estimates of net and gross carbon fluxes are derived from upscaled flux tower measurements (Beer et al., 2010; Jung et al., 2011, 2017) or remote sensing (Alemohammad et al., 2017; Köhler et al., 2018; Running et al., 2004), but both of these tend to underestimate carbon cycle IAV (Jung et al., 2020; O'Sullivan et al., 2020; Piao et al., 2020; Y. Zhang et al., 2018). Remote sensing products offer promise for diagnosing seasonal productivity responses to climate variability in the tropics and extra‐tropics (Bowman et al., 2017; Buermann et al., 2018; Li & Xiao, 2019; Liu et al., 2017; Schimel, Pavlick et al., 2015). Notably, Butterfield et al. (2020) found that while local to regional‐scale IAV in productivity was poorly correlated among remote sensing products (between 2007 and 2015), they could identify seasonal modes of variability that shared common features and environmental sensitivities. These included the (a) Amplification of the seasonal cycle of GPP, which was associated with increases in summertime soil moisture availability and positive anomalies in GPP, and (b) Seasonal redistribution of GPP that was initially driven by warmer springtime temperatures—resulting in positive GPP anomalies—but that was followed by higher‐than‐average soil moisture stress in the summer and fall—resulting in negative GPP anomalies (see also Byrne et al., 2020). Integrated over the growing season, the second mode of variability reflects climate driven shifts in the seasonal timing of photosynthesis, which has negligible changes in the annual C flux.

Given improvements to the representation of terrestrial biogeochemistry in CESM2, relative to CESM1 (Bonan et al., 2019; Danabasoglu et al., 2020; Lawrence et al., 2019; Lombardozzi et al., 2020), we investigated the representation of carbon cycle variability at interannual and seasonal timescales in the model. Part of this work was motivated by findings of Lawrence et al. (2019), who noted that metrics for terrestrial carbon cycle IAV were degraded in CLM5, especially in the tropics, relative to previous versions of the model. Second, given new observational analysis that identified dominant modes of seasonal variability in regional productivity over the satellite record (Butterfield et al., 2020; Byrne et al., 2020) we explored the predictive model's ability to simulate similar behavior in global‐scale simulations. Thus, our aims were to (a) quantify the IAV of land carbon fluxes simulated over the historical period by CESM2‐esm‐hist and their sensitivity to climate variability; (b) identify modes of seasonal variability in simulated carbon fluxes and their likely climate drivers.

Methods

Model Simulations

We analyzed simulations from the Community Earth System Model, version 2 (CESM2) that couples atmosphere, ocean, land, sea ice, land ice, and river transport components to simulate physical and biogeochemical conditions over historical and future scenarios (Danabasoglu et al., 2020). Of greatest importance for the simulations analyzed here are the atmospheric and land components, which are each briefly described below. The atmosphere model in CESM2 is the Community Atmosphere Model, version 6 (CAM6) which applies the same Finite Volume dynamical core (Danabasoglu et al., 2020). Relative to previous versions of the model, CESM2 shows an improvement in its representation of El Niño Southern Oscillation (ENSO) events and their effect on precipitation and temperature anomaly patterns (Meehl et al., 2020). The atmosphere and land models are run at a nominal 1° horizontal resolution (1.25° longitude by 0.9° latitude) and are coupled every 30 min.

The terrestrial component of CESM2 uses the Community Land Model, version 5 (CLM5; Lawrence et al., 2019). Briefly, these developments simulate transient agricultural expansion and land management (Lombardozzi et al., 2020), plant hydraulic stress (Kennedy et al., 2019), and improve the representation of plant nitrogen limitation (R. A. Fisher et al., 2019; Wieder et al., 2019). We used the historical simulation of CESM2‐esm‐hist from 1850 to 2014 that has active biogeochemical representation of terrestrial carbon and nitrogen cycles and that also simulates the prognostic evolution of atmospheric CO2 concentrations based on fluxes with oceans and terrestrial ecosystems.

As described in Danabasoglu et al. (2020), initial conditions for the land model and ocean model biogeochemical tracers in the non‐ESM piControl experiment were generated using spin‐up runs of the land and ocean models, respectively. In these spin‐up runs, the active atmospheric component was replaced with a data atmosphere that repeatedly cycled through 21 years of surface forcing that were extracted from a fully coupled CESM2 experiment. Twenty‐one years of forcing were used to capture some aspects of IAV. The land model spin‐up consisted of an accelerated decomposition (AD) mode segment and a subsequent synchronous spin‐up segment. These segments were run for 252 and 1,701 years respectively. The ocean model spin‐up was applied to biogeochemical tracers. The ocean spin‐up was run for 1,029 years and also utilized a Newton‐Krylov solver, based on Lindsay (2017) to more completely spin up a subset of the biogeochemical tracers, including the carbon pools. The esm‐piControl was initialized from the piControl experiment using an incremental coupling approach. In an intermediate experiment, which was initialized from the piControl experiment, the carbon cycle of the surface components was coupled bidirectionally to a CO2 tracer in the atmospheric model. This intermediate experiment was run for 80 years, during which the surface biogeochemical parameterizations adjusted to the prognostic atmospheric CO2. The esm‐piControl experiment was initialized with the model state from the end of this intermediate experiment, and the prognostic CO2 was coupled to the radiative computations in the atmospheric model. The esm‐hist experiments analyzed here were initialized from the esm‐piControl experiment. The CESM2‐esm historical simulations used CMIP6 forcings for anthropogenic emissions, biomass burning, and volcanic SO2 emissions from 1850 to 2014 described in Danabasoglu et al. (2020) as well as land use and land cover change described in Lawrence et al. (2019) following CMIP6 protocols outlined by Eyring et al. (2016).

Using a single ensemble member of CESM2‐esm‐hist, we focused our analysis on carbon fluxes that are simulated during the end of the historical period (1960–2014), which overlaps with atmospheric CO2 measurements. We quantified variability in carbon fluxes at interannual and seasonal time scales and correlated these fluxes with anomalies in climate drivers (temperature and moisture). The model simulates gross fluxes of GPP and Reco, with the difference between them representing net ecosystem production (NEP). Positive values for NEP represent net terrestrial ecosystem uptake of carbon. We focused on NEP instead of net ecosystem exchange or net biome production (NEE and NBP, respectively), which include fluxes from fire, land use, and land management, because the CESM2‐esm‐hist simulations have unrealistically large fire carbon fluxes from land degradation in the tropics at the end of the historical period.

Statistical Analyses

Interannual Variability

We summed monthly carbon fluxes (NEP, GPP, and Reco) that were simulated from vegetated terrestrial grid cells to calculate accumulated annual fluxes and weighted them by grid cell area and land fraction to calculate global values. We similarly calculated mean annual temperature (TBOT) and terrestrial water storage (TWS) that was simulated by the model. The CESM2‐esm results showed strong long‐term trends in relevant variables (Figure S1), so we subtracted linear trends and focused on the detrended anomalies from the climatological mean state for both monthly and annual data.

The IAV was calculated as the standard deviation of annual results simulated from 1960 to 2014. We compared the IAV in global detrended anomalies of simulated CO2 fluxes, simulated land CO2 fluxes, and NEP to those observed in the atmospheric CO2 growth rate reported by the Global Carbon Project (Friedlingstein et al., 2019). As most of the variability in global CO2 fluxes is driven by IAV in terrestrial fluxes (Figure S2), we focus on the IAV of NEP and its component fluxes.

Lagged correlations reveal teleconnections between modes of climate variability, like ENSO, that are associated with atmospheric CO2 anomalies and terrestrial CO2 fluxes (Braswell et al., 1997, see also Lindsay et al., 2014). Ecosystems, however, are responding to local climate conditions, where Bloom et al. (2020) also inferred significant ecosystem carbon cycle memory from lagged correlations with local climate anomalies. While both lagged responses are important features of terrestrial ecosystem response to climate variability, we focus here on documenting the immediate, lag‐zero response of fluxes simulated by CESM2‐esm. To characterize relationships in the data we calculated Pearson's correlation coefficients and regression statistics between detrended anomalies in carbon fluxes, moisture, and temperature. We also calculated the grid cell variance and standard deviation of annual, detrended carbon fluxes to look at covariation between simulated fluxes (as in Baldocchi et al., 2018).

Seasonal Modes of IAV

To decompose the annual cycle of carbon fluxes simulated in each terrestrial grid cell and identify modes of variability we used a singular value decomposition (SVD; Golub & Reinsch, 1971 as in Butterfield et al., 2020). The SVD decomposed the time series of detrended flux anomalies into singular vectors (SVi), the elements of which reflect the month (m) of the year (y; Figure S3). Vectors are ranked by the fraction of variance they explain in the time series. Each singular vector also receives a weight wi, one per year per singular vector, that quantifies the contribution from an individual singular vector to the observed IAV in any given year. Thus, the simulated IAV time series for a grid cell can be fully reconstructed as the weighted sum of singular vectors: Our SVD had 55 singular vectors (i, corresponding to the number of years in our analysis). Given recent publications that use and SVD to characterize seasonal variability in regional GPP fluxes (Butterfield et al., 2020; Byrne et al., 2020) we evaluated the model GPP using the same approach, but extend the analysis to also include grid‐cell variation in NEP simulated by CESM2‐esm, using the first two singular vectors to characterize the dominant modes of variability. We also calculated a redistribution metric, , as the sum of elements from a singular vector divided by the absolute values of the sum of elements from that vector (Butterfield et al., 2020) Thus, when  = 0, NEP was redistributed within the growing season without changing the annual flux. By contrast, values of  = 1 (or −1) indicate that every month had a positive (or negative) anomaly in GPP relative to the multi‐year mean

We identified the mode of variability corresponding to a seasonal amplification of NEP as the vector whose elements most strongly correlated with annual climatology of NEP. At annual time scales, the amplification vector corresponds to a net increase (or decrease) in terrestrial NEP, compared to the monthly detrended climatological mean. The other mode of variability corresponded to a seasonal redistribution of NEP, which typically has both positive and negative phases. In the northern hemisphere, this typically results from anomalously high plant productivity in the spring that is followed by lower‐than‐average plant productivity in the summer and fall (Figure S4 and results in Section 3.2; see also Butterfield et al., 2020). Thus, the redistribution vector does not necessarily result in a change in the net annual carbon flux, but it indicates changes in the seasonal timing of ecosystem production, relative to the monthly detrended climatological mean.

The values were used to confirm the appropriate identification of amplification and redistribution modes of variability in each grid cell (e.g., | |amplification > | |redistribution). To facilitate our analysis, we reversed the sign of singular vector elements, weights and values so that amplification vectors were positively correlated with the annual climatology of carbon fluxes and the redistribution vector started with a positive phase (Figure S4). For visualization, we calculated regional means for elements in the seasonal amplification and redistribution vectors across high, mid, and low latitude bands in both hemispheres (50°–80°, 20°–50°, and 0°–20°, respectively). Finally, to link modes of carbon cycle variability back to climate anomalies we calculated seasonal means for NEP, air temperature, and TWS anomalies. We looked at Pearson's correlation coefficients between these seasonal anomalies and the SVD weights generated for amplification and redistribution vectors.

Results

Interannual Variability

Detrended anomalies of terrestrial NEP that are simulated by CESM2‐esm have low variability, compared to the atmospheric growth rate of CO2 measured since 1959 and reported by the Global Carbon Project (Friedlingstein et al., 2019, Figures 1a and S2). The standard deviation of modeled NEP fluxes (0.47 Pg C y−1) is roughly half of the standard deviation in the observed atmospheric CO2 growth rate (0.95 Pg C y−1; Friedlingstein et al., 2019; Rödenbeck et al., 2018a). Notably, the standard deviation of annual detrended NEP fluxes in the tropics (23°S–23°N) simulated by CESM2‐esm is 0.24 Pg C y−1, roughly a factor of three lower than estimates from inversion models over land (90°S–25°N, 0.6–0.8 Pg C y−1; Rödenbeck et al., 2018a). Note, for convenience we inverted the sign of the atmospheric growth rate so that positive anomalies in Figure 1a show net land C uptake for both the model and observations. We also note that any temporal correlations between C flux anomalies in CESM2‐esm simulations and atmospheric observations here are unintended, because CESM2‐esm is experiencing a modeled atmosphere that does not necessarily match local, regional, or global conditions experienced during the historical record of atmospheric CO2 observations. The low IAV of net terrestrial C exchange with the atmosphere in CESM2‐esm could be driven by a number of factors that include low climate variability in the coupled model, low sensitivity of component fluxes (GPP and Reco) to climate variability, or high correlation between component fluxes at annual timescales that results in a dampening of IAV in the net carbon flux.

Figure 1

Detrended annual anomalies of global carbon fluxes, climate drivers, and their correlation. Upper panels show: (a) Atmospheric CO2 growth rate reported by the global carbon project (Friedlingstein et al., 2019) and net ecosystem production (NEP) simulated CESM2‐esm (green and black lines, respectively); (b) NEP, terrestrial water storage (TWS), and tropical air temperature anomalies simulated by CESM2‐esm (black, blue, and red lines respectively. Lower panels show correlations between simulated: (c) TWS, which is positively correlated to simulated NEP anomalies; and (d) Tropical air temperature, which is negatively correlated to simulated NEP anomalies. Note, for convenience we inverted the sign of the atmospheric growth rate in (a) so that positive anomalies in C fluxes show net land C uptake for both the model and observations.

To investigate the hypothesis that the low IAV of net terrestrial C exchange could be caused by low climate variability in in the coupled model, we plot in Figure 1b the detrended annual anomalies of NEP (as in Figure 1a) along with annual anomalies in TWS (over vegetated grid cells) and tropical air temperatures over land (23°S–23°N). The magnitude of TWS and tropical air temperature variability simulated in CESM2‐esm seems reasonable, compared to observations (σ = 0.82 Tt H2O and 0.10 K, respectively; Figure 1b; Cox et al., 2013; Humphrey et al., 2018). Consistent with the low carbon cycle variability simulated by CESM2‐esm (Figures 1a and 1b), however, the model shows weaker than observed carbon‐cycle sensitivity to global and regional climate anomalies (Figures 1c and 1d). The anomalies between NEP and TWS are statistically significant (r = 0.63, p < 0.001), but with a weaker than observed effect on carbon cycle variability (slope = 0.36 Pg C y−1 Tt H2O−1 vs. r = 0.85, slope = 1.3 Pg C y−1 Tt H2O−1 in Humphrey et al., 2018). We similarly find significant correlations between simulated anomalies of NEP and tropical temperature (r = −0.58, p < 0.001, slope = −2.6 Pg C y−1 K−1), but this relationship is also weaker than observed estimates [Cox et al. (2013), r = −0.65, slope = −5.1 Pg C y−1 K−1; Wang et al. (2013), r = −0.7, slope = −3.5 Pg C y−1 K−1] (see also Figure S5, Ballantyne et al., 2017; Rödenbeck et al., 2018a). Figure S5 shows that when forced with historical climate reanalysis land‐only simulations with CLM5 (Lawrence et al., 2019) capture the appropriate sign of terrestrial carbon cycle IAV observed in the atmospheric record, but not an appropriate magnitude of response, especially to the 1986–1987 or 1997–1998 El Niño events, resulting in low IAV of NEP fluxes (σ = 0.41 Pg C y−1). These land‐only simulations show similar IAV of TWS and tropical air temperatures as results from the coupled model (σ = 0.80 Tt H2O and 0.07 K, respectively), suggesting that the CESM2‐esm adequately represents climate variability over terrestrial ecosystems.

Observed and simulated variability in NEP is driven by variability in component fluxes, GPP and Reco. The global mean IAV of detrended GPP and Reco anomalies were 1.4 and 1.1 Pg C y−1, respectively. These values for GPP are in line with those reported by Chen et al. (2017) who looked at the IAV of GPP from ensemble of land models and upscaled flux tower estimates (1.4 and 1.5 Pg C y−1, respectively). Elsewhere, O'Sullivan et al. (2020) found little agreement in the magnitude of IAV of GPP simulated by different ensemble of land models (which included CLM4.5), upscaled flux tower observations, and remote sensing derived light‐use efficiency models, especially in highly productive tropical forests. This highlights the challenge in assessing global scale variability in plant productivity without direct observations. The challenges are even greater for ecosystem respiration, which has no remote sensing proxy and can only be indirectly inferred as the residual of other component fluxes, which are themselves highly uncertain (Basile et al., 2020; Chevallier, 2018; Liu et al., 2017).

Disaggregating the global fluxes simulated by CESM2‐esm into their gridcell contributions provides another perspective on component fluxes simulated by the model. The mean standard deviation of detrended GPP and Reco anomalies (102 and 83.8 g C m−2 y−1, respectively) in CESM2 were about 10% of the mean fluxes (950 and 910 g C m−2 y−1, respectively). Again, these estimates of GPP IAV agree with those reported from an ensemble of land models (96 g C m−2 y−1; O'Sullivan et al., 2020), but we note that the ensemble of prognostic models generally simulate lower IAV than FLUXCOM estimates and light‐use efficiency models (see also Jung et al., 2017). CESM2‐esm shows high IAV of detrended GPP fluxes in tropical savannah regions, but relatively low coefficient of variation in humid tropical forests (Figure S6). This feature is commonly seen in land models (O'Sullivan et al., 2020) and is inconsistent with observations from tropical forests (Bowman et al., 2017; Liu et al., 2017). Thus, we see the low IAV in simulated GPP fluxes as a major limitation in current land models, including CLM5, and note this as a priority area for future model developments.

Outside of the tropics, we compared the coefficients of variability (defined as the ratio of the interannual standard deviations to the seasonal amplitude of the multi‐year mean) in four ecoregions of North America as defined in Butterfield et al. (2020) from CESM2 and several remote sensing products (Table S1). The various remote sensing products show a factor of two difference in the range of variability in regional GPP estimates, and the magnitude of GPP variability simulated by CESM2 is comparable to these observationally derived estimates. This finding contradicts results reported by Wozniak et al. (2020), who found lower than observed variability in land‐only simulations conducted with CLM5 at flux tower sites.

Grid cell variance of NEP was strongly and positively correlated with the variance of component fluxes (Figure 2a). Although variance in NEP is slightly better explained by GPP variance (r = 0.79, slope = 0.77) than Reco variance (r = 0.74, slope = 0.72; Figure 2a), the anomalies of the component fluxes are highly correlated with each other (r = 0.94, slope = 0.80; Figure 2b). By contrast, observations from FLUXNET show lower correlations between GPP and Reco anomalies (r = 0.70, slope = 0.42; Baldocchi et al., 2018). We acknowledge that estimates of GPP and Reco from flux tower observations are a modeled product derived from NEE measurements. Moreover, these measurements are biased toward to relatively mesic temperature environments, creating a scale mismatch when used to evaluate results from a global model (Pastorello et al., 2017; Schimel, Pavlick et al., 2015). Regardless, the strong simulated correlation between simulated GPP and Reco likely accounts for some of the low IAV in land carbon uptake in CESM2‐esm, since years with large GPP fluxes are necessarily compensated by Reco fluxes that are of nearly the same magnitude.

Figure 2

Correlations of grid cell carbon fluxes simulated by CESM2‐esm from 1960 to 2015: (a) Variance in net ecosystem production (NEP) versus gross primary productivity (GPP) and Reco (shown in black and green, respectively); and (b) Correlation of anomalies in Reco and GPP. Pearson correlation coefficients and regression slopes from each relationship are provided in each panel, all correlations are significant (p < 0.001).

Numerous developments went into CLM5 (described in Lawrence et al., 2019), but we briefly explore the underlying causes for the high correlation between GPP and Reco fluxes in the model that are shown in Figure 2. The components of Reco include autotrophic and heterotrophic respiration. We found that annual gridcell Reco anomalies were most strongly correlated with annual anomalies in autotrophic respiration fluxes (Figure S7a). The components of autotrophic respiration include growth respiration, maintenance respiration, and a new plant respiration flux that is associated with the fixation and uptake of nitrogen model (FUN; R. A. Fisher et al., 2019; Shi et al., 2016). Of these, we find that anomalies associated with FUN are most strongly correlated with the autotrophic respiration fluxes simulated by CESM2‐esm (Figure S7b). From a mechanistic perspective, this high correlation makes sense, since for each additional increment in GPP, a corresponding increment in respiration is required to meet associated nitrogen demand associated with plant growth. This flux, however, occurs instantaneously in CLM5, and thus dampens potential IAV of net fluxes. In real ecosystems there is longer lag between photosynthesis and autotrophic respiration fluxes, notably those associated with plant nutrient uptake. Measuring and modeling these fluxes, especially related to non‐structural carbohydrates (Martínez‐Vilalta, et al., 2016) remains an outstanding challenge.

The IAV of detrended NEP anomalies simulated in CESM2‐esm is positively correlated with TWS anomalies in low and mid latitude regions (50°S–50°N), whereas high latitude ecosystems show a negative correlation between NEP and water storage anomalies (Figure 3a). These patterns are reversed for correlations between NEP and air temperature anomalies (Figure 3b). Over cold regions, and especially in boreal forests, the IAV of detrended NEP anomalies are positively correlated with air temperature anomalies. By contrast, over mid and low latitudes, especially in the Amazon, SE Asia, and N Australia, the IAV of detrended NEP anomalies are negatively correlated with air temperature anomalies. At annual time scales NEP anomalies show stronger correlations with GPP anomalies than they do with ecosystem respiration anomalies (Figures 3c and 3d). Figure S8 shows correlations of annual GPP, total water storage, and temperature anomalies.

Figure 3

Correlation coefficients between detrended annual anomalies that are simulated by CESM2‐esm from 1960 to 2014. Panels show the correlation between net ecosystem production (NEP) and (a) terrestrial water storage, (b) air temperature, (c) gross primary productivity (GPP), and (d) Reco. Only statistically significant correlations (p < 0.05, when |r| > 0.226 for 55 years of data) are shown.

Seasonal Variability

The first two vectors in the SVD explained 68% (area weighted mean) of the variance in NEP over the vegetated land surface. In general, the amplification vector explained the greatest fraction of variance in NEP, especially in arctic and boreal regions, as well as arid and temperate regions like the western North America (global area weighted mean = 39%; Figure 4a). The redistribution vector explained the largest fraction of variation over more mesic regions in the mid latitudes, especially temperate deciduous forest and agricultural regions in the eastern North America and Europe (global weighted mean = 29%; Figure 4b). Neither vector explained a large fraction of NEP variance over tropical forests, which generally showed low variability in detrended GPP anomalies (Figure S6). Figure S9 shows SVD results for GPP, which explained 75% for the total variance in GPP (45% and 29%, for amplification and redistribution vectors, respectively). The modes of variability simulated by CESM2‐esm with amplification and redistribution in western and eastern United States, respectively, are qualitatively similar to patterns observed in remote sensing estimates of GPP and flux inversions of NEE (Byrne et al., 2020).

Figure 4

Fraction of variance in detrended net ecosystem production (NEP) anomalies that was explained by (a) seasonal amplification or (b) seasonal redistribution vectors. Globally, these two vectors explained 39% and 29% of the variance in monthly NEP anomalies, respectively.

Figure 5 shows the mean annual climatology of the NEP, as well as the monthly values for amplification and redistribution vectors (gray, blue, and red lines respectively) for latitudinal bins. High and mid latitudes are characterized by a distinct annual cycle of NEP that is strongly correlated with the amplification vectors describing NEP variability (Figures 5a, 5b, and 5d). The amplification vector describes 48%, 40%, and 32% of the NEP variability in arctic, northern hemisphere temperate, and southern hemisphere temperate latitudinal zones, respectively. By contrast, the redistribution vectors in these regions explain 26%–28% of the NEP variability and are characterized by positive spring‐time anomalies that are followed by negative NEP anomalies in the summer and fall, relative to the monthly climatology. In the tropics, the seasonal cycle is more muted, but seasonal amplification and redistribution vectors describe roughly 35% and 28% of the variability in NEP, respectively (Figures 5c and 5e). Figure 6 shows analogous plots for the mean monthly vectors for the SVD of GPP anomalies, which illustrates modes of seasonal variability in GPP that are simulated by CESM2‐esm are similar to those reported from remote sensing observations (Butterfield et al., 2020; Byrne et al., 2020).

Figure 5

Zonal mean climatology of monthly net ecosystem production (NEP) and singular vectors associated with seasonal amplification and temporal redistribution of NEP (gray, blue, and red lines respectively) for the northern hemisphere and southern hemisphere (top and bottom rows, respectively). Panels show: (a) high latitude ecosystems, 50°–80°N; (b), (d) northern and southern temperate mid latitudes, 20°–50°N and 20°–50°S; and the (c), (e) tropics, 0°–20°N and 0°–20°S, respectively. The magnitude of the singular vectors is arbitrary (y‐axis). Mean fraction of variance explained and θ values, which indicate the net impact on the integrated seasonal signal of NEP for each singular vector, are also provided. Note x‐axis was shifted for southern hemisphere plots to show the climatology for austral summer.

Figure 6

Zonal mean climatology of monthly gross primary productivity (GPP) and singular vectors associated with seasonal amplification and temporal redistribution of GPP (gray, blue, and red lines respectively) for the northern hemisphere and southern hemisphere (top and bottom rows, respectively). Panels organized as in Figure 5. Note x‐axis was shifted for southern hemisphere plots to show the climatology for austral summer.

The θ values calculated in the SVD show the net impact on the integrated seasonal signal of NEP. The mean θ values associated with the amplification vector are globally positive (Figure 5). By contrast, the θ values associated with a seasonal redistribution of NEP are negative, or close to zero. This indicates that a negative, or little change in integrated seasonal signal of net carbon fluxes from seasonal redistribution of NEP. The redistribution vectors of GPP anomalies largely have θ values close to zero (Figure 6). Notably, redistribution of GPP in temperate and boreal regions leads to no net change in the annual uptake of CO2 (Figure 6). Net fluxes, however, in high latitude ecosystems show an annually integrated loss of C from land to the atmosphere with seasonal redistribution of fluxes (θ < 0, Figure 5a). These dynamics are driven by concurrent, higher than average (positive) spring anomalies in heterotrophic respiration fluxes (Figure S10) that cancel positive anomalies in GPP, resulting in negligible net C uptake in the spring (Figure 5a).

Although the seasonal redistribution of carbon variability is a major source of global carbon cycle variability in the model, it would not be evident in more aggregated metrics of variability that only look at annual times‐scales (e.g., Figures 1 and S5).

To evaluate the environmental drivers of carbon cycle variability we conducted linear regressions on the weights of the amplification (Figure 7) and redistribution (Figure 8) vectors from SVD analysis (n = 55, for each year of the simulation) with seasonal anomalies of carbon fluxes, TWS, and air temperature. We focus our analysis on GPP to isolate between the environmental sensitivities of component fluxes of NEP. The SVD weights for the amplification vector of GPP were strongly and positively correlated with the GPP anomalies during the peak of the growing season (Figure 7, left column). This is expected, since the amplification vector was identified by its correlations with the climatology of GPP. As such, correlation coefficients are highest in the summer months (JJA and DJF for northern and southern hemispheres, respectively). Strong correlations are evident at other times (e.g., negative correlations between weights and GPP anomalies across high latitudes in DJF), but the magnitude of these anomalies is small relative to the annual cycle (see also Figure 6a). Although TWS and air temperature covary (Figure S8c), high‐latitude ecosystems generally show SVD amplification weights that are more positively correlated with air temperature anomalies and negatively correlated with water storage anomalies in JJA (Figure 7, right and middle columns, Table S2). By contrast, across mid and low latitudes, SVD amplification weights for GPP are more strongly and positively correlated with wetter‐than‐average conditions and negatively correlated with air temperature anomalies.

Figure 7

Pearson correlation coefficients between singular value decomposition (SVD) weights from the amplification vector with seasonal anomalies of gross primary productivity (GPP), terrestrial water storage (TWS), and air temperature (TBOT) simulated by CESM2‐esm from 1960 to 2014. Only statistically significant (p < 0.05) correlations are shown (|r| > 0.26, two‐tailed test, n = 55).

Figure 8

Pearson correlation coefficients between singular value decomposition (SVD) weights from the redistribution vector with seasonal anomalies of gross primary productivity (GPP), terrestrial water storage (TWS), and air temperature (TBOT) simulated by CESM2‐esm from 1960 to 2014. Only statistically significant (p < 0.05) correlations are shown (|r| > 0.26, two‐tailed test, n = 55).

The SVD weights for the redistribution vector were strongly and positively correlated with the GPP anomalies during the spring, with correlation coefficients that are highest in the MAM and SON (for northern and southern hemispheres, respectively; Figure 8, left column). In the northern hemisphere, the positive phase of the redistribution vector is also more strongly correlated with warmer spring‐time air temperatures (Figure 8, right column) than TWS. Subsequent GPP anomalies in the summer and fall, however, show negative correlations with SVD weights (see also Figure 6). The summer (JJA) is still characterized by warmer, but also drier‐than‐average conditions (Figure 8, middle column; Table S3) that also have negative GPP anomalies. By fall (SON), these negative GPP anomalies are only associated with drier‐than‐average conditions. Thus, the positive phase of the redistribution vector associated with warmer‐than‐average springs is (nearly) canceled by the subsequent soil moisture stress and lower‐than‐average productivity in the summer and fall, resulting in minimal change in annual productivity.

Discussion

Our results show that the magnitude of global carbon cycle IAV simulated by CESM2‐esm is low, relative to measurements of IAV in the atmospheric CO2 growth rate (Figure 1). The CESM2 has a good representation of precipitation and temperature anomaly patterns associated with tropical Pacific sea surface temperatures (Danabasoglu et al., 2020; Meehl et al., 2020). Accordingly, the magnitude of TWS and air temperature variability agrees reasonably well with observations (Figure 1; Cox et al., 2013; Humphrey et al., 2018). This suggests that the model adequately represents global‐scale climate variability, but that this climate variability does not generate enough terrestrial carbon cycle variability in the model.

Because the timing of moisture or temperature variability in the coupled model are not expected to agree with the timing of these variations in the historical, observational record, we instead analyze relationships between carbon fluxes and climate anomalies. The model shows lower than observed sensitivity of net C fluxes to climate variability in moisture and temperature anomalies (Figures 1c and 1d; Cox et al., 2013; Humphrey et al., 2018; Wang et al., 2013). We appreciate that these studies all use slightly different time periods for their calculations, but we do not expect this to significantly alter the fundamental sensitivities of land carbon fluxes to climate variability. Including terrestrial CO2 fluxes associated with land use change and fire does not meaningfully change the magnitude of IAV of net carbon fluxes that are simulated by the model (Figure S2; excluding the year 2010). The finding that simulated C fluxes in the coupled model have low IAV is generally consistent with land‐only CLM5 simulations, which also shows weaker than observed carbon cycle variability when driven by reanalysis climate forcing data that reflects the historical drivers of the observed carbon cycle variations (Figure S5), even when including C fluxes associated with land use change and fires (Lawrence et al., 2019). Notably, these offline simulations capture the sign of terrestrial carbon cycle IAV observed in the atmospheric record, the model does not simulate the appropriate magnitude of response to the 1986–1987 or 1997–1998 El Niño events.

The IAV in net carbon fluxes results from variability in component fluxes (GPP and Reco) and their interaction (Baldocchi et al., 2018; Lasslop et al., 2010). Simulations from CESM2‐esm show strong correlations between NEP and both of its component fluxes (Figure 2a). Our results suggest that CESM2‐esm shows notably low IAV of GPP fluxes in tropical forests (Figure S6). This feature of the terrestrial C cycle simulated by CESM2‐esm (and CLM5) is not consistent with observations that suggest tropical forests respond strongly to variability in the climate system (Bastos et al., 2018; Cavaleri et al., 2017; Liu et al., 2017; Rödenbeck et al., 2018a). Ensembles of land models typically do a good job of matching the IAV of the atmospheric growth rate (Ahlström et al., 2015; Friedlingstein et al., 2019; Poulter et al., 2014), but O'Sullivan et al. (2020) recently found that land models, including CLM4.5, generally simulate low IAV of GPP fluxes in tropical forests. This suggests broader deficiencies in land models' ability to capture appropriate mechanisms for C cycle IAV, which are also evident in our CESM2‐esm simulations. Future work, therefore, should focus on understanding and improving the simulation of gross and net carbon cycle sensitivities to climate variability, especially in the humid tropics. Notably uncertainties in tropical forest IAV are related to how models like CLM simulate phenology, as well as temperature and moisture stress on photosynthesis and stomatal conductance. Indeed, O'Sullivan et al. (2020) emphasize this broader challenge in understanding the magnitude and drivers of tropical forest variability in models and observations.

Low variability in plant productivity simulated by CESM2‐esm may extend beyond the tropics. For example, Wozniak et al. (2020) found that maximum rates of GPP simulated by CLM in temperate deciduous forests were much lower than observations at a number of AmeriFlux sites. Indeed, measurements from flux towers suggest that brief periods of large photosynthetic uptake appear to be an important component of the IAV in net carbon exchange, especially in arid ecosystems (Fu et al., 2019; Kannenberg et al., 2020; Zscheischler et al., 2016). The failure of CLM5 to simulate high rates of GPP during favorable conditions suggests that the model needs parametric or structural changes in its representation of leaf‐level photosynthesis, stomatal conductance, or canopy scaling to represent photosynthetic variability. Regional analyses comparing the IAV of GPP in CESM2‐esm simulations to remote sensing estimates suggest that the variation in GPP fluxes may be appropriate for four regions of North America (Table S1; Butterfield et al., 2020). Additional work is, however, needed to evaluate the feasibility of quantifying IAV of carbon cycle components from remote sensing products (Bowman et al., 2017; Liu et al., 2017; Schimel, Pavlick et al., 2015) and comparing the result to Earth system models, especially in the tropics.

Observations suggest that site‐level variance in net carbon fluxes is more tightly correlated with GPP than Reco (Baldocchi et al., 2018), which is also true in our simulations (Figures 3c and 3d). Broadly, we note that models like CLM tends to overestimate respiration fluxes because they do not simulate other ecosystem C losses pathways that may have longer residence times (Ciais et al., 2020; Randerson et al., 2002; Trumbore, 2006). The anomalies of GPP and Reco that are simulated by CESM2‐esm are more strongly correlated than observations suggest (Figure 2b and Baldocchi et al., 2018). This high covariation between GPP and Reco offsets variance in either of the component fluxes and dampens the IAV of NEP in the model. These results also suggest that the current structure and parameterization of CESM2‐esm, which dictates the high covariance of simulated GPP and Reco, should be revised (although evaluating Reco fluxes is challenging without spatiotemporally explicit estimates for model benchmarking). We suspect that several changes made to the land model for CESM2 have exacerbated this problem.

First, CLM5 reduced the magnitude of growth respiration fluxes (Atkin et al., 2017) and reduced the total C fluxes associated with growth and maintenance respiration fluxes, relative to previous versions of the model (Lawrence et al., 2019). Second, the incorporation of the Fixation and Uptake of Nitrogen model (FUN) into CLM5 makes plants pay the carbon costs of nitrogen uptake (J. B. Fisher et al., 2010; R. A. Fisher et al., 2019; Shi et al., 2016). As currently applied in CLM5, the anomalies in FUN carbon costs are highly correlated with autotrophic respiration and, by extension, NEP anomalies (Figure S7). This follows the logic currently built into the model, whereby plants instantaneously allocate new photosynthate to build biomass or meet respiration demands and remaining carbon is used for nitrogen acquisition by FUN. In real ecosystems this expenditure of carbon to acquire nitrogen would occur below ground, with at least some of the carbon interacting with the soil matrix (e.g., through root exudates). Moreover, this transfer of fresh photosynthate belowground for nitrogen uptake does not occur instantaneously, as it does in the model. Future developments should address these deficiencies in CLM5. Finally, the parameterization for soil organic matter turnover uses a higher minimum water potential, which increases the sensitivity of heterotrophic respiration fluxes to liquid soil water availability (Carvalhais et al., 2014; Koven et al., 2017; Lawrence et al., 2019; Z. Liu et al., 2018). Independently, these changes to growth respiration, FUN, and heterotrophic respiration fluxes seem justified in their aim to represent terrestrial ecosystems more realistically, but together they likely served to reduce the IAV of net carbon fluxes that are simulated by CESM2‐esm.

The correlations of NEP and GPP anomalies with climate drivers in CESM2‐esm show strong latitudinal patterns (Figures 3, S8). Notably, correlations between carbon flux and TWS anomalies are particularly strong in many arid, semi‐arid and savannah regions (Figures 3a and S8a), a finding that is consistent with work emphasizing moisture and precipitation controls over carbon cycle variability in arid regions (Ahlström et al., 2015; Humphrey et al., 2018; Poulter et al., 2014). Concurrently, correlations between carbon flux and air temperature anomalies are stronger in arctic, boreal and temperate deciduous forests (Figures 3b and S8b), which again is consistent with observations (discussed in Section 4.2; Hu et al., 2019; Rödenbeck et al., 2018b). We recognize that inferring the relative importance of climate controls over land‐atmosphere carbon exchange remain actively discussed in the literature (Cox et al., 2013; Humphrey et al., 2018; Jung et al., 2017; Piao et al., 2020; Poulter et al., 2014). Given the spatial and temporal heterogeneity of climate anomalies and timescales of ecosystem responses (Rödenbeck et al., 2018b; X. Zhang et al., 2013) we investigated the seasonal modes of carbon cycle variability that are simulated by CESM2‐esm and their environmental covariates.

The timing of climate variations with respect to the climatological annual cycle plays an important role in the resulting seasonal and IAV of terrestrial carbon fluxes (Buermann et al., 2018). Satellite and flux tower observations in North America suggest that carbon cycle variability can be decomposed into modes of variability that are characterized by the amplification and redistribution of seasonal fluxes (Butterfield et al., 2020; Byrne et al., 2020). Likewise, the SVD analysis identified similar modes of variability in NEP and GPP simulated in CESM2‐esm (Figures 4 and S9). Amplification vectors are dominant in high latitude and arid ecosystems, while redistribution vectors that are dominant in temperate forests, boreal forests, and agricultural regions. Qualitatively, these patterns align with findings from (Butterfield et al., 2020) who found robust patterns in seasonal GPP variability from several satellite datasets that are correlated with regional anomalies of temperature and soil moisture availability. The similarity in spatial patterns of NEP and GPP variability form our SVD analyses illustrate the strong dependence of modeled NEP variability on GPP, but we recognize that seasonal changes in Reco that are also important drivers of carbon cycle variability (Basile et al., 2020; Commane et al., 2017; J. Liu et al., 2017; Z. Liu et al., 2018; Welp et al., 2007).

In grid cells where the seasonal amplification of carbon cycle fluxes characterizes most of the flux variability the second vector from the SVD corresponds to a seasonal redistribution of the fluxes (Figures 4, 5, 6). For example, in high latitude ecosystems the amplification vector describes nearly half of the variability in simulated NEP and is associated with a net increase (or decrease) in annual carbon fluxes (Figure 5a). The weights associated with the amplification vectors are most strongly correlated with summertime GPP anomalies (Figure 7, top row; Table S2), which is not surprising since we identified the amplification vector from the SVD by its correlation with the mean climatology of monthly fluxes simulated in each grid cell (Section 2.2; Figure S4). The weights from the amplification vectors in high latitude ecosystems also show a strong, positive correlation with summertime air temperature anomalies and a weaker, but still significant, negative correlation with TWS anomalies. Thus, with warmer (and drier) summertime conditions, CESM2‐esm simulated positive GPP anomalies in arctic and boreal ecosystems. Conversely, with cooler (and wetter) summertime conditions CESM2 simulated negative GPP anomalies in these regions. Our analysis cannot diagnose the proximal driver of the GPP anomalies, but given their higher correlation coefficient we assume that summertime temperature anomalies are driving the carbon cycle response, with declines in soil moisture subsequently resulting from higher evapotranspiration fluxes in warmer, more productive years.

The amplification vector also describes a high fraction of the NEP and GPP variability in lower latitudes (Figures 5 and 6). The weights from the GPP amplification vectors in mid and low latitudes shows a strong positive correlation with regional peak growing season GPP anomalies in their respective hemispheres (Figure 7). In contrast to northern high latitudes, weights from GPP amplification vectors in these regions generally show stronger correlations with TWS anomalies than they do for air temperature (Figure 7, Table S2). Thus, in mid and lower latitudes wetter (and cooler) anomalies during the growing season maxima are associated with positive GPP anomalies, whereas drier (and warmer) anomalies are associated with negative GPP anomalies. While seasonal amplification vectors do explain most of the global carbon cycle variability, some regions are better characterized by a seasonal redistribution of carbon fluxes that do not necessarily change the annual carbon flux, just its timing.

The seasonal redistribution of NEP and GPP explains roughly a quarter of global variability in these fluxes (Figures 4, 5, 6 and S9). This seasonal redistribution of net and gross carbon fluxes is the dominant mode of variability simulated by CESM2‐esm in several regions, including temperate forests and agricultural regions in Eastern North America and Europe. The spatial cohesiveness of this pattern is most notable in the northern hemisphere, where the redistribution vector has positive spring GPP anomalies that are associated with warmer‐than‐average spring temperatures (MAM; Figure 8). These positive springtime GPP anomalies are subsequently compensated by lower‐than‐average photosynthesis rates in the summer and fall that tend to be associated with increasing soil moisture stress (JJA and SON; Figure 8, Table S3). Drier summer and fall conditions could result from higher evapotranspiration in the spring, or also from increased early runoff due to earlier snowmelt during warm springs (Buermann et al., 2013). The ability of CESM2 to replicate the redistribution of NEP and GPP fluxes (Figures 4 and 8; Byrne et al., 2020) suggests that the representation of plant phenology and water stress in CESM2 are likely responding in physically and ecologically realistic ways to simulated climate variability. With seasonal redistribution, potential increases in plant productivity from an early green‐up that were facilitated by warmer spring temperatures are negated by soil moisture stress later in the growing season, leading to negligible changes in the annual rates of gross and net carbon uptake (mean values close to zero for the redistribution vector; Figures 5 and 6).

The regions where a seasonal redistribution vector dominates carbon cycle variability in CESM2‐esm (Figures 4 and S9) are also regions where native vegetation in the model use a stress deciduous phenology scheme, or they are under agricultural management, which is explicitly represented by CLM5 in CESM2 (Lawrence et al., 2019; Lombardozzi et al., 2020). Both phenology parameterizations use a growing degree day approach to simulate leaf emergence (or planting date and leaf emergence for the CLM5 crop model), so the strong correlation between air temperature anomalies and SVD weights are expected (Figure 8, Table S3). We were more surprised, however, by the negative GPP and NEP anomalies that emerge later in the growing season. These seem to be driven by drier‐ and warmer‐than‐average conditions that are consistent with satellite observations of vegetation greenness (Buermann et al., 2013, 2018). Notably, redistribution vectors diagnosed from solar induced fluorescence (SIF) derived GPP in North America are tightly linked with spring (and summer) temperature anomalies (Butterfield et al., 2020; Byrne et al., 2020; Table S1). The larger influence of seasonal redistribution in temperate latitudes that is simulated in CESM2‐esm is also consistent with observations from forests reported in Butterfield et al. (2020), but the overall importance of seasonal redistribution versus amplification on carbon cycle variability remains uncertain. Indeed, considering the relative importance of these modes of variability may be important in trying to infer appropriate sensitivities and interactions between seasonal to IAV in climate, phenology, and ecosystem carbon fluxes from both models and observations.

Conclusion

The IAV of terrestrial net carbon exchange with the atmosphere in CESM2‐esm is low. Accordingly, the model also simulates a weaker‐than‐observed sensitivity of net carbon exchange to global climate anomalies. This low variability of net carbon fluxes likely results from low IAV in photosynthesis anomalies that are simulated by the model, especially in humid tropical forests, and a high covariation GPP and Reco. Our findings emphasize the need to evaluate and revise the parameterization of carbon cycle sensitivity to climate variability, especially in the tropics, that are related to phenology, photosynthetic and stomatal conductance responses to environmental stress, the canopy scaling of these leaf level processes and potential legacy effects from extreme events. The variability in GPP that is simulated by CESM2‐esm generally shows a latitudinal gradient in climate sensitivities whereby positive NEP and GPP anomalies are driven by warmer and drier conditions in high latitude ecosystems but wetter and cooler conditions across mid and low latitudes.

Our analysis decomposes IAV in terrestrial carbon cycle anomalies into modes of variability that are characterized by seasonal amplification and redistribution of carbon fluxes. Together these modes of variability explain two‐thirds of the variability in NEP anomalies and three‐quarters of the global variability in GPP. Both CESM2 and observations show that wetter and cooler springs and summers lead to an amplification signal in GPP, especially over the high latitude ecosystems and temperate and arid regions like western North America. By contrast, the temporal redistribution of GPP anomalies is more strongly associated with variability in springtime temperatures and positive C cycle anomalies that are followed by subsequent soil moisture stress and negative carbon cycle anomalies. This seasonal redistribution of fluxes is the dominant mode of carbon cycle variability in the more mesic temperate ecosystems like the eastern North America and Europe. The seasonal redistribution of carbon cycle variability is notable because although it is not apparent in more aggregated (annual) measurements of IAV, it does seem widespread in both the model and observations. Thus, while the total magnitude of net and gross terrestrial carbon flux variability simulated by CESM2 may be too low, the simulated interannual and seasonal variability does qualitatively capture patterns of carbon sensitivities to climate variability. More broadly, we contend this kind of analysis is useful in diagnosing strengths and weaknesses in biogeochemical models in comparison to observational data.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Supporting Information S1

Click here for additional data file.