Process-based analysis of terrestrial carbon flux predictability

Despite efforts to decrease the discrepancy between simulated and observed terrestrial carbon fluxes, the uncertainty in trends and patterns of the land carbon fluxes remains high. This difficulty raises the question to what extent the terrestrial carbon cycle is predictable, and which processes explain the predictability. Here, the perfect model approach is used to assess the potential predictability of net primary production (NPPpred) and heterotrophic respiration (Rhpred) by using ensemble simulations conducted with the Max-Planck-Institute Earth System Model. In order to assess the role of local carbon flux 5 predictability (CFpred) on the predictability of the global carbon cycle, we suggest a new predictability metric weighted by the amplitude of the flux anomalies. Regression analysis is used to determine the contribution of the predictability of different environmental drivers to NPPpred and Rhpred (soil moisture, air temperature and radiation for NPP and soil organic carbon, air temperature and precipitation for Rh). Global NPPpred is driven to 62 and 30% by the predictability of soil moisture and temperature, respectively. Global Rhpred is driven to 52 and 27% by the predictability of soil organic carbon and temperature, 10 respectively. The decomposition of predictability shows that the relatively high Rhpred compared to NPPpred is due to the generally high predictability of soil organic carbon. The seasonality in NPPpred and Rhpred patterns can be explained by the change in limiting factors over the wet and dry months. Consequently, CFpred is controlled by the predictability of the currently limiting environmental factor. Differences in CFpred between ensemble simulations can be attributed to the occurrence of wet and dry years, which influences the predictability of soil moisture and temperature. This variability of predictability is caused 15 by the state dependency of ecosystem processes. Our results reveal the crucial regions and ecosystem processes to be considered when initializing a carbon prediction system.

limited understanding of carbon flux variability, models are not able to fully reproduce the spatiotemporal patterns of the ter-25 restrial carbon cycle. This is reflected in the poor representation of soil organic carbon (SOC) in Earth System Models (ESM), the inability to adequately model gross primary production (GPP) from eddy covariance flux tower sites (Luo et al., 2015) and the difficulty to detect the efforts taken in emission reduction due to internal variability of atmospheric CO 2 variability (Spring et al., 2020). In order to produce more realistic predictions, efforts in model development have been directed towards using observations to constrain model parameters (Zeng et al., 2014;Bloom et al., 2016;Mystakidis et al., 2016;Chadburn et al., 30 2017; Tziolas et al., 2020), and to refine model structure to incorporating more processes and interactions (Krull et al., 2003;Stockmann et al., 2013;Xu et al., 2014;Luo et al., 2016). While efforts in model development are continuing to narrow the gap between the simulated and observed carbon cycle, the lack of progress in improving the predictive ability of the models raises the question to what extent the terrestrial carbon cycle is predictable at all (Luo et al., 2015).
The potential predictability of a system can be estimated by using the perfect model framework. Ensemble simulations 35 are initialized along a control run with each member of the ensemble having slightly perturbed initial conditions. The upper limits of predictability are then derived by analysing the divergence of the ensemble simulations. This method assumes (a) perfect model physics which are able to reproduce the full spectrum of natural variability and (b) perfect knowledge of the modelled system and a model whose representation of the real world is "perfect enough" (Boer et al., 2013). Séférian et al. (2018) used the perfect model framework to assess the potential predictability of terrestrial carbon fluxes (CFpred) at annual 40 time-steps. They estimated the predictive horizon of terrestrial carbon fluxes to be two years globally and up to three years in northern latitudes. The high variability of predictability among different initializations suggests a state-dependence of CFpred, but no further mechanisms of predictability were investigated therein. Multiple processes can being considered as the sources of CFpred. Due to the high sensitivity of the terrestrial carbon cycle to climate, climate predictability provides carbon fluxes with a basic prediction horizon. The main contributor to climate predictability is El Niño-Southern Oscillation (ENSO), which Previous studies that focus on the mechanisms of CFpred investigated the role of various land processes and how they contribute to the overall CFpred. Weiss et al. (2014) found increased predictability of evaporation and to some extent temperature 60 due to a dynamic simulation of leaf area index (LAI), which would also extend CFpred. The role of land surface initialization on CFpred was studied by Zeng et al. (2008) and Lovenduski et al. (2019). Zeng et al. (2008) isolated the fraction of CFpred which is based solely on initial conditions and compare fully coupled dynamic simulations with statistical models. Lovenduski et al. (2019) quantified the degree to which CFpred improves when the land surface is initialized. They also assessed the relative importance of the individual land surface processes for the variability of terrestrial carbon fluxes and found that CFpred 65 depends on the correct initialization of vegetation carbon biomass and soil moisture rather than temperature. These studies have shown the significant advantage of dynamic forecasting systems, suggesting CFpred extends beyond the predictability of the forcing variables due to land surface processes. However, these studies were not focused on contributions of individual drivers of carbon fluxes to CFpred or on processes responsible for maintaining CFpred.
Here, we use perfect model simulations conducted with an ESM to investigate the structure and mechanisms of the CFpred.

