Compensatory effects conceal large uncertainties in the modelled processes behind the relationship between the El Niño–Southern Oscillation (ENSO) and CO₂

A large fraction of the interannual variation in the global carbon cycle can be explained and predicted by the impact of the El Niño–Southern Oscillation (ENSO) on net biome production (NBP). It is therefore crucial that the relationship between ENSO and NBP is correctly represented in Earth system models (ESMs). In this work, we look beyond the top-down ENSO–CO₂ relationship by describing the characteristic ENSO–NBP pathways in 22 Coupled Model Intercomparison Project Phase 6 (CMIP6) ESMs. These pathways result from the configuration of three interacting processes that contribute to the overall ENSO–CO₂ relationship: ENSO strength, ENSO-induced climate anomalies, and the sensitivity of NBP to climate. The analysed ESMs agree on the direction of the sensitivity of global NBP to ENSO but exhibit very high uncertainty with regard to its magnitude, with a global NBP anomaly of −0.15 to −2.13 Pg C yr⁻¹ per standardised El Niño event. The largest source of uncertainty lies in the differences in the sensitivity of NBP to climate. This uncertainty among the ESMs increases even further when only the differences in NBP sensitivity to climate are considered. This is because differences in the climate sensitivity of NBP are partially compensated for by ENSO strength. A similar phenomenon occurs regarding the distribution of ENSO-induced climate anomalies. We show that even models that agree on global NBP anomalies exhibit strong disagreement with regard to the contributions of different regions to the global anomaly. This analysis shows that while ESMs can have a comparable ENSO-induced CO₂ anomaly, the carbon fluxes contributing to this anomaly originate from different regions and are caused by different drivers. These alternative ENSO–NBP pathways can lead to a false confidence in the reproduction of CO₂ by assimilating the ocean and the dismissal of predictive performance offered through ENSO. We suggest improving the underlying processes by using large-scale carbon flux data for model tuning in order to capture the ENSO-induced NBP anomaly patterns. The increasing availability of carbon flux data from atmospheric inversions and remote sensing products makes this a tangible goal that could lead to a better representation of the processes driving interannual variability in the global carbon cycle.

1 Introduction

The relationship between the El Niño–Southern Oscillation (ENSO) and atmospheric CO₂ observations at Mauna Loa was first reported by Bacastow (1976). Altered atmospheric circulation patterns during El Niño events cause warm and dry conditions across the tropics, leading to a reduction in net biome production (NBP) due to reduced net primary productivity (NPP) and increased or decreased heterotrophic respiration (Rh) (Qian et al., 2008; Bastos et al., 2018). ENSO-induced climate anomalies have a significant impact on the gross primary productivity (GPP) of 32 % of the vegetated land area and can explain up to 26 % of the interannual variation in global GPP (Zhang et al., 2019). Some El Niño events can be severe enough to turn the Amazon Basin, a carbon sink of global importance, into a net carbon source (Tian et al., 1998).

However, ENSO not only explains a large fraction of the NBP variability but also is the main source of seasonal-to-decadal predictability in the Earth system (Manzanas et al., 2014; Zeng et al., 2008; Spring and Ilyina, 2020; Li et al., 2022). Tropical carbon flux anomalies lag behind ENSO by 3 to 6 months (Zhu et al., 2017; Wang et al., 2016), meaning that even without further knowledge of the evolution of ENSO, near-term carbon flux variability can be anticipated based on present ENSO conditions. On top of this lag effect, Earth system model (ESM) simulations starting in winter can predict ENSO conditions for up to 1 year (Barnston et al., 2019).

Further predictability is added to the system by the land surface, which prolongs the ENSO-induced climate anomalies. The larger the anomaly, the longer it takes for soil moisture and conditions to return to normal, and ENSO years are often among the most extreme years of variability (Holmgren et al., 2001). Even longer predictability mechanisms might be triggered through vegetation dynamics (Holmgren et al., 2001). This can happen in dry years through the lasting impacts of defoliation and tree mortality (Wigneron et al., 2020; Santos et al., 2018) or through wildfires, which require decades of recovery (Silva et al., 2018). Wet events, on the other hand, can provide long-term predictability as these events drive plant recruitment in semi-arid ecosystems (Holmgren et al., 2001). Extreme events play a crucial role in the vegetation dynamics of these ecosystems, where the establishment of trees and shrubs needs sustained wet conditions (Chang-Yang et al., 2016; López et al., 2006).

