Vegetation modulates the impact of climate extremes on gross primary production

primary production Milan Flach1,2, Alexander Brenning2, Fabian Gans1, Markus Reichstein1,3, Sebastian Sippel4,5, and Miguel D. Mahecha1,3,6 1Max Planck Institute for Biogeochemistry, Department Biogeochemical Integration, P. O. Box 10 01 64, D-07701 Jena, Germany 2Friedrich Schiller University Jena, Department of Geography, Jena, Germany 3German Centre for Integrative Biodiversity Research (iDiv), Leipzig, Germany 4Norwegian Institute of Bioeconomy Research, Ås, Norway 5ETH Zürich, Institute for Atmospheric and Climate Science, Switzerland 6University of Leipzig, Remote Sensing Center for Earth System Research, Germany Correspondence: Milan Flach (milan.flach@bgc-jena.mpg.de)


Data
To detect hydrometeorological extreme events we use 2-m air temperature, incoming shortwave radiation (both from ERA5, original resolution 0.25 • , Copernicus Climate Change Service (C3S) (2017)), and surface moisture (v3.2b, original resolution 55 0.25 • from the GLEAM model-data integration framework, (Miralles et al., 2011;Martens et al., 2017)). We consider surface moisture as a hydrometeorological variable due to its importance for drought detection although it is influenced by vegetation.  2010)) We group the available land cover classes in forest ecosystems (land cover classes containing "forest"), agricultural ecosystems (containing "crop"), and, all remaining other land cover types.

Preprocessing and anomaly detection
We compute deviations from a smoothed median seasonal cycle in the hydrometeorological variables, which we denote as anomalies. For detecting extreme events, we apply a multivariate anomaly detection procedure described in detail in (Flach et al., 2018). It (i) accounts for seasonal changes in the variance of the anomalies using a moving window technique, and (ii) uses climatic similarities to obtain more robust thresholds for extreme event detection via spatial replicates as proposed by Mahecha et al. (2017) (for more details see the B).