70
Initialized ensemble simulations are created from a range of ENSO states. Analysed are the carbon fluxes with the highest contribution to the interannual variability of the land-atmosphere CO 2 exchange. These are NPP with an interannual SD of 0.99 PgCyr −1 and heterotrophic respiration (Rh) with an SD of 0.29 PgCyr −1 (Wang et al., 2016). The potential predictability of NPP (NPPpred) and Rh (Rhpred) is derived from rate of divergence within the ensemble members. We evaluate the predictability data to find how NPPpred and Rhpred differ in their spatiotemporal patterns and variability. Lastly, we identify the 75 key drivers of NPP and Rh and determine their contribution to NPPpred and Rhpred. We use this framework to explain the attained spatiotemporal patterns of CFpred and identify the underlying land system processes producing these patterns.
The potential predictability is assessed by using a correlation-based and a distance-based metric. The Anomaly Correlation 90 Coefficient (ACC) is a commonly used metric to measure forecast skill (Jolliffe and Stephenson, 2012) which calculates the correlation between predicted and observed anomalies as where j and t are grid cell and lead time, cov is the covariance and f and o the forecast and validation anomalies. Similar to (Collins and Sinha, 2003;Becker et al., 2013), the noise in the ACC is reduced by averaging over several ACC values. This is 95 achieved by taking all eleven ensemble members as the validation in turn, while the mean of the remaining ensemble members serves as the forecast. Although the ACC is an intuitive metric which is calculated from all initializations and thus provides a robust estimation of the predictability, it does not allow to investigate the variability of predictability between initializations.
The comparison of predictabilities between initialization is achieved by the use of a distance based metric which is computed for all initializations individually. The distance based metric used here is the normalized ensemble variance (V (t)) based on the 100 method proposed by Griffies and Bryan (1997). Predictability is defined as the ensemble variance normalized by the variance of the climatology as where t is lead time, M the number of ensemble members, X i the ith member, X the ensemble mean and σ 2 the variance of the control simulation. In this study, the complement of the normalized ensemble variance is used as: resulting metric indicates perfect predictability at a value of 1, and an ensemble spread that exceeds the climatological variance for values below zero.
While ACC and V c allow the estimation of regional predictability, these metrics are not suitable to evaluate the impact of local predictabilities on the predictability of the global carbon cycle. This is due to the disregard of the flux amplitude in the calculation of the metrics. Both of the metrics are prone to producing above average predictabilities in regions where carbon 110 fluxes are generally low or even close to zero, such as subtropical deserts. Here we propose a weighted predictability metric that allows to assess local predictabilities with regard of their impact on the predictability of the global carbon cycle. V c is weighted by using an approach similar to risk assessment, which is calculated as the product of likelihood and impact. Here a weighted predictability wV c is calculated by multiplying V c with the absolute carbon flux anomaly of the ensemble mean:

Decomposition of predictability
In order to investigate the drivers of CFpred, the V c of NPP and Rh are decomposed into components contributing to the predictability of these fluxes similar to an approach used by Jung et al.
with V c F LU X being the complementary normalized ensemble variance of NPP or Rh, a DRI k coefficient of the kth driver (for example TEMPpred), V c DRI the predictability of the kth driver and the residual error term. Grid cell, lead time and initialization are denoted by the indices j, t and i. The regression coefficients are calculated by using non-negative least squares (Mullen and van Stokkum, 2012) for every grid cell and lead time by using the data from all initializations. After fitting the 135 regression model to the data, the individual components of CFpred are calculated as where V c F LU X DRI describes the amount of predictability of F LU X that can be contributed to the driver k.

Results and discussion
Out of the 35 ensemble simulations initialized along the control run, 7 simulations are part of the El Niño and 8 simulations 140 part of the La Niña group (Fig. 1). The El Niño simulations peak between the September before initialization and January with peak values between 2.2 and 3.6 • C (three month running mean Niño 3.4 SST anomaly). They show a fast decline of the anomaly with most models having a negative anomaly in December of the first year and evolving into a La Niña event in the second year. Peaks of the La Niña simulations fall between September and June and, while their relative peak anomalies are smaller (-1.6 to 3.0 • C), the negative anomaly can be sustained well into the second year. aspects. While the ACC of NPP has a slower temporal decline with values above 0.8 for 2 to 3 months around the equator, the ACC of Rh drops below 0.5 within the first two months for most latitudes. However, Rh shows much higher long term predictability, especially in the second year of the simulation where Rhpred is much higher than NPPpred.
While both predictability patterns show signs of a seasonal cycle, they are out of phase with Rhpred distinctly following the wet season and NPPpred appearing to be higher in the dry seasons of the first year. This has a large role on the comparability of  These findings highlight the importance of correctly simulating the ESNO process. Especially the localization of ENSO related rainfall patterns is crucial, since they provide a sustained and predictable anomaly in water availability.
Many of the identified spatial patterns of CFpred can be discovered in similar studies. Most models agree on the Amazon 175 basin as the global hotspot of CFpred (Zeng et al., 2008;Ilyina et al., 2021), and some reflect the increased predictability in Southeast Asia and southern Africa (Zeng et al., 2008), but the comparison of predictability horizons remain difficult due to the use of different predictability metrics.
The results reveal different areas in which an operational NPP forecast can be used to increase food security. The high NPPpred of the Sahel and Kalahari savanna ecosystems (Fig. 3) could be used to plan stocking rates in order to avoid grassland 180 Figure 3. ACC of NPP and Rh. The color scale is cropped at zero. Only values above the 95% confidence interval are shown.
degradation due to overgrazing in dry years (Tews et al., 2006). Other promising regions are Northeast and central Brazil. The high NPPpred in these areas could be used to select of crop varieties which are more or less drought tolerant depending on the given forecast.  by SOILpred has a spatial extent that broadly covers all regions of high NPPpred, TEMPpred is concentrated to certain areas.