Because ENSO plays such a large role in the variability and predictability of NBP, the correct representation of related processes in ESMs is especially important. Three key processes that explain the relationship between ENSO and NBP can be arranged hierarchically. At the highest level of this hierarchy is the strength of the ENSO events. Despite considerable advancements in our understanding of ENSO dynamics, there remains a wide range of simulated ENSO amplitudes in ESMs (Beobide-Arsuaga et al., 2021). The amplitude of ENSO can be measured as the standard deviation (SD) of sea surface temperature anomalies (SSTAs) in the Niño3.4 region (5° N–5° S, 170–120° W) and ranges from 0.4 to 1.4 °C in models from the sixth phase of the Coupled Model Intercomparison Project (CMIP6) (Brown et al., 2020; Beobide-Arsuaga et al., 2021; Cai et al., 2022).

The second process that explains the impact of ENSO on NBP involves ENSO-induced climate anomalies. Most ENSO teleconnections are caused by a reorganisation of tropical convection patterns (Perry et al., 2020). These patterns create temperature and precipitation anomalies in northern South America, Southeast Asia, and northern Australia. However, ENSO-induced changes in the upper atmosphere can create Rossby waves that propagate polewards and lead to climate anomalies in the midlatitudes. ENSO teleconnection strengths show high uncertainty among CMIP5 models, with an average correlation with observed teleconnection patterns of 0.7 for temperature and 0.46 for precipitation (Perry et al., 2020). Although the representation of the relationship between ENSO and tropical precipitation improved from CMIP5 to CMIP6, there are still considerable deviations from the observed relationship (Yang and Huang, 2022). The impact of ENSO on eastern Asian summer rainfall, for example, can only be captured by 11 out of 20 CMIP6 models (Fu et al., 2021).

The last stage within the hierarchy of processes that shape the ENSO–NBP relationship involves differences in biogeochemistry – specifically, the sensitivity of NBP to local climatic anomalies caused by ENSO. Due to limitations in carbon flux observations and the covariability of climatic conditions, the contributions of temperature and moisture to driving the carbon cycle remain a debated topic in the literature (Piao et al., 2020). This uncertainty can be demonstrated by the sensitivity of the atmospheric CO₂ growth rate to tropical temperatures from reanalysis data, where determined sensitivities differ by a factor of 2 (Piao et al., 2020). Regional differences in the sensitivity of carbon fluxes to climate depend on the ecosystem type and climate. Semi-arid ecosystems and tropical forests, which cover most of the land area affected by ENSO, show the highest sensitivity to climate variability (Bastos et al., 2013; Poulter et al., 2014; Ahlström et al., 2015; O'Sullivan et al., 2020). However, these biomes are also where carbon flux sensitivities have the highest inter-model spread and bias relative to observations (O'Sullivan et al., 2020; Koirala et al., 2022).

While the land carbon sink variability in ESMs can be constrained using the sensitivity of the entire tropical carbon sink to ENSO-driven climate anomalies (Cox et al., 2013; Zechlau et al., 2022), we aim to test whether a spatially explicit constraint is feasible. The aim of this study is to look beyond the relationship between ENSO and the global atmospheric CO₂ growth rate in ESMs and reveal the sources of uncertainty in the ENSO–NBP relationship. We quantify the specific ENSO–NBP pathways that describe the locations and drivers of the ENSO-induced NBP anomalies. These pathways are characterised by three main processes that shape the ENSO–NBP relationship: ENSO strength, ENSO-induced climate anomalies, and biogeochemistry. We quantify how much the uncertainties in these three processes contribute to the uncertainty in the ENSO–NBP relationship among the ESMs and compare the ESMs with observations.

Ziehn et al. (2020)Shi et al. (2020)Wu et al. (2019)Swart et al. (2019)Danabasoglu et al. (2020)Lovato et al. (2022)Cherchi et al. (2019)Séférian et al. (2019)Golaz et al. (2019)Döscher et al. (2022)Dunne et al. (2020)Orbe et al. (2020)Volodin et al. (2018)Boucher et al. (2020)Watanabe et al. (2021)Hajima et al. (2020)Mauritsen et al. (2019)Yukimoto et al. (2019)Bethke et al. (2021)Seland et al. (2020)Seland et al. (2020)Sellar et al. (2019)

Table 1 The Earth system models used in this study, as well as the experiment runs, the simulation lengths, and whether or not fire emissions are provided.

2 Methods

2.1 Data

We measure the interactions between ENSO, temperature, precipitation, and NBP across 22 CMIP6 ESMs and compare them with observation-based data sources. In order to have a large sample size of ENSO conditions, we use unforced pre-industrial control (piControl) simulations. The analysed ESMs, their variant labels, and the simulation lengths are listed in Table 1. We decompose the NBP anomalies into the following components:

$$\begin{array}{}\text{(1)}& \text{NBP}\approx \text{NPP}-\text{Rh}-\text{fire}.\end{array}$$

