Quantifying Carbon Cycle Extremes and Attributing Their Causes Under Climate and Land Use and Land Cover Change From 1850 to 2300

The increasing atmospheric carbon dioxide (CO2) mole fraction affects global climate through radiative (trapping longwave radiation) and physiological effects (reduction of plant transpiration). We use the simulations of the Community Earth System Model (CESM1‐BGC) forced with Representative Concentration Pathway 8.5 to investigate climate‐vegetation feedbacks from 1850 to the year 2300. Human‐induced land use and land cover change (LULCC), through biogeochemical and biogeophysical processes, alter the climate and modify photosynthetic activity. The changing characteristics of extreme anomalies in photosynthesis, referred to as carbon cycle extremes, increase the uncertainty of terrestrial ecosystems to act as a net carbon sink. However, the role of LULCC in altering carbon cycle extremes under business‐as‐usual (continuously rising) CO2 emissions is unknown. Here we show that LULCC magnifies the intensity, frequency, and extent of carbon cycle extremes, resulting in a net reduction in expected photosynthetic activity in the future. We found that large temporally contiguous negative carbon cycle extremes are due to a persistent decrease in soil moisture, which is triggered by declines in precipitation. With LULCC and global warming, vegetation exhibits increased vulnerability to hot and dry environmental conditions, increasing the frequency of fire events and resulting in considerable losses in photosynthetic activity. While most regions show strengthening of negative carbon cycle extremes, a few locations show a weakening effect driven by declining vegetation cover or benign climate conditions for photosynthesis. Increasing hot, dry, and fire‐driven carbon cycle extremes are essential for improving carbon cycle modeling and estimation of ecosystem responses to LULCC and rising CO2 mole fractions. Moreover, large aberrations in vegetation productivity represent potential and growing threats to human lives, wildlife, and food security.

Prior studies (Bonan, 2015;Reichstein et al., 2013;Zhu et al., 2016) have found that the combined effect of CO 2 fertilization, increasing temperature, nitrogen deposition, and lengthening of growing seasons lead to increased vegetation productivity and strengthening of carbon sinks with a negative feedback to climate change. However, coupled carbon-climate models have large uncertainties in future projections of ecosystem responses and feedback strength (Frank et al., 2015;Hoffman et al., 2014;Ichii et al., 2005;Reichstein et al., 2013). Hoffman et al. (2014) found persistent atmospheric CO 2 biases in Coupled Model Intercomparison Project 5 (CMIP5) models because of uncertainties in biological and physical processes related to carbon accumulation. While most Earth system model-based climate change studies analyze projections till year 2100, these projections may miss physical-biogeochemical feedbacks that arise later from the cumulative effects of climate warming (Moore et al., 2018). The negative sensitivity of terrestrial carbon cycle to rising temperature will likely have increasing adverse implications on carbon cycle extremes over time (Frank et al., 2015;Hubau et al., 2020). Understanding the direction and strength of these feedbacks is essential for estimating long-term CO 2 concentrations and predicting and mitigating the impact and extent of climate change. These limitations could alter the assessment of the rate of increase of atmospheric CO 2 and intensity of associated feedbacks with the terrestrial biosphere.
While the effects of increased warming due to GHGs are spatially extensive, the LULCC effects are more regional (Pitman et al., 2012). The land-use history reconstruction used in this study estimated the proportion of land surface impacted by human activities to be 60% of total vegetated area, mainly due to conversions from primary vegetation to managed vegetation by 2100 (Hurtt et al., 2011). The conversion of land from natural to managed ecosystems reduces the carbon sink and its capacity to uptake anthropogenic CO 2 and influences climate by modifying biogeophysical and biogeochemical processes (Bonan, DeFries, et al., 2012). Changes in the plant functional type (PFT) at any grid cell modifies the distribution of above and below ground carbon (Oleson et al., 2010), response to light and energy (Bonan et al., 2011), distribution of soil organic matter and nutrients (Koven et al., 2013). Human activities, such as the conversion of forests and grasslands to agricultural land and urbanization, alter net radiation, sensible and latent heat partitioning, biogeochemical cycles, and the hydrologic cycle. Reduction of temperate vegetation cover by deforestation increases the albedo of the surface which decreases the net radiation that drives surface cooling and reduces evapotranspiration that may result in declines in precipitation; but tropical deforestation for pastures decreases the total atmospheric column heating and atmospheric vertical motion which leads to a warmer and drier climate (Bonan, 2015;Bonan, DeFries, et al., 2012). Since interannual variability (IAV) in GPP is strongly influenced by interannual variations in radiation, temperature, and precipitation (Ichii et al., 2005), the impact of LULCC in addition to CO 2 on carbon cycle extremes will likely increase over time.
Recent studies have investigated the characteristics of extreme anomalies in GPP due to climate change until the year 2100 (Flach et al., 2020;Frank et al., 2015;Reichstein et al., 2013;Zscheischler et al., 2014) and a few concluded that losses in carbon uptake due to negative extremes in GPP are compensated by increased CO 2 fertilization (Reichstein et al., 2013;Zscheischler et al., 2014). However, to our knowledge, no study has examined the extreme anomalies in carbon cycle beyond 2100 and the role of human LULCC in modifying carbon cycle extremes.
Rising CO 2 and LULCC impact many components of the terrestrial carbon cycle, such as total ecosystem carbon, net biome productivity, net ecosystem productivity, net primary productivity, and GPP. The overarching goal of this study is to investigate the role of CO 2 and LULCC in modifying the characteristics of extremes in GPP, referred to as carbon cycle extremes in this paper, and attribute changes to individual and compound effects of climate drivers. We hypothesize that (a) rising CO 2 emissions will lead to larger increases in the intensity, frequency and duration of negative carbon cycle extremes than positive extremes; and (b) LULCC activities in addition to CO 2 emissions will further increase GPP interannual variability and magnitude of carbon cycle extremes though total GPP will reduce. We performed analysis to (a) examine the magnitude, duration, frequency and spatial distribution of negative and positive carbon cycle extremes; (b) investigate the lagged response of antecedent climate drivers (and their compound effect) that drive carbon cycle extremes; (c) analyze the climate conditions that trigger long duration temporally continuous extremes; and (d) inspect regional changes in climate-carbon feedbacks for the Central and South American tropics.