Composition of predictability
TEMPpred is high in a band spanning from the Amazon basin to northern South America, in southern Africa and in Southeast Asia. The largest contributor to Rhpred is SOCpred (52%) followed by TEMPpred (27%). Similar to NPPpred, the temperature In order to facilitate a system for operational NPP prediction, a network of sensors could be installed to gather data on the initial condition of the land surface. The patterns of the role of soil moisture in predicting NPP (Fig. 4) reveal the areas on which the efforts in establishing such a network should be focused on to maximize the impact.
There are more variables that are regarded as key drivers of NPP variability and could have been considered as predictors 200 in the regression models. Most importantly, LAI and humidity play an important role in NPP variability (Schaefer et al., 2002). Several studies show the role of a dynamical simulation of LAI in extending the predictability of land surface processes (Zeng et al., 1999;Wang et al., 2010Wang et al., , 2011Weiss et al., 2012Weiss et al., , 2014. Here, the inclusion of LAI as a predictor is rejected because of the susceptibility of regression models to correlated predictors. The changing concentration of atmospheric CO 2 is causing trends in NPP as global atmospheric levels are rising (Winkler et al., 2021), however, we assumed that the interannual 205 variability of CO 2 fertilization is below a meaningful contribution to overall variability. Although clay content plays a major role on carbon turnover rates in soil (Coleman et al., 1997), it is not considered in the JSBACH Rh submodel (Tuomi et al., 2009) and was not included in this study.

Seasonality
The seasonal patterns of NPPpred revealed in the ACC data ( Fig. 2 and 3) are reproducible by the decomposed predictability 210 metric V c (Fig. 5). They show the reemergence of predictability in the dry season at various locations, and reveal that this phenomenon can not be contributed to a single factor. The largest pattern is a reemergence in July to November at 1 • to 4 • South and can be associated with the high NPPpred in the southern Amazon ( Fig. 3 NPP September 1st year). This pattern is due to increased TEMPpred throughout the dry season, which is extended by high deepSOILpred in September, and even reoccurs in the second year of the simulation. Another pattern explains the high NPPpred in southern Africa between August 215 and October, which is due to deepSOILpred.
These cases of high dry season NPPpred in the tropics are most likely due to the seasonally changing limitations of NPP.
During the productive wet season, plant growth is limited by incoming radiation (Wang et al., 2010) which has little variability and poor predictability. Instead, most of the interannual variability of NPP can be explained by dry season variability. One study found over 80% of western Amazon NPP variability to take place between July and September (Wang et al., 2011). The 220 water limitation of NPP during the dry season (Tian et al., 2000) does not only introduce higher variability as compared with the energy limited wet season, but the coupling of NPP to soil moisture also lends NPP the high predictability of soil moisture.
Although the seasonality of Rhpred shows a reversed tendency to NPPpred with higher predictability in the wet season, the mechanisms explaining the seasonality are similar. The seasonally varying Rhpred can be explained by the inherently different predictability of the seasonally dominant limiting factor of Rh (Fig. 6). During the dry season the limiting factor of 225 Rh is precipitation, which has a generally low predictability. The absence of precipitation for several weeks will inhibit soil respiration completely. There is a sharp increase in Rh variability in the dry-wet transition because the onset of precipitation is difficult to predict. As precipitation increases, the moisture constraint is asymptotically lifted and approaches zero. At this point, Rh becomes limited by substrate availability, which has a much higher predictability than climatic variables. The high SOCpred is due to the persistence of SOC anomalies because of the low decomposition rates and the pause of decomposition during dry seasons. Although TEMPpred is higher than PRECIPpred, it only plays a minor role in tropical Rhpred, because tropical Rh has relatively low temperature sensitivity (Meir et al., 2008).
These pronounced seasonal patterns of Rhpred hinge on the implementation of the precipitation sensibility function in MPI-ESM. The shape and parameterization of the rate modifying function of decomposition to moisture sets Rh to be more sensible to precipitation in the dry than in the wet season. However, the relationship between Rh and moisture in the tropics is the highly 235 debated subject of various studies coming to different conclusions. These studies suggest a parabolic or no relationship with soil moisture (Meir et al., 2008) or linear increase with precipitation (Tian et al., 2000).

Interannual variability
Using the distance based predictability metric V c also allows to evaluate the variability of predictability between different initializations. Among the regions with the highest interannual variability of NPPpred are the southern Amazon basin (Box in 240 Fig. 7), with a mean V c of 0.24 and an SD of 0.32, and northwestern Australia, with a mean of 0.16 and an SD of 0.60 (23 • S, 122 • W). Figure 7 shows how the interannual variability of NPPpred is affected by initial soil moisture. The majority of regions with a high NPPpred ( Fig. 3 and 4) have a higher predictability in years initiated from wet states. Exceptions to this trend are India and northwestern Australia, where NPPpred is higher in dry years. The strongest difference in NPPpred is in the Amazon basin, where overall NPPpred and interannual variability of predictability are also at the global maximum.

245
To determine the mechanisms responsible for this difference in predictability we focus on the composition of the NPPpred in the southern Amazon basin (box in Fig. 7). To represent wet and dry years, a composite analysis is used based on the ENSO states (The El Niño years are the driest extremes at initialization while soils are often saturated at the beginning of La Niña years).
The different composition of NPPpred within the southern Amazon basin is shown in Fig. 8. La Niña years have an overall 250 higher NPPpred, which even lasts throughout the second year of the simulations. However, the drivers causing the difference of increased La Niña predictability are changing over time. At the start of the growing season, which is between December and July, midSOILpred contributes largely to the increased La Niña predictability while deepSOILpred gains in importance around June, when topsoils begin to dry out. An increase in TEMPpred explains a large fraction of increased La Niña predictability throughout the first year.
255 Figure 7. Difference in NPP predictability (Vc) based on the initial soil moisture. The mean NPP predictability of the first year from the 20% driest initializations are subtracted from the 20% wettest initializations for every grid cell. Red color means higher NPP predictability in wet years and blue color a higher predictability in dry years. Soil moisture from 19 -78 cm depth is used to determine initial conditions. A large fraction of years included in the initializations are ENSO years, where the initial anomaly is further extended through persisting oceanic forcing. The black box and yellow triangle stand for regions examined in the main text.
The increase of midSOILpred during the growing season can be explained by the relationship between precipitation and the change in soil moisture in spring ( Fig. 9 (a)). Although the variability of precipitation is comparable between the ENSO states, there is little change in soil moisture in the La Niña years while the relationship between precipitation and soil moisture change is more pronounced in the El Niño years. The difference in this covariance between the ENSO states is linked to the initial water content (Fig. 9 (b)). The El Niño year is initialized at a depleted state, and precipitation is used to recharge midSOIL.

260
This leads to the translation of the variability in precipitation to a variability in midSOIL. Since midSOIL is saturated at the initialization of the La Niña year, it is hardly affected by the variability of precipitation and the excess water leaves the system as runoff or drainage.
The same mechanism is responsible for the difference in deepSOILpred. As midSOIL dries out during the summer months, NPP is increasingly coupled to deepSOIL. Every ensemble member of the La Niña simulation receives enough precipitation to 265 saturate deepSOIL, thereby reducing its variability, while none of the members in the El Niño year can recharge the soil water deficit. Increased NPPpred in wet years due to TEMPpred can have multiple reasons which are difficult to disentangle. As soil moisture and surface temperature are coupled through evapotranspiration, a reduced variability in soil moisture suggests a reduced variability in temperature as well. Contributing to this effect is the nonlinear mechanism controlling evaporation. At 270 the wet end of the spectrum, evaporation is not limited by soil moisture, meaning that a small variability in soil moisture of a wet soil does not affect evaporation. A counteractive process that might increase predictability in dry years is described by Koster et al. (2011). They suggested that in ecosystems which are generally at the wet end of the spectrum (which is the case for the Amazon basin) land-atmosphere coupling is stronger in dry years when evaporation is limited by soil moisture. This increased coupling can extend TEMPpred by linking it to soil moisture. However, their study was conducted on North America,

275
where land-atmosphere coupling is generally stronger than in tropical rainforest (Guo and Dirmeyer, 2013).
To investigate processes behind the difference in temperature variability per ENSO state, we analyzed the key elements of the surface energy balance. Almost all processes have a continuously higher variability in the El Niño years (Fig. 10).
The strongest difference in variability is in net longwave radiation, but this is most likely an effect of increased variability of surface temperature and not the cause. The SD of net shortwave radiation and ground heat flux are evenly increased by around 280 0.4 W m −2 across the first year. Except for some winter and spring months, the latent and sensible heat fluxes have an increased  variability in the El Niño years. At the peak, difference in variability in August is mostly due to an increased variability in the latent heat flux. As mentioned above, there are also certain regions with an inverse relationship between wetness and NPPpred. These are predominantly in arid regions like northwestern Australia, India, northern Caucasus and the western US (Fig 7). The mecha-285 nisms explaining the increased NPPpred in dry years are exemplified on two initialzations from the dry and wet spectrum in northwestern Australia at 23 • South, 122 • East (yellow triangle in Fig. 7).
This higher NPPpred can be contributed to less variability in deepSOIL and PAR (Fig. 11). The predictability providing mechanism of deepSOIL is comparable with the process in the Amazon basin. With soil moisture dynamics frequently operating at extreme ends of the water holding capacity, the variance can be minimized by all ensemble members being pushed 290 against the boundaries of the system. As opposed to the Amazon basin, in northwestern Australia the ensemble members are are clustered at the dry end of the water holding capacity ( Fig. 11 (a) dry years), while any introduction of soil moisture will increase the variability.
Another difference in NPPpred is caused by a differing variability of PAR ( Fig. 11 (b)). Most dry years have little cloud cover and no restriction of incoming radiation. However, in wet years it is difficult to predict the extent of precipitation and 295 cloud cover, which increases the variability of PAR.
The relationship between initial soil moisture and climate predictability is noted by others. Koster et al. (2011) have determined that, depending on the region, the direction of this relationship can go either way. This asymmetry of predictability is present in areas of high land-atmosphere coupling and is caused by the nonlinear relationship of evaporative fraction with soil moisture. Another study has investigated the predictability of European summer heat and find different weather regime 300 frequencies in initially dry and wet conditions (Quesada et al., 2012). This study adds to the view that predictability is not a mere function of location, but depends on the state of the system and predictability therefore has a strong temporal variability.

Conclusions
In this study, we take a closer look at spatiotemporal patterns of terrestrial CFpred, identify the climatic and environmental sources of predictability, and the feedback mechanisms prolonging the memory of the system. We propose a metric of CFpred 305 weighted by the amplitude of carbon flux anomalies. This metric allows to evaluate the role of different regions and processes to the predictability of the global carbon cycle.
We find that the spatiotemporal patterns of NPPpred and Rhpred are determined by (a) the predictability of the carbon flux drivers, (b) the climatic anomalies caused by low-frequency climate modes as ENSO, (c) the seasonal change of limiting factors and (d) threshold processes and nonlinearity of ecosystem responses.

310
On the global average, NPPpred is to 62% explained by SOILpred, and to 30% by TEMPpred. Rhpred is explained by SOCpred and TEMPpred (50 and 27%) predictability. Decomposing the predictability signal shows there is a high spatiotemporal variability in the drivers of predictability. SOILpred and SOCpred are distributed across all areas of high CFpred, while TEMPpred is mostly to be found in the Northern Amazon basin for CFpred and Southern Africa, North America and Southeast Asia for NPPpred. Rhpred can outlast NPP predictability because SOC, its main drivers, has a much higher anomaly persis-315 tence than the drivers of NPP. On the other hand, NPP is more directly affected by climatic drivers and is therefore able to benefit from the predictability of persisting climatic anomalies like the effects of ENSO. Intraannual variability of CFpred is controlled by the seasonally specific limiting factor of NPP and Rh. This leads to NPP gaining predictability in the dry season, when soil moisture replaces PAR as the limiting factor, while Rhpred has its peak in the wet season, when SOC drives the carbon fluxes instead of precipitation in the dry season. This change in limiting factors is due to the nonlinear relationships 320 of transpiration to soil moisture and Rh to precipitation. Both of these relationships describe a saturation point, at which the variability of moisture (precipitation) becomes insignificant to carbon fluxes. Lastly, interannual variability of NPPpred reveals an asymmetry of predictability driven by initial soil moisture and subsequent precipitation. This effect is caused by ecosystems operating at the boundary conditions of the soil moisture regime. The ensemble members of predominantly wet ecosystems are harmonized in wet years when precipitation exceeds the water holding capacity and excess water is removed through runoff 325 and drainage. The reversed effect applies for ecosystems operating at the dry end of the spectrum. These processes reduce the covariance between precipitation and NPP.
Our results highlight the sources of CFpred and can be used for model development to improve the representation of the terrestrial carbon cycle. Further research could be directed towards the simulation of the ENSO imprint in climate models and the relationship between soil moisture and terrestrial carbon fluxes.