Because some ESMs simulate types of disturbances other than fires, Eq. (1) does not lead to an exact reproduction of NBP but provides an approximate value for the contribution of the large carbon fluxes to the land–air CO₂ exchange. Fire emissions are only available for 13 of the analysed ESMs. Although CMCC-CM2-SR5 simulates fire emissions, the data are not available online. Instead, we calculate fire emissions for CMCC-CM2-SR5 by inverting the mass balance in Eq. (1).

The carbon fluxes of the ESMs are compared with two observation-based data sets, referred to as observations here. These are land–atmosphere carbon fluxes from the Copernicus Atmosphere Monitoring Service (CAMS) (Chevallier et al., 2005), an atmospheric inversion product, and net ecosystem exchange (NEE) from upscaled flux tower measurements (FLUXCOM version: ANN.CRUNCEPv6; Jung et al., 2019). Although NEE includes land–atmosphere carbon fluxes that are not part of NBP (Ciais et al., 2022), such as geological CO₂ release and lateral fluxes, these fluxes only play a minor role in the global carbon cycle and are unlikely to significantly affect the sensitivity to ENSO (Canadell et al., 2021). We also compare the carbon fluxes to NBP from the TRENDY v6 ensemble of land surface models (Sitch et al., 2015). This data product contains the results of land surface models that are also components of the ESMs used in this work. Due to this dependence, we do not treat the TRENDY results as observations. To account for uncertainty in meteorological data, we use different reanalysis products. We use sea surface temperature (SST) reanalysis data (Hadley Centre Sea Ice and Sea Surface Temperature (HadISST); Rayner et al., 2003); temperatures from the ERA-Interim reanalysis (1979–2018; Dee et al., 2011); version 4 of the Climatic Research Unit gridded Time Series (CRU TS v4 (1959–2019); Harris et al., 2020); the Japanese 55-year Reanalysis (JRA-55 (1958–2019); Kobayashi et al., 2015); version 2 of the Modern-Era Retrospective analysis for Research and Applications (MERRA-2 (1980–2021); Gelaro et al., 2017); the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) 40-Year Reanalysis Project (1948–2021; Kalnay et al., 1996); bias-adjusted ERA5 reanalysis data (WFDE5 (1979–2018); Cucchi et al., 2020); and precipitation data from the ERA-Interim reanalysis, JRA-55, MERRA2, the NCEP–NCAR 40-Year Reanalysis Project, and the Rainfall Estimates on a Gridded Network (REGEN) data set (1950–2016; Contractor et al., 2020).

Figure 1 The 14 Intergovernmental Panel on Climate Change (IPCC) climate reference regions with the largest ENSO-induced NBP anomalies: SE Asia (SEA), northern South America (NSA), NE South America (NES), eastern southern Africa (ESAF), the South American monsoon (SAM) region, northern Australia (NAU), NW South America (NWS), central Africa (CAF), South Asia (SAS), central Australia (CAU), eastern Australia (EAU), southern Central America (SCA), western southern Africa (WSAF), and SE South America (SES).

2.2 Processing and analysis

We calculate the annual anomalies in Niño3.4 SST, temperature, precipitation, and carbon fluxes for the analysis. In this study, we base annual averages on a time window from July to June of the following year, rather than on a window from January to December. We use this definition of a year to better capture distinct ENSO events. ENSO SSTAs usually peak during boreal winter, and warm El Niño events are often followed by a cold La Niña event (An and Kim, 2017). Years starting in January are, therefore, not centred around event peaks and could thus contain the tail of an El Niño event and the beginning of a La Niña event. Annual anomalies in the data are obtained by subtracting the climatology from the ESM data and removing the linear trend from the observational products. ENSO strength is calculated as the SD of annual SSTAs for the Niño3.4 region. NBP and climate anomalies are aggregated into regional averages using the boundaries of the updated Intergovernmental Panel on Climate Change (IPCC) reference regions (Iturbide et al., 2020). For this analysis, we only focus on the regions with the strongest ENSO–NBP relationships. These regions are identified by averaging the ENSO-induced NBP anomalies across all data sources and selecting the 14 regions with the highest absolute ENSO-induced NBP anomalies (Fig. 1). We calculate the ENSO-induced climate and carbon anomalies using the coefficients of a linear regression model with a zero intercept:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{\Delta }& {\displaystyle }{\text{Temp}}_{ijp}={\mathit{\beta }}_{{\mathrm{ET}}_{ij}}\times \mathrm{\Delta }{\text{ENSO}}_{jp},\text{(3)}& {\displaystyle }\mathrm{\Delta }& {\displaystyle }{\text{Precip}}_{ijp}={\mathit{\beta }}_{{\mathrm{EP}}_{ij}}\times \mathrm{\Delta }{\text{ENSO}}_{jp},\text{(4)}& {\displaystyle }\mathrm{\Delta }& {\displaystyle }{\text{NBP}}_{ijp}={\mathit{\beta }}_{{\mathrm{EN}}_{ij}}\times \mathrm{\Delta }{\text{ENSO}}_{jp},\end{array}$$

