Natural decadal variability of global vegetation growth in relation to major decadal climate modes

The ongoing climate change can modulate the behavior of global vegetation and influence the terrestrial biosphere carbon sink. Past observation-based studies have mainly focused on the linear trend or interannual variability of the vegetation greenness, but could not explicitly deal with the effect of natural decadal variability due to the short length of observations. Here we put the variabilities revealed by remote sensing-based global leaf area index (LAI) from 1982 to 2015 into a long-term perspective with the help of ensemble Earth system model simulations of the historical period 1850–2014, with a focus on the low-frequency variability in the global LAI during the growing season. Robust decadal variability in the observed and modelled LAI was revealed across global terrestrial ecosystems, and it became stronger toward higher latitudes, accounting for over 50% of the total variability north of 40°N. The linkage of LAI decadal variability to major natural decadal climate modes, such as the El Niño–Southern Oscillation decadal variability (ENSO-d), the Pacific decadal oscillation (PDO), and the Atlantic multidecadal oscillation (AMO), was analyzed. ENSO-d affects LAI by altering precipitation over large parts of tropical land. The PDO exerts opposite impacts on LAI in the tropics and extra-tropics due to the compensation between the effects of temperature and growing season length. The AMO effects are mainly associated with anomalous precipitation in North America and Europe but are mixed with long-term climate change impacts due to the coincident phase shift of the AMO which also induces North Atlantic basin warming. Our results suggest that the natural decadal variability of LAI can be largely explained by these decadal climate modes (on average 20% of the variance, comparable to linear changes, and over 40% in some ecosystems) which also can be potentially important in inducing the greening of the Earth of the past decades.