Data
We used simulations of the Community Earth System Model (version 1) with biogeochemistry enabled, CESM1(BGC), at approximately 1° × 1° resolution to analyze climate-driven extreme anomalies in total photosynthetic activity. CESM1(BGC) is a fully coupled global climate model composed of land, atmosphere, and ocean components (Lawrence et al., 2012;Lindsay et al., 2014). The atmospheric CO 2 forcing time series consisted of the historical , Representative Concentration Pathway 8.5 (RCP 8.5;2006-2100, and Extended Concentration Pathway (ECP 8.5; 2101-2300) mole fractions, which increased from 285 ppm in 1850-1962 ppm in 2,300 ( Figure 1a). Analysis of the CMIP5 models using the International Land Model Benchmarking (ILAMB) (Collier et al., 2018) show (Table S1 in Supporting Information S1) that CESM1(BGC) is one of the best performing model in terms of overall and IAV benchmark scores when compared to observational benchmarks.
We analyzed terrestrial carbon cycle extremes (or GPP extremes) using two simulations, namely, with and without LULCC. In the simulation without LULCC, the land cover was fixed at pre-industrial (year 1850) values. In the with LULCC simulation, transient land cover was prescribed over the historical and RCP 8.5 period (1850-2100) (Hurtt et al., 2011) and consists of the prescribed spatial distribution of PFTs (Lawrence et al., 2012). Land-use conversion is assumed to stop at the year 2100, and the distribution of PFTs remains constant at year 2100 level for the period 2100 -2300, while wood harvest is maintained at a constant rate over that period (Mahowald et al., 2017). The reactive nitrogen deposition followed the spatially variable time series from 1850 to 2100 (Lamarque et al., 2010) and was subsequently held constant from 2100 to 2300. Biogeochemical processes on land and in the ocean respond to the historical and prescribed RCP 8.5 and ECP 8.5 atmospheric CO 2 concentration forcing (Figure 1a). Increasing CO 2 fertilization, water-use efficiency, and lengthening of growing seasons enhance total photosynthesis and gross primary production (GPP). These processes will be further altered by the changes in PFT distribution under with LULCC. Figure 1c shows the 5 years moving average of global annually integrated GPP for simulations with and without LULCC, both forced with the same atmospheric CO 2 concentration trajectory (Figure 1a). The magnitude of the average GPP with LULCC is less than the average GPP without LULCC potentially due to the conversion of primary vegetation to managed vegetation.

Methods
A significant deviation from the mean value is called an extreme. Extreme events in GPP signify large variations in photosynthetic activity, with positive extremes representing increases in carbon uptake, while negative extremes being indicative of loss of carbon uptake (Zscheischler et al., 2014). While extremes in climate have been extensively studied (Ban et al., 2015;Flach et al., 2020;Frank et al., 2015;Reichstein et al., 2013;Seneviratne & Hauser, 2020), few studies have focused on extremes in GPP and their underlying drivers (Flach et al., 2020;Xu et al., 2019;Zscheischler et al., 2014). In this study, we identify extremes in global carbon uptake by computing percentile-based thresholds (Seneviratne et al., 2012). Described in detail in Section 3.1, a positive (or negative) extreme in GPP is defined as anomalies in GPP that are greater (or less) than a selected percentile-based threshold of GPP anomalies. The thresholds were computed by selecting 1st and 99th percentile anomalies in GPP, calculated at every land grid cell for each consecutive 25 years from the year 1850 through 2300. The GPP extreme events are then defined as values of GPP anomalies above q (positive extremes) or below − q (negative extremes). The schematic Figure 2 illustrates the steps to compute GPP extremes.
Extremes in GPP could also be categorized as isolated versus temporally continuous extreme events. The isolated and temporally continuous extremes in GPP are analogous to hot days and heatwaves. Extreme events in GPP that are continuous in time represent a significant cumulative impact on carbon uptake. Akin to the definition of temporally continuous heatwaves by Baldwin et al. (2019), we define a temporally continuous extreme (TCE) event ( Figure S1 in Supporting Information S1) in GPP as (a) three or longer months (equal to a season length) of GPP extreme anomaly occurring consecutively, (b) and are considered as single continuous event through any gaps of 2 months or shorter in duration (assuming that ecological recovery is unlikely in less than a season) beyond 3 months. A GPP TCE event that occur after a season, that is, 3 months or more, is considered a separated TCE. Section 3.2 illustrates the method for attribution of GPP TCEs to individual and compound effects of climate drivers.

Identification and Detection of GPP Extreme Events
The time series of GPP at any grid cell consists of a trend, seasonality, and anomalies ( Figure 2) components. We used singular spectrum analysis, a non-parametric spectral estimation method based on embedding a time series in vector space (Golyandina et al., 2001), to extract signals with specific frequencies. The trend in GPP at any grid cell captures long term change in mean GPP, which is influenced by long term changes in climate drivers, atmospheric CO 2 concentration, and LULCC. Since El Niño-Southern Oscillation (ENSO), and to a lesser extent other large-scale drivers of climate variability, enhance the variability in terrestrial photosynthesis, the nonlinear trend at each grid cell was calculated by adding all the frequencies of 10 years and longer. Hence, carbon cycle extremes in our study include the GPP anomalies induced by ENSO which is one of the largest modes of climate variability and peaks about every three to 7 years (Chylek et al., 2018) and exerts strong regional effects on the terrestrial carbon cycle (Poveda et al., 2001). The seasonality in GPP follows a periodic cycle of 12 months. The conventional way to compute the annual cycle is to determine the mean climatology over a period, however, the climatology does not reflect the intrinsic nonlinearity of the climate-carbon feedback, especially under external forcing (Wu et al., 2008) because of the increased modulation of the annual cycle in GPP under business-asusual rising CO 2 emissions. We calculated the modulated annual cycle, which allows the annual cycle to change from year to year and consists of signals with a frequency of 12 months and its harmonics. The anomalies at each grid cell were calculated by subtracting the trend and modulated annual cycle from the GPP time series (Figure 2). Hence, GPP anomalies comprised of the high-frequency signals (<12 months) and the interannual variability (>12 months and <10 years). We examined the probability distribution function (PDF) of GPP anomalies for all land grid cells for every 25-year time window from 1850 through 2300. While the left tail of the PDF of GPP anomalies represents large losses in carbon uptake, the right tail portrays large gains in carbon uptake. We chose the 1st and 99th percentile of all GPP anomalies as thresholds in the 25-year time windows to identify large losses and gains in carbon uptake.
The percentile-based thresholds (1st and 99th) for every 25-year time window yield the negative (Th−) and positive (Th+) threshold trajectories of GPP anomalies. The extreme anomalies in GPP range from 425 GgC/ month for the period 1850-1874 to 840 GgC/month for 2275-2299, as shown in Figure 3. Since GPP anomalies are a subset of GPP, increasing CO 2 fertilization, water use efficiency, and lengthening of growing season lead to increases in GPP ( Figure 1c) and the thresholds of GPP anomalies. The simulations with and without LULCC, show higher magnitudes of Th− than Th+, indicating that negative anomalies in GPP are stronger than positive for the same percentile. Also, the values of the thresholds for simulation with LULCC were higher than without LULCC. To enable the comparison of simulations with and without LULCC, we selected one threshold trajectory, without LULCC Th+ (q), to apply for calculations of all positive and negative extremes in the current study. Thus, the positive GPP extremes are defined as GPP anomalies greater than q and negative GPP extremes are the GPP anomalies less than −q for both simulations, with and without LULCC.
Integral of negative (or positive) GPP extremes over land grid cells represents the total global loss (or gain) in carbon uptake per month. The time series of frequency and extent (area) under GPP extremes was computed by integrating the count and area of grid cells under GPP extremes.

Attribution to Climate Drivers
Human activities, through fossil fuel emissions and land cover changes, modify the climate and climate-carbon feedbacks. To attribute significant changes in carbon uptake and GPP to climate drivers, we computed linear regression of temporally continuous GPP extremes with anomalies in atmospheric precipitation (Prcp, composed of atmospheric rain and snow), precipitation minus evapotranspiration (P − ET), soil moisture (Soilmoist, up to 1-m depth), monthly maximum daily temperature (T max ), monthly minimum daily temperature (T min ), monthly mean daily temperature (T sa ), and fire (Fire, total column level carbon loss due to fire) ( Table 1).
The impact of climate drivers on GPP often has a lagged response because the terrestrial ecosystem has ingrained plasticity to buffer and push back effects of climate change (Zhang et al., 2014). The controls of different climate drivers on GPP and its extremes is also dependent on location, timing, soil type and moisture, and vulnerability of land cover type (Frank et al., 2015). Therefore, at every location of GPP TCEs, we used linear regression to compute the correlation of TCEs with the cumulative lagged response of anomalies of every climate driver at time-lags of 0-12 months, for every 25 years time window from 1850 to 2300. However, lagged responses beyond four months were insignificant and thus we only report on one, two, and three lags. Anomalies from past months were included in the computation of lagged correlations. For instance, to compute a climate-carbon lag response of 3 months, the GPP extreme anomalies at month t were correlated with the average of climate driver anomalies at t − 3, t − 2, and t − 1 month. We also investigated the GPP TCEs driven by triggers (climate drivers during onset of TCEs) and persistent climate drivers (considering entire duration of TCEs). We identified the most dominant climate driver at any grid cell under simulations with and without LULCC, as the climate driver with the highest absolute correlation coefficient (significance value, p < 0.05). The percent global distributions of dominant drivers for every time window were computed to inspect the changing patterns of dominant drivers at various lags and over time.
A GPP TCE event could be driven by one or a combination of climate drivers. The study of compound events consisting of multiple different climate drivers leading to a GPP extreme improves our understanding of interactive effects of climate drivers (Zscheischler et al., 2018). While the climate extremes and carbon extremes often do not occur concurrently, the compound effects of climate drivers that are not climate extremes can nevertheless impact the carbon cycle (Pan et al., 2019). We studied the compound effects of water (Prcp, P − ET, and Soilmoist)-, temperature (T max , T sa , and T min )-and fire (Fire)-driven TCEs in GPP.

Detection and Identification of GPP Extremes
The 1st and 99th percentiles of the PDF of GPP anomalies for all land grid cells in 25-year time windows were used to compute the threshold trajectories  forwith and without LULCC simulations. Figure 3 shows trajectories of positive and negative GPP extremes for simulations with and without LULCC. The rising atmospheric CO 2 concentrations and increasing CO 2 fertilization, water use efficiency and lengthening of growing seasons (G. Bonan, 2015; P. J. Lawrence et al., 2012) lead to an increase in the global GPP and GPP anomalies. Consequently, Th+ and Th− increased in both simu lations; for the simulation without LULCC Th+ increased from 420 GgC/month during 1850-1874 to 664 GgC/month during 2275-2299; the corresponding values for Th+ with LULCC were 426 and 739 GgC/month, for Th− without LULCC were −465 and −755 GgC/month, and for Th− with LULCC were −469 and −813 GgC/month. The increasing magnitude of thresholds for GPP extremes over time highlights the intensification of GPP extremes over time. For the same percentile, the magnitude of negative thresholds were larger compared to positive thresholds. Hence, negative anomalies in GPP or losses in carbon uptake were much larger than gains in carbon uptake, primarily due to negative impact of increasing anthropogenic CO 2 and LULCC on carbon cycle anomalies.
Higher thresholds for the simulation with LULCC despite lower GPP (Figure 1c) compared to the simulation without LULCC, highlights that increasing magnitude of GPP anomalies were likely due to LULCC and wood harvest. Higher anomalies GPP is potentially due to large reductions in the area of tree PFTs of primary vegetation (forests) in the RCP 8.5 LULCC scenario, −3.5 × 10 6 km 2 from 1850 to 2100. The decrease in primary vegetation is associated with large increases in crop and grass areas, leading to a reduction of ecosystem carbon of −49 PgC from 1850 to 2100. The global wood harvest carbon flux increased to 4.2 PgC year −1 (Lawrence et al., 2012) in the year 2100, and then was kept at a constant harvest rate from 2100 to 2300. The legacy effects of human land cover change and continued wood harvest becomes more visible beyond 2100 when the difference in annual GPP of both simulations widens ( Figure 1c). As a result of enhanced variability in GPP with LULCC, we saw a larger magnitude of negative and positive thresholds for the simulation with LULCC ( Figure 3).
The global time series of the intensity of losses and gains in carbon uptake, calculated by integrating negative and positive extremes in GPP, are shown in Figures 4a and 4b for simulations without and with LULCC, respectively. Compared to the simulation without LULCC, the total additional loss of global carbon uptake due to LULCC was −46.53 PgC for period 1850-2100 and −141.76 PgC for 2101-2300. The respective difference in total global GPP (with LULCC minus without LULCC) was −676 PgC for 1850-2100 and −1416 PgC for 2101-2300; and relative to the total global GPP, the additional losses in total carbon uptake due to LULCC increased from 6.9% (1850-2100) to 10% (2100-2300). Hence, LULCC impacts global carbon cycle by reducing the total global GPP and increasing losses in carbon uptake during GPP extremes. The rates of increase in the intensity of positive extremes in GPP for the simulation without LULCC were 529 and 377 MgC/month for the periods 1850-2100 and 2101-2299, respectively, and for the simulation with LULCC were 863 and 692 MgC/month, respectively. The corresponding rates of increase in the intensity of negative extremes in GPP for the simulation without LULCC were −647 and −680 MgC/month, respectively, and for the simulation with LULCC were −1092 and −866 MgC/month, respectively. The changes in the rates of the intensity of positive extremes in GPP for the simulation with LULCC were analogous to the simulation without LULCC. However, the intensity of negative extremes in GPP for the simulation with LULCC shows a decrease beyond 2100, possibly due to non-increment of wood harvest rate and a constant PFT distribution at the year 2100 values for the period from 2100 through 2300 which decreases the relative variability in GPP after 2100. The larger intensity of negative extremes in GPP for the simulation with LULCC than without LULCC could result from increased wood harvest and land conversions to agriculture, pastures, and urban areas. The rate of increase in the intensity of positive extremes in GPP was higher in the simulation with LULCC than without LULCC probably due to large re-growth of secondary forests in the regions of the Eastern U.S., Europe, Africa, and South America (Hurtt et al., 2011). In the simulation with LULCC, by the year 2100 under RCP 8.5, the LULCC transitions resulted in high-density croplands in the U.S., Europe, and South East Asia; high-density pastures in the Western U.S., Eurasia, South Africa, and Australia (Hurtt et al., 2011). Primary forests were present in high northern latitudes and parts of Amazonia. Hence, LULCC forcing coupled with CO 2 forcing under RCP 8.5 and ECP 8.5 intensified both losses and gains in carbon uptake during extremes in GPP, with net losses in carbon uptake dominating the net climate-carbon response.
The rate of increase of intensities of positive and negative extremes in GPP was stronger for the simulation with LULCC, although the total GPP was less than the simulation without LULCC (Figure 1c). Stronger intensities of GPP extremes with LULCC were due to the larger IAV in GPP with LULCC compared to without LULCC ( Figure 5). LULCC alters the climate by modifying the biogeochemical and biogeophysical processes, which further affects the carbon cycle. For example, in semiarid climates, loss of vegetation cover increases the surface albedo, which increases the reflected solar radiation and cools the surface climate that weakens the boundary layer, reducing the probability of precipitation that creates the dry climate and reduces plant productivity. Biogeochemical processes include uptake of carbon during photosynthesis at an increased atmospheric concentration of CO 2 and loss of carbon during respiration in a warmer climate. Biogeophysical and biogeochemical processes do not occur in isolation and depend on the hydrologic cycle (Bonan, 2015). With increases in atmospheric CO 2 concentration, the climate becomes warmer and increases the intensity of precipitation and accompanying precipitation extremes (Ban et al., 2015;O'Gorman, 2015). Clearing of land and deforestation has led to and will lead to cooling in high latitudes, warming in tropics, and uncertain changes in mid-latitudes (Lawrence et al., 2016). As a result, human LULCC increases the regional heterogeneity in vegetation that alters the climate drivers (Ichii et al., 2005), further increasing the global interannual variability and anomalies of GPP as well as the spatially heterogeneous distribution of GPP extremes. Therefore, the effect of both CO 2 and LULCC forcing, represented in with LULCC, increases the variability of biogeochemical and biogeophysical feedbacks, thus resulting in increased IAV of GPP in with LULCC than without LULCC.
The regional atmospheric circulation and climate change lead to spatial variations in the distribution of GPP extremes and the intensities of GPP extremes. Since the definition of GPP extremes is based on the threshold trajectory of first percentile global anomalies of GPP for consecutive 25-year time windows in the simulation without LULCC (Figure 3), the total number of grid cells under positive GPP extremes were constant at around 64,000 (for every 25-year time windows) while they vary for all other scenarios. The total number of grid cells and area affected by negative extremes in GPP for the simulation with LULCC were largest, possibly because of increased negative feedback of climate variability on the carbon cycle due to the cumulative CO 2 and LULCC forcing. Relative to the frequency of positive extremes in GPP for the year 1850 (without LULCC), the frequency of positive extremes (with LULCC) increased by 17% and 28% for the periods 1850-2100 and 2101-2300 respectively; and the respective growth rates of negative extremes (with LULCC) were 13% and 19%. For with LULCC, growth rates for the area affected by positive GPP extremes were 16% and 28%, and for negative GPP extremes at 12% and 20% during 1850-2100 and 2101-2300 respectively. Higher growth rates of frequency and area-affected by positive extremes in GPP were probably due to increases in secondary forest cover; however the losses in the expected carbon uptake due to LULCC were larger than gains.

Attribution to Climate Drivers
The attribution of climate-driven extremes in GPP were performed for GPP TCE events, which are time-continuous GPP extremes meeting the criteria described in Section 3. Variability in terrestrial carbon cycle intensified the climate-driven GPP TCEs under the combined forcing of human-induced LULCC and anthropogenic CO 2 emissions. The mean duration (Figure 6a) and standard deviation ( Figure S2 in Supporting Information S1) of negative and positive TCEs in GPP for the simulation with LULCC were greater than the simulation without LULCC. In addition, the duration of TCEs lengthened over time, with more long-duration TCEs and fewer short-duration TCEs as time progressed. Figure 6b shows the density plot of the count of negative TCEs in GPP versus the total duration of GPP TCEs in 25-year windows for the simulation with LULCC. The increasing mean duration of the negative TCEs in GPP in the future likely causes a larger reduction in the carbon uptake, which has significant implications for carbon storage and the carbon budget, and it can negate the positive feedback of increasing CO 2 fertilization under RCP 8.5 and ECP 8.5 CO 2 scenarios.
Extremes in carbon cycle and climate extremes often do not occur simultaneously (Pan et al., 2019). In this study, we first detect TCEs in GPP and then attribute the climate drivers using linear regression. Pearson's correlation coefficients and corresponding significance values (p-values) were computed between every climate driver and extreme anomalies in GPP during TCEs. The dominant climate driver at every land grid cell was defined by the maximum absolute correlation coefficient of climate drivers (p < 0.05). To consider the response of the ecosystem to prevailing climatic conditions, linear regressions were performed at multiple lags between extreme anomalies in GPP and climate driver anomalies. Figure 7 (and Figure S3 in Supporting Information S1) show the percent distribution of dominant climate drivers for every 25-year time window for both simulations, with LULCC and without LULCC, respectively, at lags of 1, 2, and 3 months. The relationship of water availability (soil moisture) with TCEs in GPP was globally dominant during most time windows (Figure 7). Reduction in soil moisture leads to an anomalous reduction in photosynthesis (Frank et al., 2015). Soil moisture is the most dominant driver of negative carbon cycle extremes, highlighting a strong positive correlation of plant productivity with soil moisture, as shown in Figure S4b in Supporting Information S1. Temperature (T max and T sa ) and Fire have a negative correlation with TCEs in GPP (spatial distribution shown in Figures S4d and S4e of Supporting Information S1), where an anomalous increase in these drivers leads to loss in carbon uptake. T max was dominant at a lag of 2 months through 2100, but Fire shows dominance at higher lags, especially beyond 2100 in the simulation without LULCC. The percent count of dominant temperature drivers, T max and T sa , increased over time, especially after 2100. Hence, carbon extremes driven by hot conditions will have the largest increase in a warming climate, especially after 2100. The lagged response of anomalously hot air temperatures and lack of soil moisture creates hot and dry conditions suitable for fire events. Fire events driven by hot and dry conditions increase at a higher rate at all lags in the simulation with LULCC than the simulation without LULCC. After 2100 in the simulation with LULCC, fire stands out as the dominant driver at all lags, highlighting that the human-induced LULCC will make ecosystems more vulnerable to fire events. The results of our study are consistent with the findings of Williams et al. (2014) and Zscheischler et al. (2014), suggesting that the declines in GPP with warm temperature extremes are due to dependence of GPP on soil moisture, and to the strong negative correlation between soil moisture and temperature. The rising temperature under global warming will create a warmer environment in the future, increasing the risk associated with heatwaves and their impacts on the ecosystem. A decline in precipitation and soil moisture compounded with warm temperature may cause an unprecedented increase in loss of carbon uptake and potentially reduce the terrestrial carbon sink.
The combined effect of climate drivers often has a larger impact on extremes in carbon cycle than the simple addition of individual climate drivers (Flach et al., 2020;Frank et al., 2015;Zscheischler et al., 2014Zscheischler et al., , 2018. To capture the compound effect of climate drivers (Table 1) on GPP TCE events, the climate drivers were grouped under fire, temperature, and water-driven extremes (Section 3.2). Figures 8 and S5 show the spatial distribution of climate drivers at a lag of 1 month for simulations with and without LULCC, respectively. We observed an increase in the number of TCEs in GPP and changes in the spatial distribution of GPP TCEs and climate drivers for time periods 1900-1924, 2000-2024, 2100-2124, and 2200-2224 in the simulation with LULCC (Figure 8).
Most locations experienced TCE events due to the combined effects of all the climate drivers.
We computed all combinations of water (dry or wet), temperature (hot or cold), and fire-driven climatic conditions that cause GPP TCEs. For brevity, we report results for one 25-year time window one per century. The fractions (larger than 0.05) of total negative TCEs in GPP driven by mutually exclusive climate conditions are shown in Figure 9. More than half of the total TCEs in GPP occurred when the environmental conditions are hot and dry, resulting in fire events and account for large reductions in GPP (Figures 8 and 9). Similar to other studies (Frank et al., 2015;Zscheischler et al., 2014), a stronger correlation between warm and dry climate has a substantial impact on the terrestrial carbon cycle. The increased fraction of TCEs attributed to hot climate over time is possibly due to combined attribution of T max and T sa , as shown in Figure 7. Figure S6 of Supporting Information S1 represents the total fraction of TCEs in GPP driven by mutually inclusive climatic conditions. As global warming increases the temperature of the planet, ecosystems will become more vulnerable to the hot and dry climate and at increasing risk of fire events (Figure 9).
Here we highlight the temporal changes in the distribution of carbon cycle extremes and their climate drivers for the simulation with LULCC. East Asia experiences a large number of TCEs during the 21st century (Figure 8) driven by dry climate and fire, which however declines in the 22nd century due to an increase in precipitation and available soil moisture. The Savannas near Congo in Africa display an increased vulnerability of vegetation to fire and hot temperatures. The impact of LULCC in Savannas is prominent after 2100 ( Figures S5d and 8d in Supporting Information S1) when the pattern of dominant drivers change from temperature-driven (green color) to a compound effect of increase in dry, hot, and fire conditions (gray color). The contiguous United States (CONUS) experiences more GPP TCEs, especially in the 23rd century (Figure 8d), driven by hot temperature and water limitation (highlighted in cyan color). Indonesia and its neighboring islands show increasing GPP   The fractions are mutually exclusive, that is, events driven by hot and dry climate is not counted in either hot or dry climate driven events. Any location could be affected by one or compound climatic conditions. For example, a carbon cycle extreme could be driven by any combination of hot or cold, dry or wet, and with or without fire. We only show the combination of driving climate drivers that have total fraction of more than 0.05. The combined effect of hot and dry climate accompanied by fire leads to most negative TCE events in GPP.
TCEs attributed to a combination of hot temperatures and fire. The number of negative TCEs in South America are approximately the same for periods in Figures 8a-8c except for the period 2200-2224 (Figure 8d) when the extremes show a large increase. During 2200-2224, South America shows a drastic increase in negative TCEs i.e., attributed primarily to dry, hot, and fire conditions (represented by gray color). Africa experiences an increase in the number of GPP TCEs attributed to hot climate, and, in the far future (2101-2300), witnesses an increased frequency of fire events. Australia's east coast also exhibits an increase in GPP TCEs because of water limitation, dry climate, and they are accompanied by fires in the far future. The South and Southeast Asia (composed of India, Myanmar, Thailand, and Cambodia) experiences an increase in TCEs in GPP over time.
These Asian regions where primary forest were converted to cropland (Hurtt et al., 2011), saw an increase in fire-driven extreme events in GPP, highlighting a rising vulnerability to fire events due to LULCC.

Regional Analysis of Climate Change Impact on Climate Drivers and Negative Carbon Cycle Extremes
Most climate models show an increase in the interannual variability of land-atmosphere CO 2 exchange over time (Keenan et al., 2012). The increasing atmospheric CO 2 concentration influences climate through its radiative effect (i.e., Greenhouse effect) and indirectly through physiological effect (reduced plant transpiration) (Cao et al., 2010). The increasing radiative effects cause changes in circulation, impacting precipitation, soil moisture, and global scale water cycle. Increase in precipitation and CO 2 fertilization and reduction in stomatal conductance could lead to an increase in GPP; however, droughts, heatwaves and fires curtail plant photosynthesis (Swann et al., 2016). Under climate change, while most regions experience an increase in mean GPP, some regions exhibit a decline.
Reduction in growth rate of expected terrestrial carbon uptake below carbon emissions could lead to weakening of carbon sink capacity ( Figure S7 in Supporting Information S1). Increase in total GPP and IAV in GPP over time contribute to the increase in magnitude of negative carbon cycle extremes, defined here as negative extremes in GPP. However, some regions experience weakening of negative carbon cycle extremes often due to decline in total GPP, and reduction in IAV of GPP and climate drivers ( Figure S8 in Supporting Information S1). Driven by large scale circulation changes, increasing temperature and decreasing precipitation and soil moisture in Central America (CAM) lead to forest mortality, decline in GPP and weakening of negative carbon cycle extremes. The decrease in precipitation in the Northern South American Tropics (NSA) and Central Amazon Basin (CAB) and increase in precipitation in the Southwest Amazon Basin (SAB) is due to plant physiologic response. Langenbrunner et al. (2019) investigated this phenomenon in the regions of the Amazon and the Andes and found that under increased atmospheric CO 2 , stomatal conductance decreases, reducing water loss through transpiration, which decreases the convective activity and causes a reduction in rainfall over the Amazon Basin (CAB, NSA) and increases moisture advection by low-level westward jets toward the Andes, leading to increases in the rainfall over the Andes (SAB). Reducing variability and increasing precipitation over SAB increase available soil moisture increasing regional GPP, however with reduced IAV, resulting in weakening of negative carbon cycle extremes. Although CAB witnesses a decline in precipitation it demonstrates an increase in GPP likely due to slight reduction in soil moisture and large CO 2 fertilization effect which compensates for the reduction in soil moisture. Hence, in CAB both GPP and IAV in GPP increases over time causing strengthening of negative carbon cycle extremes. These regional patterns of changes in GPP, negative carbon cycle extremes, and climate drivers are seen across the globe. For example, the regional changes in Northern Guinea, Southern Guinea, Central Guinea and Indonesia resemble the processes in CAM, SAB, CAB and CAM, respectively. Additional feedbacks due to LULCC -as tropical forests are replaced with crops or grasslands that transpire less than forests, further reducing rainfall over the Amazon -causes hotter and drier climate with an increasing risk of fire and larger losses of carbon uptake (Figures S9 and S10 in Supporting Information S1).
Our robust attribution analysis based on successive time windows captures the changes in the anomalies of climate drivers and carbon cycle. During 1900During -1924During and 2000During -2024, most of the negative TCEs in GPP that occurred in SAB (Figure 8), were due to water limitation (highlighted by the blue color in the RGB maps). With increased water availability in SAB (Langenbrunner et al., 2019), the effects of heat associated with temperature increases are mitigated and the number of negative carbon cycle extremes reduces. The decrease in precipitation in NSA led to a reduction in soil moisture, resulting in a hot and dry climate, thus making the vegetation prone to fire at high temperatures (the compound effect is highlighted by the gray color). The reduced GPP in NSA over time further leads to a fewer number of negative carbon cycle extremes (Figure 8). Large IAV and anomalies in GPP in CAB lead to increased concentration of negative TCEs in GPP driven possibly due to combined effect of increased hot, dry, and fire events (shown in gray color).

Drivers and Triggers of Carbon Cycle Extremes
A grid cell could experience any number of extreme events of any length during successive 25 years time windows from 1850 through 2300. Our methodology helps identify long temporally continuous extremes which represent large magnitude carbon extremes and have higher significant regression coefficient. We use this methodology to identify the prevailing climatic conditions that act as triggers for an extreme and the conditions that cause an extreme to persist. We demonstrate the findings from Chaco Province in Argentina where the major plant functional type is broadleaf deciduous tree (43%).
Under the simulation without LULCC and the time period 2000-2024, the total number of positive TCE events in GPP were seven with a total duration of 53 months, which were greater than the total five negative TCEs in GPP with a total duration of 40 months (Table S2 in Supporting Information S1). The total gain in carbon uptake was 35.17 TgC, which was greater than the total loss in carbon uptake of −28.15 TgC, resulting in a net gain in carbon during TCE events. Under the simulation with LULCC, during the same period of 2000-2024 the total number of negative GPP TCEs (6 events, 57 months) were greater than positive TCEs (5 events, 35 months). Thus, the total loss in carbon uptake (−35.86 TgC) with LULCC was higher than the gain in carbon uptake (23.87 TgC). While this region had a net gain in carbon uptake of 7.02 TgC for the simulation without LULCC, and it had a net loss in carbon uptake of −11.99 TgC for the simulation with LULCC. The increase in negative extremes in GPP with net losses in carbon uptake demonstrates the role of human-induced LULCC in negatively affecting carbon uptake capacity.
We performed a qualitative investigation of every negative and positive TCEs in GPP ( Figure S11 in Supporting Information S1) using normalized time series of GPP, GPP anomalies and anomalies of major climate drivers (Prcp, Soilmoist, T max , and Fire). Most positive TCEs in GPP were driven by increase in precipitation, followed by rise in soil moisture and decline in T max . Most negative extremes in GPP were driven by decline in precipitation and rising T max , which cause high evapotranspiration and loss of soil moisture, eventually creating dry and hot conditions that cause fire events. Some of the negative extremes in these regions were also found to be driven by hot and dry events without fire.
To identify the triggers of carbon cycle extremes, we define the onset period of GPP TCEs as first-quarter of every TCE. The regression onset period of GPP TCEs and anomalies of climate drivers were computed with consideration of the lagged response of climate drivers on GPP TCEs. The lagged soil moisture anomalies were highly correlated (p < 0.01) with entire duration of TCE events (persistence, Table S3 in Supporting Information S1), and the lagged precipitation anomalies were highly correlated with the onset (trigger, Table S4 in Supporting Information S1) of TCE events in the Chaco Province.
For the simulation with CO 2 only forcing (i.e., without LULCC) during time period 2000-2024, the correlation coefficient, for the whole duration of GPP TCEs, was highest for soil moisture anomalies at a lag of 1 month. The precipitation anomalies had high positive correlations with GPP TCEs (i.e., reduction in precipitation was correlated with a decrease in GPP) at lags of 3 and 4 months. T max and Fire had large negative correlations with GPP anomalies (i.e., increase in temperature correlates with a decrease in GPP or increase in carbon uptake loss) at a lag of 2 months. Hence, for a negative GPP TCE event, the compound effect of a decline in precipitation followed by increase in temperature caused a reduction of soil moisture, resulting in hot and dry conditions, increasing the probability of occurrence of fire and causing a long negative TCE event (Prcp > T max > Soilmoist > TCE). For the same region under simulation with LULCC, the dominant trigger was T max (Table S5 in Supporting Information S1) followed by precipitation. Therefore, human-induced LULCC coupled with increasing CO 2 levels tends to enhance the vulnerability and loss in vegetation due to hot and non-dry climate conditions, resulting in more negative extremes in GPP. The analysis illustrates the strength of our methodology in identifying the evolution of carbon cycle extremes for multiple conditions and at fine resolution.

Limitations of CESM1(BGC)
The ILAMB benchmarking scores (Collier et al., 2018) of the carbon fluxes and climate drivers indicate that CESM1(BGC) is among the best CMIP5 models (Table S1 in Supporting Information S1). However, the Community Land Model version 4 (CLM4), the land model used in the CESM1(BGC) (Oleson et al., 2010), had some limitations which could potentially impact some of our findings. The simulated GPP by CLM4 had a positive bias across the globe compared to FLUXNET eddy covariance tower estimates due to lack of colimitation of GPP to canopy scaling and parameterization of leaf photosynthesis kinetics (G. B. Bonan et al., 2011). With improvements in the model parameterizations of radiative transfer, leaf photosynthesis and stomatal conductance, canopy scaling of leaf processes, inclusion of multilayer canopy model, and updated maximum rate of electron transfer parameters (G. B. Bonan et al., 2011;, the positive bias in GPP was reduced. With inclusion of vertically resolved CN model in CLM4.5 (Koven et al., 2013), the positive bias in GPP was further improved and terrestrial carbon storage increased consistently with observations. In the current study focused on the patterns of carbon cycle extremes, the positive bias of GPP was likely captured by trend and removed for calculation of GPP anomalies. However, any associated increase in the IAV of GPP corresponding to positive bias of GPP would increase the magnitude of GPP anomalies and both negative and positive GPP extremes. Thus, there is a potential to overestimate the magnitude of both negative and positive GPP extremes but the relative comparison is insightful and most likely not affected by this limitation of CLM4.
CLM4 lacks representation of dynamic crop, thus the cooling effects in irrigated lands, changes in the sensible and latent heat, and soil carbon change are not well represented in CLM4.0 but were improved in the CLM5.0 (Lombardozzi et al., 2020). The simulated climate-carbon feedbacks in CLM4 assume an instantaneous response of plant physiological processes to changes in temperature (Lombardozzi et al., 2015). With inclusion of temperature acclimation, though the ecosystem carbon storage pool grew (mostly in higher latitudes), the effect on photosynthesis of tropical regions was minimal. Since most of the detected extremes in the tropical forest and other high biomass regions, the lack of representation of agriculture and temperature acclimation on GPP are not expected to have significant effect on the analysis of extremes.

Conclusions
Using the fully-coupled Earth system model, CESM1(BGC), we analyzed the development of extreme events in GPP and attributed those carbon cycle extremes to climate drivers from the year 1850 through 2300 for simulations with and without LULCC. The changes in land cover directly modify the biogeophysical and biogeochemical feedbacks of terrestrial vegetation and indirectly through the physiological responses of vegetation on climate drivers, leading to increase in interannual variability in GPP and higher magnitude of GPP anomalies. Relative to the simulation without LULCC, the simulation with LULCC exhibited increased variability in GPP and higher intensity, duration, extent, and frequency of extremes in GPP. These characteristics were greater for negative extremes in GPP than positive extremes, implying larger than expected losses in carbon uptake than carbon gains. Although the total GPP for the simulation with LULCC was less than the total GPP without LULCC, the simulation with LULCC showed a sharper decline of carbon uptake of −6.9% in the near future (1850-2100) and −10% in the far future (2100-2300) with respect to total GPP without LULCC.
Increasing atmospheric CO 2 concentrations drive growth in vegetation photosynthesis or GPP due to carbon fertilization, and reduction in stomatal conductance. Due to circulation changes, a few regions witness a decrease in precipitation, creating a drier climate that when supplemented with warmer temperatures lead to decline in GPP, such as the regions of the eastern Amazon and CAM. Reducing magnitudes of GPP over time and analogous decreases in the interannual variability of GPP produce weakening of negative carbon cycle or GPP extremes in the Amazon. Moreover, the weakening of other negative carbon cycle extremes could be a result of less variable and benign climatic conditions for vegetation productivity, as seen around the Andes. Hence, it is imperative to inspect trends in GPP and negative carbon cycle extremes simultaneously to understand the nature of changes in extremes and their implications at regional scales.
We found that the duration of GPP TCEs increased with higher CO 2 concentration. The duration and impact of negative TCEs in GPP are enhanced when human-induced LULCC was considered, resulting in increased loss in carbon uptake. We illustrated this with the case study of the Chaco Province in Argentina that had a net gain in carbon uptake (7.02 TgC) during 2000-2024 in the simulation without LULCC and while a net loss in carbon uptake  in the simulation with LULCC.
The dominant climate driver was soil moisture that had highest correlations with GPP extremes (p < 0.05) and the correlations were mostly positive, indicating that anomalous decrease in soil moisture or drier climate conditions cause anomalous reductions in GPP. Other major individual drivers were hot temperatures and fire. Fire was a dominant climate driver, especially after 2100 in the simulation with LULCC, highlighting the increased vulnerability of ecosystems to fire events due to the impact of human activities on ecosystems. We also found that the decline in precipitation triggers a negative carbon cycle TCE event, but the reduction in soil moisture or water limitation was the dominant driver for those negative carbon cycle TCE events to persist. The compound effects of climate drivers were also analyzed, and an increasing number of regions under carbon cycle extremes were attributed to hot and dry conditions. The largest fraction of negative carbon cycle extremes were driven by the compound effects of hot, dry, and fire events. Warmer conditions under climate change increases the risk of occurrence of fire events and their impacts on the carbon cycle and extremes of the future.
This study presents a detailed analysis of detection and identification of carbon cycle extreme events, how these extremes evolve from 1850 through 2300, and how human-induced LULCC alters them using a fully coupled Earth system model forced with the RCP 8.5 and ECP 8.5 CO 2 concentration scenarios. It also attributes the climate drivers of such extremes in carbon cycle, and analyzes the changing patterns and dominance of climate drivers, under with and without LULCC scenarios. Study provides new insights into the contribution of human activities in altering carbon cycle extremes and the vulnerability of terrestrial vegetation and associated ecosystem services that could present increasing risks to human lives, wildlife, and food security.

Disclaimer
The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy. gov/downloads/doe-public-access-plan).

Data Availability Statement
The CESM1(BGC) model output used for detection and attribution of carbon cycle extremes in the study are available at https://doi.org/10.5281/zenodo.5548153. Data analysis was performed in Python, and all analysis codes are publicly available on GitHub at https://github.com/sharma-bharat/Codes_Carbon_Extremes_2300 and archieved at https://doi.org/10.5281/zenodo.6147120.