where ΔX_ijp represents the annual anomaly, ΔENSO_jp represents the annual mean Niño3.4 SSTA, the regression coefficient (${\mathit{\beta }}_{{\mathrm{EX}}_{ij}}$) corresponds to the region ($i=\mathrm{1},\mathrm{\dots },\mathrm{46}$, i.e. all IPCC regions with land surface) and data source ($j=\mathrm{1},\mathrm{\dots },\mathrm{25}$, i.e. 22 ESMs and three observational NBP products using an ERA-Interim climate), and p represents the year. Although El Niño and La Niña events do not produce entirely symmetrical responses in the atmosphere and the land system, we use this method for the sake of simplified results, assuming that El Niño and La Niña events are similar.

We describe the ENSO–NBP pathways using three distinct processes: ENSO strength, ENSO-induced climate anomalies, and NBP sensitivity to climate (biogeochemistry). To quantify these pathways, we calculate the ENSO-induced NBP anomalies due to the differences in each of these processes (ΔNBP^ENSO, ΔNBP^CLIM, and ΔNBP^BIO). This is done by considering the differences in only one of the processes at a time while using the mean conditions across all ESMs for the other processes. To compare the effects of similar ENSO events across the ESMs, we calculate the 90th percentile of Niño3.4 SSTAs for each ESM and for HadISST. These standardised ENSO events are referred to as “ENSO-90” here and have a return interval of around 11 years. We calculate ΔNBP^ENSO by multiplying the mean global NBP sensitivity to ENSO (${\mathit{\beta }}_{{\mathrm{EN}}_j}$) with the ENSO-90 SSTA (ΔENSO_j90) values for each model:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }{{\mathrm{NBP}}^{\mathrm{ENSO}}}_j={\mathit{\beta }}_{{\mathrm{EN}}_j}\times \mathrm{\Delta }{\text{ENSO}}_{j\mathrm{90}}.\end{array}$$

The differences due to ENSO-induced climate anomalies and NBP sensitivity to climate are assessed by fitting a multiple linear regression model of NBP (MLR_NBP) for each region and data source. These models predict NBP based on annual temperature and precipitation anomalies as follows:

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }{\text{NBP}}_{ijp}={\mathit{\beta }}_{{\mathrm{NT}}_{ij}}\times \mathrm{\Delta }{\text{Temp}}_{ijp}+{\mathit{\beta }}_{{\mathrm{NP}}_{ij}}\times \mathrm{\Delta }{\text{Precip}}_{ijp},\end{array}$$

where ΔNBP_ijp represents the NBP anomaly and its sensitivity to temperature (precipitation) is denoted by ${\mathit{\beta }}_{{\mathrm{NT}}_{ij}}$ (${\mathit{\beta }}_{{\mathrm{NP}}_{ij}}$). To assess the differences due to ENSO-induced climate anomalies, we use MLR_NBP with temperature and precipitation anomalies from an ENSO-90 event for each ESM and averaged sensitivity values across all ESMs:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathrm{\Delta }{{\mathrm{NBP}}^{\mathrm{CLIM}}}_{ij}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{22}}}\sum _{j=\mathrm{1}}^{\mathrm{22}}\left({\mathit{\beta }}_{{\mathrm{NT}}_{ij}}\right)\times {\mathit{\beta }}_{{\mathrm{ET}}_{ij}}\times \mathrm{\Delta }{\text{ENSO}}_{j\mathrm{90}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{22}}}\sum _{j=\mathrm{1}}^{\mathrm{22}}\left({\mathit{\beta }}_{{\mathrm{NP}}_{ij}}\right)\times {\mathit{\beta }}_{{\mathrm{EP}}_{ij}}\times \mathrm{\Delta }{\text{ENSO}}_{j\mathrm{90}}.\end{array}\end{array}$$

Conversely, the differences due to biogeochemistry are calculated using temperature and precipitation anomalies from an ENSO-90 event averaged across all ESMs with the following model-specific NBP sensitivities to climate:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathrm{\Delta }{{\mathrm{NBP}}^{\mathrm{BIO}}}_{ij}& ={\mathit{\beta }}_{{\mathrm{NT}}_{ij}}\times {\displaystyle \frac{\mathrm{1}}{\mathrm{22}}}\sum _{j=\mathrm{1}}^{\mathrm{22}}({\mathit{\beta }}_{{\mathrm{ET}}_{ij}}\times \mathrm{\Delta }{\text{ENSO}}_{j\mathrm{90}})\\ & +{\mathit{\beta }}_{{\mathrm{NP}}_{ij}}\times {\displaystyle \frac{\mathrm{1}}{\mathrm{22}}}\sum _{j=\mathrm{1}}^{\mathrm{22}}({\mathit{\beta }}_{{\mathrm{EP}}_{ij}}\times \mathrm{\Delta }{\text{ENSO}}_{j\mathrm{90}}).\end{array}\end{array}$$