Introduction
The terrestrial biosphere plays an important role in the global carbon cycle by absorbing about one-third of the current anthropogenic carbon dioxide (CO 2 ) emissions (Friedlingstein et al 2022). Understanding the global natural vegetation response to climate change is essential to evaluate its ability as a carbon sink (Tagesson et al 2020), which is critical for future sustainable ecosystem services and management (Keenan and Williams 2018). Decades of satellite remote sensing records make it feasible to investigate how global vegetation responds to the changing climate on a global scale.
There has been considerable interest in monitoring and understanding the linear trend of global and regional vegetation greenness using observationbased evidence (see a review by Piao et al (2020)). For instance, the general greening trend of the Earth from 1982 to 2009 was reported, mainly explained by CO 2 fertilization effects (Zhu et al 2016). The contribution of climate change to the greening was principally attributed to anthropogenic warming and regional linear trends in precipitation (Piao et al 2020).
Interannual variability (IAV) of global vegetation greenness has also been a focus, which serves as an important indicator of ecosystem vulnerability or resilience to external disturbances, in addition to long-term trends and averages (see a review by Thornton et al (2014)). Earlier work (e.g. Keeling et al 1995;Zeng et al 2005) highlighted the terrestrial processes in driving the global carbon cycle. Braswell et al (1997) pioneered investigating the response of global vegetation dynamics to temperature anomalies on interannual timescales, and biome-specific and response time-specific relations were carefully considered. Follow-up studies usually include longer records of vegetation index derived from satellite remote sensing (e.g. Fensholt et al 2012) and more climate variables such as precipitation (e.g. Zeng et al 2013) and vapor pressure deficit (He et al 2022); or focus on the IAV of the global vegetation greenness itself: it was found to be intensified during 1982-2015 (Chen et al 2019a). Other studies offer insight into the link between IAV of ecosystems and climate modes through extremes (Reichstein et al 2013, Liu et al 2017, Bowman et al 2020, Byrne et al 2021. Natural IAV is mainly contributed by the El Niño-Southern Oscillation (ENSO), the most prominent year-to-year source of IAV that takes place in the tropical Pacific (figure S1(a)) (Timmermann et al 2018, McPhaden et al 2020 and impacts the global climate through atmospheric teleconnections (Yeh et al 2018).
Climate's natural decadal or multi-decadal variability is attracting increasing attention (see a review by Deser et al (2010)). Most of these climate modes arise in the Pacific and the Atlantic, notably the ENSO decadal variability (ENSO-d) (e.g. An and Wang 2000, Sun and Yu 2009, Park et al 2020, the pan-Pacific mode of the Pacific decadal oscillation (PDO) (Mantua et al 1997), and the Atlantic multidecadal oscillation (AMO) (Kerr 2000). The ENSO-d results from the decadal modulations of ENSO's characteristics, such as amplitude, frequency, skewness, and flavor (spatial pattern). The PDO is characterized by anomalous cooling in the central and western North Pacific and horseshoe-shaped warming along the west coastline of North America in its warm phase (figure S1(b)). The PDO phase change was suggested to exert considerable impacts on marine ecosystems and climate conditions (Newman et al 2016). The AMO is characterized by slower sea surface temperature (SST) oscillations in the North Atlantic (figure S1(c)), and it was suggested to modulate temperature and rainfall over much of the Northern Hemisphere, in particular in North America and Europe during summer (Zampieri et al 2017). These decadal climate modes oscillate at longer timescales than seasonal or interannual climate variability. Given the fundamental impacts on the global climate by these climate modes, it is reasonable to assume that they further lead to decadal timescale variations in global vegetation growth, exerting lasting effects. In particular, the AMO shifted from a cold to a warm phase in the 1990s, which continued into the 2010s (figure S2(c); e.g. Frajka-Williams et al 2017), and its effects can be mixed with other long-term climate change signals. Previous attempts have been made to elucidate decadal climate forcings on global-scale terrestrial ecosystem processes, but are limited by ecosystem modeling approaches, and the coincidence of phase changes of these climate modes and anthropogenic climate changes were not taken into account (e.g. Ito 2011, Zhu et al 2016, Zhang et al 2018, Park et al 2020. In this study, using an observation-based global leaf area index (LAI) spanning the period from 1982 to 2015, we aim to systematically examine the natural decadal variability (ENSO-d, PDO, and AMO) in global vegetation and its linkage to the major decadal/multidecadal climate modes that originated in oceans. This period is characterized by widespread vegetation greenness that has been increasing globally. The results are compared with a set of historical (1850-2014) Earth system model (ESM) simulations to subtract the impacts of the linear climate changes on vegetation greenness and to estimate the robustness of the natural variability and the contribution to changes in global vegetation by the modulation of climate modes.

Global datasets of LAI
We used the LAI3g product with 1/12 • spatial resolution and half-month temporal resolution from 1982 to 2015 (Zhu et al 2016). This product is derived from the third-generation Global Inventory Modeling and Mapping Studies normalized difference vegetation index (NDVI3g) using an algorithm based on feedforward neural networks. LAI3g has been widely used and contributes to studies of the greening/browning trend across the globe (Zhu et al 2016, Chen et al 2019b as well as vegetation activities related to IAV (Piao et al 2014).
We used integrated growing season LAI anomalies to represent vegetation growth influenced by decadal climate modes. The growing season LAI anomaly is (c), (d) the contribution (fraction) of decadal variability to total variability (interannual + decadal), for LAI3g data and ESM simulation ensemble mean, respectively. Growing season is defined as the 6 months with maximum LAI in the mean seasonal cycle, or 12 months if minimal LAI is larger than 3 (considered as evergreen plants). estimated by subtracting the mean LAI seasonal cycle from the LAI time series for each grid cell. The growing season is defined as the 6 month maximum LAI in the mean seasonal cycle or the 12 month maximum LAI if there are evergreen plants in the grid cell (considered when the minimal climatological  LAI is greater than 3).

Ensemble ESM simulations
Since observation-based LAI3g covers 34 years, which captures relatively few phase changes of the decadal climate modes, we used ESM simulations to represent larger sample numbers of these decadal oscillations. The model we used is a state-of-the-art ESM EC-Earth, version 3.3 (Zhang et al 2021, Döscher et al 2022, which is part of the Coupled Model Intercomparison Project Phase 6. EC-Earth 3.3 comprises model components of the atmosphere, ocean, sea-ice, land surface, and dynamic vegetation. The vegetation component LPJ-GUESS (Smith et al 2014) is integrated into EC-Earth, and simulated vegetation parameters are determined by intrinsic vegetation dynamics and meteorological forcings such as temperature, precipitation, and radiation from the coupler. In the ESM configuration (EC-Earth3-Veg), which is ideal for simulating climate variability and its interactions with the global vegetation, we expect that the model can reveal the statistical robustness of such physical relationships at decadal and longer timescales.
The model has a standard resolution of ∼0.7 • (∼85 km) for the atmosphere, land, and vegetation and ∼1 • (∼100 km) for the ocean and sea ice. The atmosphere and ocean have 91 and 75 vertical layers, respectively. The simulation was initialized from a pre-industrial equilibrium state and then forced with transient boundary conditions from 1850 to 2014 (EC-Earth Consortium 2019). There is a total of eight model ensembles with the same setup. Each ensemble simulation can be regarded as capable of producing decadal climate variability (figure S3; figures 1(d) and 3(a)).

Decadal variability of LAI
We use standard deviation to quantify the LAI variability at different timescales. Changes in the standard deviation can reflect the changes in the amplitude of the climate mode in growing season LAI variations. For each grid cell, the LAI variability is calculated as the standard deviation of growing season LAI anomalies , so it is the standard deviation of 34 year values. The LAI anomalies are derived from the original LAI monthly times series removing the mean LAI seasonal cycle (over the whole period), then averaged for the growing season (to annual resolution). Then a running mean filter (7 years) is applied to separate the interannual (shorter than 7 year) and decadal (longer than 7 year) timescale variability.
We examined both detrended and original growing season LAI anomalies from LAI3g data because the data period coincided with the phase change of AMO. Thus the long-term anthropogenic global warming effect and AMO-induced warming effect over this period are mixed. For EC-Earth3-Veg model data that covers 100 year (1915-2014 is used) climate and global vegetation evolution, we always used detrended LAI time series and climate mode indices.

Climate variability mode indices
Three major decadal/multidecadal climate variability modes were considered, ENSO-d, PDO, and AMO. These climate modes are indicated by simple annual indices averaged over monthly data. The following method was applied to SST observations (extended reconstructed sea surface temperature (ERSST) v5 (Huang et al 2017)) and simulated SST to retrieve these indices. For ENSO-d, a standard ENSO index was first calculated as the leading empirical orthogonal function (EOF) of detrended monthly SST anomalies in the Niño-3.4 region (5 N-5 S, 170 W-120 W). Then, a running mean filter of 7 years was applied to the ENSO index to remove interannual ENSO variability (generally regarded as 1.5-7 years (Timmermann et al 2018)). The rest of the time series is labelled as ENSO-d. This index describes the tropical Pacific Ocean decadal oscillations between El Niño-like and La Niña-like states at a decadal time scale. For PDO, we used a common definition, which is the leading EOF time series of monthly SST anomalies over the North Pacific (20 N-70 N, 120E-110 W) after removing the global mean SST anomaly (Deser et al 2010). The AMO was derived from monthly SST anomalies after removing the global mean SST anomaly in the North Atlantic (20 N-70 N, 80 W-0 W) (Zhang et al 2021). All climate mode indices were averaged to an annual resolution to be consistent with that of the growing season LAI.

Attribution of decadal climate modes (multiple regression)
Multiple regression (equation (1) is used to attribute decadal timescale LAI variability to ENSO-d, PDO, and AMO (in bold in equation (1). The regression is done for both detrended and original growing season LAI3g LAI anomalies and climate mode indices, and only for detrended modelled LAI anomalies and climate mode indices. These variables are z-scored before the regression (1)

Notes on data-model comparison
We used a group of ensemble simulations to isolate the internal model variability from the forced climate response induced by the boundary conditions (Kay et al 2015), such as greenhouse gas levels, aerosol concentration, and land-use changes. It is difficult for the model to reproduce consistent spatial patterns (climate modes themselves as well as the atmospheric teleconnections) and to have a case-by-case comparison (e.g. for positive/negative phases of climate modes in the same period of 1982-2015) due to model internal variability. Therefore, our quantitative data-model comparison of the global vegetation decadal variability focused on LAI standard deviation, and its relative change, rather than its absolute value. We also discussed the decadal fraction of the total (interannual + decadal) LAI variability, or the overall contribution of decadal LAI variability from climate modes (figures 1(d), 3(a) and 2(d)-(f)). ESM simulation ensemble mean being capable of reproducing the meridional distribution of decadal LAI variability (figure 3(a)) lends credibility to the use of simulated climate forcing processes on the vegetation dynamics by the decadal intrinsic climate modes. The advantage of the ensemble historical simulations is the length (100-year output can contain five complete 20-year oscillations), and the possibility of an estimate of model internal variability for robust statistical results using the ensemble mean and spread.
We rely on the simulations to identify decadal variability in simulated growing season LAI anomalies and to quantify the forcing effects of simulated decadal climate modes associated with ocean processes. The regression coefficients and R2 showed a noticeable ensemble spread (figure not shown). Nevertheless, the explained variance of multiple regression (its ensemble mean) provides an estimate of the pure effects of decadal climate modes on vegetation greenness in the model, since both linear trends and IAV are removed. In principle R2 of ENSO-d, PDO and AMO cannot be separated, either for observation/reanalysis or ESM simulations. But a direct comparison of the two enables the estimate of R2 of AMO. By subtracting R2 derived from detrended observation data (LAI3g, including ENSO-d, PDO effects; figure 2(e)) from R2 of the ESM ensemble mean (ESM LAI, including ENSO-d, PDO & AMO effects), the contribution of AMO can be quantified (figure 2(f)); by subtracting R2 derived from undetrended observation data (LAI3g, including decadal as well as linear climate forcing effects) from R2 of ESM ensemble mean, the contribution of linear climate change can be quantified (figure 2(d)).

Composite analysis of climate variables based on climate modes
We analyze the impacts of surface air temperature (National Centers for Environmental Prediction (NCEP) reanalysis 2; Kanamitsu et al 2002) and precipitation (Global Precipitation Climatology Project (GPCP) data; Adler et al 2018) forcing associated with decadal climate modes on global vegetation growth by composite analysis (figure 4). We composite temperature and rainfall anomalies during the vegetation growing season over each land grid cell by positive and negative phases of each climate mode. These positive and negative phases are defined by using the 80th and 20th percentiles of the undetrended (otherwise the AMO phase shift is removed) ENSO-d, PDO and AMO indices (see figure S2) to ensure more reprehensive phase changes of climate modes. The relative importance of temperature and precipitation forcing effects can be quantified by comparing how close the climate variable follows the influences of climate modes using their correlations. To take ENSOd as an example, if the log of correlation coefficient ratio of corr(temperature anomaly, ENSO-d index) and corr(precipitation anomaly, ENSO-d index) is positive, the former correlation is higher which indicates potentially stronger temperature forcing effects.

Identification of decadal signals in LAI
Both interannual and decadal timescale variabilities of growing season LAI anomalies in the satellite observations are prominent across the global land ecosystems (figures 1(a) and (b)). Larger variability generally appears in areas with high vegetation productivity, such as tropical forests. The peak of variability in the tropics decreases toward higher latitudes in the north and south, reaching minimums in the midlatitude deserts and then rebounding in the temperate and boreal forest regions (figures 1(a) and (b)). The decadal contribution to total LAI variability shows that decadal variability is relatively more important at the mid-and high-latitudes (figures 1(c) and (e) red line). On average, the decadal variability (standard deviation) contributes to over 40% of the total variability, and over 50% for the majority of regions north of 40 • N (figures 1(c) and (e) red line). In particular, the high share of decadal-scale variability in the total variability is dominated by linear trends over the northern Eurasian continent and in South and East Asia (figures S4 and S5).
The relatively stronger decadal variability towards higher latitudes is also pronounced in the model ensembles (figures 1(d) and 3(a) black line). This feature is in line with the centers of action of the major decadal climate modes of PDO and AMO in the North Pacific and North Atlantic. Interestingly, however, simulated decadal variability is notably stronger than satellite observations (detrended) in the Northern Hemisphere north of 40 • N (figure 3(a) the black line vs. the red line; figure S4(c)), but in a better agreement at Northern Hemisphere (NH) higher latitudes when the linear trend was kept in the observations (figure 3(a) the black line vs. the blue line; figure S4(d)). Indeed, there is a robust greening trend in this region as seen in the observations which can dominate the decadal variability if the linear trend is included ( figure S5(c)). This greening is associated with long-term impacts (can be considered as a linear trend, e.g. CO 2 fertilization effects, anthropogenic warming, land-use changes; figure S6) or the AMO. Coincidently, the AMO shifted from a negative phase to a positive phase throughout the study period of 1982-2015 (figure S2(c)), so when the detrending was applied to the satellite data, the impacts of AMO on LAI variability north of 40 • N may be diminished. This situation is different for low-latitudes and Southern Hemisphere vegetation decadal-scale changes as they have a closer agreement with detrended data (figure 3(a) the black line vs. the red line). This suggests that they are more closely related to ENSO-d and PDO, and are thus not affected by the detrending.

Investigation of the linkage of LAI decadal variability to climate modes
The global IAV in the terrestrial biosphere (Buermann et al 2003) can be attributed to the remote impacts (in contrast to local impacts in the tropical Pacific) of ENSO that perturb the global climate through atmospheric teleconnections (Liu and Alexander 2007). It is also influenced by photosynthesis responses to changes in precipitation between years (Thornton et al 2014, Ahlström et al 2015 or by a smaller magnitude, biogeophysical soil processes such as legacy effects (Braswell et al 1997, Wong et al 2021 that last a few years. However, for periodic variability with longer timescales, ocean dynamics is likely the dominating cause of decadal climate variations (Gill and Adrian 1982). Next, we analyzed the decadal LAI variability to its specific triggers and drivers in the atmosphere and ocean.
We used multiple regression analysis to reveal the global impacts of three major decadal climate modes (figure 2). The LAI anomaly pattern regressed on the ENSO-d index depicts an overall reduction in LAI during positive El Niño-like SST anomalies in the tropics stretching across the Amazon, equatorial Africa, and the maritime continent ( figure 2(a)).
The vegetation response appears more robust in the tropics (negative LAI anomalies in East Amazon, equatorial Africa, and the maritime continent in figures 2(a) and (e)). This is consistent with El Niño impacts at an interannual timescale reported in a previous study (Betts et al 2020). Regression coefficients also suggest that positive El Niño-like climate anomalies lead to positive LAI anomalies in Central North and South America, East Asia, and eastern Australia.
The decadal vegetation response to positive PDO modulation is the weakest among the studied climate modes. It is relatively robust in the circum-Pacific and equatorial African regions, with positive responses in western Amazon, central America, central Africa, and the maritime continent and negative responses in the western US, eastern Amazon, South Africa, and East Australia. Note that the differences in the regression of undetrended ENSO-d or PDO are not remarkable (figures S2(a) and (b)).
The spatial pattern of the regression coefficient of the LAI decadal anomalies versus AMO is heavily dependent on the detrending of the AMO index (figure not shown). As mentioned above, the AMO shifted to a positive state during 1982-2015, inducing North Atlantic basin warming; so, it could be mixed with the global warming linear trend during this study period. If the index was not detrended and the phase change of AMO was retained, the regression pattern suggests favorable conditions for vegetation growth in Europe, eastern North America, Sahel, South Africa, and India under positive AMO, although these effects can be partly contributed by the linear trend in climate forcings (figures 2(c) and S6).
The combined effects from these climate modes as represented by simple indices can be separated from those long-term (linear) impacts with the help of ESM simulations (section 2.6). It is found the decadal climate modes (figures 2(e) and (f)) and the linear trend (figure 2(d)) have roughly comparable effects on LAI decadal variability (figures 3(b) and (c)). Together they explained more than 50% of the decadal variability (i.e. all variability longer than IAV) in East Amazon, South Africa, East Australia, and Indonesia, approximating 20% in all regions with relatively large LAI decadal variability (figures 2(e) and (f)). In East Asia, Indonesia, East Amazon, North America, and Europe, LAI decadal variability is more dominated by the decadal climate modes. In Sahel and India, the long-term linear trend plays a more important role. Most of these regions were identified as greening hotspots in recent decades (Piao et al 2020).
Earlier studies considered lagged responses of the ecosystem to climate change at interannual timescales in correlation and regression analyses (e.g. Braswell et al 1997, Zhou et al 2003. We also performed multiple regression with a time lag, assuming the decadal climate modes led the vegetation response by 1 year. However, global R2 generally decreased ( figure S7). This implies that the LAI decadal variability is dominated by a zero-lag modulation of the climate modes. Previous studies reported that optimal lag ranges from 0 to 24 months between various climate conditions and the vegetation or carbon flux response (e.g. Ito 2011, Zhu et al 2017. In our analyses, because of the annual-resolved time series based on the growing season average, which already involved potential lag effects of 1 year or less, the effects of further delayed responses are not robust and will not be further considered.

Understanding the forcing mechanisms of climate modes
To understand how different phases of climate modes affect vegetation growth, we conducted a composite analysis for the surface air temperature and precipitation anomalies during the vegetation growing season over each land grid cell, and quantify the relative importance of temperature and precipitation (figure 4). The composite temperature and precipitation anomalies overall demonstrate distinctive opposite changes during positive and negative phases of these decadal climate modes (figure S2), highlighting the driving effects of decadal timescale climate oscillations. Note that other potentially important climate indicators that might be more representative of regional climate and have strong impacts on specific ecosystems such as vapor pressure deficit impacts on tropical forests (Barkhordarian et al 2019, Worden et al 2021 are not considered in this study. Similarly, complex interactions and feedbacks between terrestrial ecosystems and climate (e.g. Wang et al 2020) involving ecosystem disturbances (such as climate extremes, wildfires, land-use) are beyond the scope of this study.
The significant negative tropical LAI anomalies can be attributed to droughts associated with the positive phase of ENSO-d (figures 2(a), 4(b1), (b2) and (c1)), mostly robust in the eastern Amazon and maritime continent. In comparison, the temperature effects show a smaller impact, as both cooling and warming were associated with lower LAI, such as in equatorial Africa and western Amazon. Increased precipitation in positive ENSO-d also contributed to enhanced growing season LAI in central North and South America, East Asia, and northwestern Australia (figures 4(b1), (b2) and (c1)).
Temperature effects were more marked and consistent in PDO modulation of mid/high latitude LAI. For instance, cooling (warming) in Southeast Australia (western coast of Canada) could shorten (extend) the growing season, leading to negative (positive) LAI (figures 4(a3), (a4) and (c2)). In the tropics (e.g. equatorial Africa), the higher LAI can be induced by cooler temperatures and reduced heat stress (Kim et al 2017), while the impacts of dry conditions were not obvious.
For AMO, temperature and precipitation composites show sub-continental-scale changes (figures 4(a5), (a6), (b5) and (b6)). However, they were superposed by the long-term linear changes in temperature and precipitation during 1982-2015 (figures S5(a) and (b)) due to the coincident phase shift of AMO. Model simulations suggest that there might be less distinctive decadal variability associated with temperature changes compared to precipitation changes (figures S8(e) and (f)). Earlier studies documented decreased summer (which can be approximately regarded as the growing season here) precipitation and warmer temperatures in the United States and increased summer precipitation in Europe during positive AMO phases (Sutton and Hodson 2005, Lyu and Yu 2017, O'Reilly et al 2017, which is partly consistent with our study. The highly positive impacts of warm AMO phases on vegetation growth, which are robust in western and eastern Europe, the Sahel, South Africa, and India, can be attributed to more rainfall in these regions, and warmer temperatures in eastern Canada (figure 4(c3)).

Summary and discussion
Our results indicate the importance of decadal timescale changes in global LAI, with the increasing magnitude of anomalies toward higher latitudes. These changes were modulated by climate modes embedded in natural slow ocean processes. Quantified by simple indices, the ENSO-d, PDO, and AMO together explained on average 20% of the decadal variability in LAI and ∼40% for parts of northern North America, the eastern Amazon, Europe, South Africa, Northeast Australia, India, and Indonesia. The oscillations of the climate modes exerted overall distinct climate conditions across the globe and induced vegetation greening/browning. The relative importance of climate forcing from temperature and precipitation anomalies depends on the specific climate mode and the local climate. For ENSO-d, large parts of tropical land are characterized by a precipitation-vegetation relationship. For PDO modulation, cooling effects in the tropics and extra-tropics result in opposite vegetation responses, attributed to the reduced negative impact of warming (e.g. heat stress and shortened growing season growth). The AMO effects are more complex as they are mixed with the long-term linear trends in our 34 year study period, and increased precipitation in Canada, western and eastern Europe, the Sahel, South Africa, and India induced vegetation greening in these regions. Studies on the impacts of ongoing climate changes on global vegetation growth should consider not only the linear trend and IAV but also the natural decadal climate modes. Our study improves the understanding of the linkages between these climate modes and the terrestrial ecosystems and has useful implications for long-term ecosystem management and services.

Data availability statement
All EC-Earth3-Veg historical simulations are available from ESGF (Earth System Grid Federation) nodes. More details of these simulations can be found at the following URL/DOI: https://doi.org/10.22033/ESGF/ CMIP6.4706.