70
The extreme event detection algorithm itself is applied to the set of hydrometeorological anomaly time series and returns anomaly scores computed by kernel density estimation. Kernel density estimation showed good performance among other possible methods and accounts for nonlinearities in the data (Flach et al., 2017). The resulting anomaly scores can be interpreted as a univariate index of deviation from the general multivariate pattern. We consider the highest 5% of the anomaly scores to be extreme events (95th percentile), which is within the typical range of percentiles defining extreme events (McPhillips et al., .

Framework for extracting event-based statistics
We use the extracted binary information (extreme / non-extreme) to compute statistics based on the spatio-temporal structure of the extreme events similar to (Lloyd-Hughes, 2011;Zscheischler et al., 2013;Mahecha et al., 2017;Chen et al., 2019).
Extreme voxels are considered to belong to the same extreme event if they are connected within a 3 x 3 x 3 (long x lat x 80 time) cube. Note that this definition includes connections over edges. We compute event-based statistics from the 1000 largest extreme events globally as introduced also for the Russian heatwave (Flach et al., 2018). Specifically, we calculate affected volume, centroids, mean and integral of GPP separately for positive and negative anomalies, as well as the distance between the centroids of the positive and the negative anomalies of GPP during the event. We consider an event to be predominantly a   relative drought (relative heatwave) if more than 50% of the surface moisture (temperature) values during the extreme event are 85 beneath (exceed) the 5th (95th) percentile of the variable. We select drought (n = 98) and heat (n = 44) events and combined drought-heat events (n = 71), which are taking place during the growing season (total n = 213), i.e. the centroid of the event is within the half year encompassing the seasonal GPP maximum. Our statistics account for the spherical geometry of the Earth by weighting with the cosine of latitude.
Furthermore, we evaluate if the positive and negative anomalies in GPP during the event predominantly have a spatial or 90 temporal component. Therefore, we split the event in parts with enhanced and parts with reduced productivity. Between those two parts, we compute the spatio-temporal distance between the centroids of each part. We consider positive and negative GPP anomalies to occur predominantly spatially if the temporal distance of the centroids is almost simultaneous, i.e. less than one time step in the data (eight days). GPP anomalies are considered to be predominantly temporally changing if the spatial distance of the centroids is less than 110 km (approximately one degree at the equator). Both, spatial and temporal components 95 can be found for centroids which are more than 110 km and more than eight days away. . Bar sizes are proportional to the affected volume of the identified events. Forests tend to be associated with enhanced productivity rates, while agricultural ecosystems tend to be associated with reduced productivity.

Statistical model of GPP during extreme events
As we detect heatwaves and droughts relative to the mean seasonal patterns, positive or negative GPP anomalies during the droughts and heatwaves may additionally be influenced by differences in the conditions in the hydrometeorological variables during the extreme event, differences in background climate in which the vegetation is growing, or duration and affected area of 100 the event. We use gradient boosting machines (Friedman, 2001) to predict average GPP anomalies during the event as a function of mean surface moisture, mean temperature, mean radiation during the event, duration, affected area, land cover class, and mean climate during the growing season, i.e. mean temperature and surface moisture during all growing seasons between 2003 and 2018. We tune model parameters (shrinkage parameter, depth of the trees, bag fraction, minimal number of observations per node) by following a workflow described in Elith et al. (2008) using a hyper grid search from 100 different random 105 initialisations of splitting the data into training (75%) and testing (remaining 25%). We compute uncertainty of the variable importance measure described in (Friedman, 2001) from each of the 100 best models of the hyper grid search. Additionally we use an approach based on Local Interpretable Model-agnostic Explanations (LIME), which tries to predict each single observation in a black box model based on locally weighted regression (Ribeiro et al., 2016). Here, this approach helps to understand (1) the effect of specific land cover classes, and (2) the direction of the effect.

Results
Our analysis based on a 5% threshold in the multivariate anomaly scores leads to a detection of 213 events (98 relative droughts, 44 relative heatwaves, 71 compound drought-heatwaves) between 2003 and 2018.
If we only discriminate forest and agricultural ecosystems, we find substantial differences in the direction of the GPP anomalies during extreme droughts and heatwaves in the growing season. In agricultural and other non-forest land-cover types, GPP 115 was reduced during the identified events (agricultural land-cover types: 64% (56-72%) reduction, Figure 1 (a); other ecosystems 60% (53-67%), Appendix Figure A1). In forested areas, instead, a majority of 71% (63-78%, 95% confidence interval) of events shows enhanced productivity (Figure 1  The events analyzed here are based on relative radiation, heat and water availability anomalies (see Methods). To better understand the role of absolute climate conditions we show the reported GPP anomalies in the terms of absolute temperatures 125 and surface moisture levels in Figure 3(a). The figure shows that reduced rates of GPP tend to coincide with very low surface moisture and high temperature (eight-daily averages). Delineating different ecosystems within this space shows that they are arranged along decreasing surface moisture values. Most extreme events in forests tend to occur under slightly higher surface moisture conditions compared to agricultural and other ecosystems (Figure 3(b)). Forests are hit less frequently critical dry conditions for which we predominantly observe reduced productivity. In contrast, we observe reduced productivity during the 130 events for agricultural ecosystems, which experience frequently critical hot and dry conditions (Figure 3(b)).
While Figure 3(a) shows that temperature and soil moisture have some effect on the direction of the impacts, they are insufficient to explain the observed patterns in detail. To unravel the importance of land cover type and other factors we predict average GPP anomalies using gradient boosting machines (R 2 = 0.43, Friedman (2001) Section 2.4) and explore their relative variable importance. Growing season temperature, event duration, and land cover type are, in decreasing order, are the most 135 important variables in the statistical model (Figure 4(a)).
Apart from identifying important variables that explain the GPP anomalies during drought and heat anomalies, we disentangle the direction of each factor's effect in the model, and, in particular for specific land cover classes. Whereas growing season temperature and duration show a negative model coefficient, i.e. a longer duration and a warmer climate are associated with a stronger impact, as expected, productivity in different land cover types is influenced in contrasting ways: Land cover  types including forests and woody savannas show increased average GPP during the extreme events. In contrast, agricultural ecosystems (land cover types including crops), grasslands, savannas, and other land cover types reduce average GPP anomalies (Figure 4(b)). Warmer growing season climates and higher temperatures during the event are associated with negative GPP anomalies. In contrast, greater availability radiation and higher surface moisture have a positive influence on the impact.
We showed that land cover type is one of the major factors influencing the GPP anomaly during the event. A single hy-145 drometeorological extreme event can affect two or more adjacent land cover types simultaneously with potentially contrasting impacts (spatial contrasting anomalies), but enhanced productivity can also be observed earlier than reduced productivity (or vice versa, temporally contrasting anomalies). To explicitly quantify the role of spatial vs. temporal effects on the GPP anomalies during extreme events we split each event in parts with enhanced and reduced GPP anomalies and compute the centroidal distance in space and time. In fact, positive and negative GPP anomalies mostly co-occur simultaneously in adjacent spatial 150 regions (116 events of 213 events in total within ± 8 days, Figure 5). Especially for large scale events (large volume), a considerable distance of the anomalies can be observed in space and time. However, taking only the temporal distance into account, we have more events with enhanced productivity before the reduced productivity (temporal distance < −8 days, n = 44) than after (> 8 days, n = 33).  Contrasting responses of ecosystems to climate extremes, e.g. in the US in 2012 (Wolf et al., 2016) or in Russia in 2010(Flach et al., 2018, are not singular cases but are shown to be frequent phenomena in response to hydrometeorological extreme events at the global scale. Within the same extreme event, reduced and enhanced productivity can be observed simultaneously in adjacent spatial regions. This finding complements previous studies on temporal (Wolf et al., 2016;Sippel et al., 2017a;Buermann et al., 2018) or spatial contrasting responses (Jolly et al., 2005;Zaitchik et al., 2006;Lewińska et al., 2016).

160
This study provides evidence that the impacts of extreme drought or heat anomalies on GPP during growing seasons is, firstly a function of event duration and long-term climate, but secondly, also depends on the affected land cover type. In particular the tendency towards positive vs. negative responses seems to be controlled by tree cover (similar to the results of limits are not violated. Yet, the prevalence of certain land cover types is partly controlled by climatic gradients, and therefore land cover cannot really be considered independently of the mean climatological conditions that likewise play a role ( Figure   3(a)). Climate conditions also lead to adaptation of physiological processes. For instance, forests in dry ecosystems may be characterized by a more conservative water use strategy (Teuling et al., 2010;van Heerwaarden and Teuling, 2014;Ramos 175 et al., 2015) and adapted to drought compared to analogous land cover types whose biogeographic history experienced colder and more moderate conditions (Doughty et al., 2015). Moreover, forests have access to deeper soil water compared to other ecosystems (Yang et al., 2016;Fan et al., 2017). The degree of isohydricity may further differentiate the response of forests, as it differs between tree species (Roman et al., 2015;Ruehr et al., 2015;Yi et al., 2017).
Our study only reports on GPP responses during the climatic anomaly without considering the legacy of the events. Re-180 sponses may emerge with some time lag between weeks to months (Schwalm et al., 2012;Ruehr et al., 2015), or even at longer time scales (years) (Saatchi et al., 2013;Anderegg et al., 2015). Hence, finding enhanced productivity of forests during some heat event does not exclude increased mortality in the long-term. Forest ecosystems are known to potentially respond much delayed to environmental stress, which can trigger strong secondary impacts like insect outbreaks (Hicke et al., 2006;Rouault et al., 2006;Allen et al., 2010), or fires (Brando et al., 2014. In contrast, agricultural systems are known to be very directly 185 vulnerable to droughts (De Keersmaecker et al., 2016;Bachmair et al., 2018). We choose the growing season as time period of interest, which is notably different than summer for some regions, e.g. in the Mediterranean where more positive responses to warm anomalies in the cold season may be expected (Sippel et al., 2017b), and also impacts of droughts may be less than during the dry season (Huang et al., 2018). Note that due to complex interactions between GPP and ecosystem respiration no direct translation of the results into net ecosystem exchange is expected (Richardson et al., 2007). 190 Another aspect to discuss is data quality. Gross primary productivity from FLUXCOM-RS may inherit errors from the underlying remote sensing products; these have, in particular, been discussed for tropical forests (Asner et al., 2004;Asner and Alencar, 2010;Wu et al., 2018). Recently, Stocker et al. (2019) showed at the global scale that remote sensing retrieved GPP underestimates drought impacts due to soil moisture effects on light use efficiency. Comparing our estimates of GPP impacts to published data from eddy covariance stations for two case studies (US 2012, (Wolf et al., 2016), and Europe 2003(Ciais et al., 195 2005Reichstein et al., 2007)) indicates that we do indeed underestimate GPP impact. Thus, we suspect that in addition to the GPP estimates used by Stocker et al. (2019), also FLUXCOM-RS GPP underestimates the impacts of climate extreme events specifically for forest ecosystems. As FLUXCOM-RS exhibits a good agreement for forests globally with GPP estimates based on solar-induced fluorescence (Walther et al., 2019), the lack of sensitivity to drought and heat impacts in forest ecosystems may be a more general issue in remote sensing data.

Conclusions
We conclude that a more differentiated consideration of the role of land cover reveals firstly major differences between forest and agricultural ecosystems. These differences may originate from a different (micro-)climate or different water management strategies including the access to deeper soil water or point to more strongly lagged impacts in forest ecosystems. However, the lack of sensitivity of forest ecosystems to droughts and heatwaves is stronger than we would expect it to be. Thus, we think 205 that our results also point towards deficiencies in FLUXCOM-RS derived GPP which are potentially a more general issue in remote sensing derived indices of vegetation activity. These deficiencies call for the development of new global GPP products with a higher sensitivity to droughts and heatwaves, which can unravel the role of forest ecosystems in a more frequently hot and dry future climate.  . by Flach et al. (2018). In summary, the used approach defines climatically and phenologically similar regions by using the leading principal components (here: three) of the seasonal cycles of the hydrometeorological variables (temperature, surface moisture, radiation) in addition to the vegetation proxy (gross primary productivity). Similar cycles appear in the same region of the obtained principal component space ( Figure B1). Thus, a simple classification can be obtained by dividing the principal 220 component space into equally sized cubes. Here we use 25 breaks for each of the first three principal components, which leads to 814 classes globally of similar climate and phenology. For each pixel, we sample four random spatial replicates from each region to efficiently run the following anomaly detection workflow globally (previously the procedure was used for Europe only). 14 https://doi.org/10.5194/bg-2020-80 Preprint. Discussion started: 27 March 2020 c Author(s) 2020. CC BY 4.0 License.