We compare the contributions of the three processes to the overall uncertainty by measuring the spread of the ENSO–NBP relationship among the ESMs as the coefficient of variation (CV), i.e. the ratio of the standard deviation to the mean (expressed as a percentage).

Figure 2 ENSO-induced NBP anomalies across ESMs (green) and observational products (purple) for (a) an El Niño event with a 90th-percentile intensity (CV of 48 % (39 % without GFDL-ESM4)), as well as the contributions of three processes to differences in the ENSO–NBP relationship: (b) differences in ENSO strength (CV of 32 %), (c) differences in ENSO-induced climate anomaly patterns (teleconnections (CV of 27 %)), and (d) differences in biogeochemistry (CV of 56 %). Results from the TRENDY ensemble are denoted as diamonds.

Figure 3 Regional ENSO-induced NBP anomalies in ESMs and observational products from an El Niño event with a 90th-percentile intensity. Note the differences in the x-axis ranges.

3 Results

3.1 Global ENSO–NBP relationship

The NBP anomalies of a standardised ENSO-90 event range from −0.15 to −2.13 Pg C yr⁻¹, with a CV of 48 %, which we use as the reference value for the uncertainty in the ENSO–NBP relationship (Fig. 2a). The mean ESM ENSO-90 NBP anomaly has a value of −0.88 Pg C yr⁻¹, which is between the CAMS and TRENDY values. The mean ENSO-90 NBP anomaly from the observations is pulled down by FLUXCOM, which is reported to underestimate the interannual variability in carbon fluxes (Jung et al., 2019). There is strong disagreement regarding the regional contribution to global NBP anomalies (Fig. 3). This disagreement is exemplified in SE Asia (SEA) and northern South America (NSA), the two regions contributing the most to global NBP anomalies. The combined NBP anomalies of these two regions account for 48 % of the global NBP anomalies across all ESMs. However, the global contribution of SEA and NSA ranges from −23 % in MIROC-ES2H (unlike the other ESMs, MICRO-ES2H has positive NBP anomalies in SEA during El Niño events), to 70 % in NorESM2-LM. There is also little agreement regarding the ratio of SEA to NSA anomalies. Although the mean NBP anomaly from SEA is 21 % larger than the mean anomaly from NSA, half of the 22 ESMs have larger anomalies in NSA than in SEA.

Although the method used here does not take the asymmetry of ENSO events into account, we found this effect to be negligible on the global scale for most ESMs (Fig. S1 in the Supplement).

3.2 ENSO strength

ENSO strength varies between 0.35 °C (INM-CM5-0) and 1.39 °C (CMCC-ESM2), while the ESM mean (0.86 °C) is slightly higher than that of the HadISST reanalysis (0.76 °C; Fig. 4). Despite the wide variation in ENSO strength, these differences are partially offset by the sensitivity of global NBP to ENSO. ESMs with strong ENSO values tend to have a lower NBP sensitivity to Niño3.4 SSTAs.

Figure 4 ENSO strength and the sensitivity of global NBP to ENSO in 17 ESMs, FLUXCOM, CAMS, and TRENDY. The strength of the ENSO events is given as the SD of the mean annual sea surface temperatures in the Niño3.4 region (x axis). NBP sensitivity to ENSO is given as the global NBP anomaly relative to a 1 °C anomaly in the Niño3.4 region (y axis). A correlation of 0.46 between ENSO amplitude and NBP sensitivity to ENSO compensates for some of the differences between ESMs. IAV: interannual variability.

To single out the effect of the differences in ENSO strength, we multiply the ENSO-90 SSTAs by the mean NBP sensitivity (−0.86 Pg C yr⁻¹ °C⁻¹). The resulting ENSO-induced NBP anomalies range from −0.39 to −1.64 Pg C yr⁻¹. Considering only the differences in ENSO strength leads to a CV of 32 %, which is a 33 % reduction compared to the overall ENSO–NBP CV (Fig. 2b).

3.3 ENSO-induced climate anomalies

We compare the regional patterns of ENSO-induced climate anomalies between ESMs and observations and assess how these differences affect the ENSO–NBP relationship. The ESMs generally capture the sign and strength of ENSO-induced temperature and precipitation anomalies (Figs. S2 and S3). There are, however, regional differences in the uncertainty in climate anomalies among the ESMs. The spread in ENSO-induced temperature anomalies is especially high in central Australia (CAU), northern Australia (NAU), and NSA, and there are high uncertainties in the ENSO-induced precipitation anomalies in SEA and NW South America (NWS).

We use a statistical model with uniform NBP sensitivities to determine the effect of differences in ENSO-induced climate anomaly patterns on NBP (Figs. 2c and S4). Differences in ENSO-induced climate anomalies contribute less to the uncertainty in the ENSO–NBP relationship than the other two processes (CV of 27 %). Most of the global NBP anomalies fall within a similar range, except those from MPI-ESM1-2-LR. The above-average NBP anomalies in MPI-ESM1-2-LR are caused by strong ENSO-induced temperature anomalies in several regions (Fig. S2).

The comparison of the data sources used has some limitations because the reanalysis-based observation data include climate forcing, while the ESM data from piControl runs do not. However, the measured changes to ENSO teleconnection patterns remain weak and are only expected to alter by the mid-21st century (Yeh et al., 2018).

Figure 5 Differences in regional NBP anomalies based on ENSO-induced climate anomalies from ESMs and observations.

We find some notable biases in the ENSO-induced climate anomalies among the ESMs. Figure 5 shows the ENSO-90 NBP anomalies from the MLR_NBP model using mean values for NBP sensitivity and ENSO-induced climate anomalies from ESMs and observations. The strongest biases are found in SEA, where the observed ENSO-induced climate anomalies result in NBP anomalies almost twice as strong as those in the ESM climate. This bias is mostly due to stronger ENSO-induced precipitation anomalies in the observations, especially in MERRA2 (Fig. S3). Other areas with biases include eastern southern Africa (ESAF) and SE South America (SES), where ENSO-induced NBP anomalies are stronger in the observed climate than in the ESM climate, as well as the South American monsoon (SAM) region and central Africa (CAF), where the ESM climate creates stronger NBP anomalies than the observed climate.

3.4 Biogeochemistry

Differences in biogeochemistry lead to a CV of 56 %, representing a 17 % increase compared to the differences in the overall ENSO–NBP relationship (Fig. 2d). This makes biogeochemistry the largest source of uncertainty in the ENSO–NBP relationship. Indeed, considering only differences in biogeochemistry produces a larger uncertainty than the apparent uncertainty in the ENSO–NBP relationship. This highlights the compensatory effects that offset some of the differences among the ESMs.

We decomposed the ENSO-90 NBP anomalies into NPP, Rh, and fire emissions to reveal the factors driving the differences in biogeochemistry (Fig. 6). Global ENSO-induced NPP anomalies are relatively consistent, except in the MIROC ESMs. Fire emissions make up 43 % and 32 % of the ENSO-induced NBP anomalies in E3SM-1-1 and NorCPM1, respectively, while in 7 out of the 13 ESMs (with fire emissions), fire explains fewer than 4 % of the ENSO-induced NBP anomalies. ENSO-induced fire emissions originate mostly from SEA. Notable deviations among the NBP anomalies are due to uncertainties in the sign of Rh. Rh anomalies can either increase or dampen the effect of reduced NPP. In particular, the high NBP anomalies in GFDL-ESM4 result from increased Rh, most of which originates from SEA. However, there is no consistency in the role of Rh with respect to ENSO-induced NBP anomalies. While GFDL-ESM4 and UKESM1-0-LL have comparable NPP anomalies for SEA, NBP anomalies are 20 times higher in GFDL-ESM4 because the NPP anomalies in UKESM1-0-LL are offset by Rh. This demonstrates that while the NBP sensitivity to climate is suitable for describing CO₂ dynamics, it fails to capture the underlying processes.

Figure 6 The decomposition of ENSO-induced NBP fluxes (diamonds) into net primary production (NPP), heterotrophic respiration (Rh), and fire for global fluxes, Southeast Asia (SEA), and northern South America (NSA). The anomalies represent the carbon fluxes of an ENSO event with a 90th-percentile intensity. Negative values indicate reduced NPP and increased Rh and fire emissions.

3.5 SEA and NSA

Since the majority of the ENSO-induced NBP anomalies originate from SEA and NSA, we will examine the cause of these anomalies in more detail. Figure 7 shows the ENSO-induced precipitation anomalies and the sensitivity of NBP to precipitation for these two regions. The NBP anomalies in SEA are less constrained than those in NSA due to both the uncertainty in the ENSO-induced precipitation anomalies and the sensitivity of NBP to precipitation. While the differences in the ENSO-induced climate anomalies have a small impact on the uncertainty in the global ENSO–NBP relationship, we find large uncertainties on a regional scale. The ENSO-induced precipitation anomalies range in value from 100 to −210 mm yr⁻1 in SEA and from −33 to −180 mm yr⁻1 in NSA.

Although the ENSO-induced NBP anomalies from IPSL-CM6A-LR for SEA reflect the mean NBP anomalies from the other ESMs (Fig. 3), they result from the compensation of two anomalous behaviours. In IPSL-CM6A-LR, unlike in most ESMs, ENSO creates positive precipitation anomalies in SEA. However, this atypical behaviour is cancelled out by the negative sensitivity of NBP to precipitation.

Figure 7 Compositions of ENSO-induced NBP anomalies with respect to SE Asia (SEA) and northern South America (NSA). The x axis shows the ENSO-induced precipitation anomalies, and the y axis shows the sensitivity of NBP to precipitation.

These results highlight the need to describe the complete ENSO–NBP pathway to fully understand the relationship between ENSO and CO₂. The pathways describe not only the strength but also the where and why of the ENSO–NBP relationship. We demonstrate the extent of these differences qualitatively by comparing the ENSO–NBP pathways of the four ESMs that are most similar in terms of apparent ENSO–NBP strength (Table 2). The ENSO-90 NBP anomalies of these four ESMs fall within a very narrow range of −0.97 to −0.99 Pg C, which encompasses the TRENDY anomaly (−0.98 Pg C). Within these four ESMs, we focus on the role of SEA as this is the region with the strongest ENSO-induced NBP fluxes. The share of SEA in the global NBP anomaly is 9 % in ACCESS-ESM1-5, 12 % in MPI-ESM1-2-LR, 19 % in CanESM5, and 56 % in NorCPM1. We further use the MLR_NBP model to separate the SEA NBP anomalies into components caused by either temperature or precipitation. Temperature anomalies explain 74 % and 86 % of the NBP anomalies in MPI-ESM1-2-LR and ACCESS-ESM1-5, respectively. While the high NBP anomaly in NorCPM1 is exclusively caused by precipitation, the positive precipitation anomaly in CanESM5 even mitigates the overall NBP anomaly caused by temperature. Lastly, we break down whether the temperature-related NBP anomalies are caused by strong ENSO-related temperature anomalies or by a high sensitivity of NBP to temperature. Although CanESM5 has by far the highest temperature-related NBP anomalies, the actual temperature anomalies are 2 and 4 times smaller than those in ACCESS-ESM1-5 and MPI-ESM1-2-LR, respectively. The temperature-related NBP anomalies in CanESM5 can be attributed to the high sensitivity of NBP to temperature. The remaining ESMs, ACCESS-ESM1-5 and MPI-ESM1-2-LR, are the most similar in terms of temperature-related NBP anomalies. However, this is because the 2.5-times-higher ENSO-induced temperature anomalies in MPI-ESM1-2-LR are somewhat compensated for by differences in the sensitivity of NBP to temperature.

Table 2 Differences in the ENSO-induced NBP anomaly pathways for four ESMs. All models have an ENSO-induced NBP anomaly of around −0.98 Pg C yr⁻¹ per standardised El Niño event. SEA NBP refers to the ENSO-induced NBP anomaly originating from Southeast Asia (SEA), T_frac represents the fraction of NBP anomalies in SEA that can be attributed to temperature, ΔTemp denotes the ENSO-induced temperature anomaly, and β_NT represents the sensitivity of NBP to temperature.

4 Discussion

We compared the strength and characteristics of the ENSO–NBP relationship across 22 CMIP6 ESMs. The largest source of uncertainty in the simulated ENSO–NBP relationship is due to differences in biogeochemistry. Although differences in ENSO-induced climate anomalies are strong at the regional scale, these errors cancel out globally.

The sensitivity of NBP to Niño3.4 SSTAs is still poorly constrained, with SSTA values ranging from −0.13 to −2.00 Pg C °C⁻¹ (with an SD of 0.44 Pg C °C⁻¹). A possible explanation for the negative relationship between ENSO strength and the sensitivity of global NBP to Niño3.4 SSTAs could be an artefact of model tuning. Although we could not find direct evidence that the climate sensitivity of the terrestrial carbon cycle is directly tuned to ENSO, it is possible that this happens indirectly by constraining the overall variability in atmospheric CO₂. Some of the spread in the sensitivity of NBP to Niño3.4 SSTAs can be explained by differences in individual processes, such as ENSO, or the overall sensitivity of carbon fluxes to climatic drivers (Padrón et al., 2022). Although this type of error leads to a large uncertainty in the ENSO–NBP relationship, the errors introduced by a single process do not compromise the consistency of the results. Differences in ENSO strength, for example, could be balanced out with a single scaling factor.

Another type of error is based on differences in the ENSO–NBP pathways. The combined differences in ENSO-induced climate anomaly patterns and biogeochemistry lead to high uncertainty in the processes behind the ENSO–NBP relationship. Our exemplified descriptions of four ENSO–NBP pathways show that there is little agreement regarding the origin and drivers of ENSO-induced NBP anomalies, even if they result in similar CO₂ growth. The disagreements in the pathways primarily affect the estimation of regional carbon dynamics. However, the differences in the ENSO–NBP pathways can also distort the prediction of atmospheric CO₂ in initialised prediction systems. This results from the interaction of initial conditions and ESM-specific ENSO teleconnection patterns. Although ESMs might exhibit comparable relationships between ENSO and global NBP, a specific ENSO event can still result in different NBP anomalies. This can be exemplified by the pathways of ACCESS-ESM1-5 and NorCPM1 (Table 2). While both ESMs exhibit a similar average response to ENSO, the NBP anomaly is almost exclusively from SEA in NorCPM1; however, SEA does not play a large role in ACCESS-ESM1-5. Consequentially, the initial conditions and large-scale weather patterns that influence SEA interact with the ENSO-induced climate anomalies in SEA and co-determine the global NBP anomalies. An initial positive water storage anomaly in NorCMP1, for example, can mitigate the impact of reduced precipitation in SEA and limit the reduction in global NBP.

The main challenge in improving the representation of the ENSO–NBP relationship is addressing the errors in the biogeochemistry. The differences in the sensitivity of NBP to climate drive the uncertainty in the overall ENSO–NBP relationship. A large portion of this uncertainty can be attributed to the partitioning of NBP. While there is some deviation in the climate sensitivity of NPP, there is no consensus on the sign of Rh sensitivity. Under normal circumstances, NPP and Rh are positively correlated (Baldocchi et al., 2018) because, firstly, similar climatic conditions favour both types of fluxes and, secondly, the organic material required for Rh is provided through NPP. An exception to this behaviour was proposed by Zeng et al. (2005) and results from a “conspiracy” between ENSO-induced climate anomalies and plant and soil physiology. Increased temperatures and reduced precipitation during El Niño events in the tropics can limit NPP while enhancing Rh. However, this additive effect of NPP and Rh is only observed in SEA in GFDL-ESM4 and EC-Earth3-CC. In most other ESMs, the ENSO-induced reduction in NPP is partially compensated for by the negative Rh anomaly.

A smaller but still significant portion of the uncertainty can be attributed to ENSO-induced climate anomalies. While the extent of global climate anomalies does not vary much among the ESMs, there is high uncertainty in the spatial distribution of the ENSO-induced anomaly patterns. Although ENSO-induced climate anomalies affect several regions across the globe, they induce the strongest NBP anomalies in SEA and NSA. Important steps in model development include reflecting the observed balance of ENSO-induced climate anomalies for these regions and reducing the strong bias in climate anomalies in SEA.

The uncertainty among observed carbon fluxes is mostly due to the low interannual variability in NBP in FLUXCOM. Although FLUXCOM data are not suitable for estimating absolute ENSO-induced NBP anomalies, they can still be used to assess the relative contributions of individual regions to the global ENSO–NBP signal. The ENSO-induced NBP anomalies in CAMS are also below the mean ESM values. This generally lower response of land carbon fluxes to ENSO in inversion data, compared with ESMs, has been reported by others (Bastos et al., 2018, 2020). Although there are differences between inversion products based on atmospheric transport models, assimilated observations, prior carbon fluxes, and fossil fuel emission data (Gaubert et al., 2019), the differences among inversion products are small compared to the differences among ESMs (Bastos et al., 2020).

5 Conclusions

Although ESMs are able to reproduce the relationship between ENSO SSTAs and CO₂, there is little agreement with regard to the processes behind this relationship. While some of the ESMs and the multi-model mean reflect the observed ENSO–NBP pathways well, individual ESMs deviate substantially. This is because uncertainties on regional scales are balanced out when fluxes are aggregated globally. Consequentially, the correct reproduction of atmospheric CO₂ variability in assimilation runs forced by SSTAs does not necessarily indicate a good representation of atmospheric circulation patterns and biogeochemistry. It could also result from an alternative ENSO–NBP pathway that does not reflect observable processes. We attribute this high uncertainty in the ENSO–NBP relationship to our limited understanding of the sensitivity of terrestrial carbon fluxes to climate. This ongoing challenge is due to the low availability and quality of carbon flux observations, limiting the ability of ESMs to reproduce interannual variability in the terrestrial carbon cycle. This is, however, where the ENSO–NBP relationship provides untapped potential. Instead of tuning ESMs with local observations of carbon flux data, the models could be optimised to reproduce the large-scale ENSO-induced patterns. This could be a favourable alternative to the traditional approach due to the availability of high-accuracy data from atmospheric CO₂ measurements and continental carbon flux anomalies from inversion products. Improving the reproduction of the regional response of terrestrial carbon fluxes to ENSO-induced climate anomalies is not only a tangible goal but also one that may lead to ESMs with both a better ability to simulate interannual variability in global carbon fluxes and improved predictability of the global carbon cycle.