Regional asymmetry in the response of global vegetation growth to springtime compound climate events

Compound climate events can strongly impact vegetation productivity, yet the direct and lagged vegetation productivity responses to seasonal compound warm-dry and cold-dry events remain unclear. Here we use observationally-constrained and process-based model data and analyze vegetation productivity responses to compound events of precipitation and temperature in spring and summer across global mid-to-high latitudes. We find regional asymmetries in direct and lagged effects of compound warm-dry events. In high-latitudes (>50°N), compound warm-dry events raise productivity. In contrast, in mid-latitudes (23.5–50°N/S), compound warm-dry events reduce productivity and compound warm-dry springs can cause and amplify summer droughts, thereby reducing summer productivity. Compound cold-dry events impose direct and lagged adverse impacts on productivity in mid-to-high latitudes, exceeding the impacts from individual cold and dry events. Our results highlight the benefits of a multivariate perspective on vegetation vulnerability as precipitation and temperature often covary and jointly drive vegetation impacts. The direct and seasonally-lagged effects of compound weather and climate events in spring on vegetation productivity vary with latitude and can amplify the effects of individual weather events, according to observationally-constrained estimates and process-based terrestrial ecosystem models.

C limate extremes can strongly impact vegetation productivity 1 . Vegetation impacts arising from individual climate extremes such as droughts or heatwaves have been well documented [1][2][3][4][5] , including vegetation resistance and resilience to extreme drought 6,7 and lagged summer plantproductivity responses to warm spring 8 . Given that vegetation is synergistically affected by variability of both temperature and precipitation 9 , compound anomalies in temperature and precipitation, such as simultaneous dry and hot conditions, can strongly impact vegetation productivity 10,11 . In fact, at the regional scale there is evidence for direct impacts on vegetation productivity from co-occurring dry and hot extremes [12][13][14][15][16] . A well-known example is the 2003 European dry and hot summer that decreased vegetation productivity by 30%, thereby offsetting four years of net uptake of atmospheric CO 2 over this region 15 . Similarly, the 2010 European drought and heatwave caused record-breaking reductions in plant-productivity due to heavy precipitation deficit and strong heat stress 12,16 . Notably, while several studies analyzed and projected the impact of simultaneous dry and hot conditions on vegetation productivity 3,17,18 , little attention has been devoted to investigating the lagged response of vegetation productivity to compounding climate conditions such as summertime vegetation productivity responses to springtime compound events, and regional asymmetries of such responses.
Cold spell is another form of extreme event that can cause vegetation damage by leaf frostbite, shortening growing season, and reducing photosynthetic carbon uptake 19,20 . The vegetation productivity responses to compound events that include cold extremes have received little attention so far. Despite ongoing global warming, understanding the impact of cold-related extreme events to vegetation productivity is important, especially because cold extremes still occur frequently 21 and may even increase in frequency in the near future in mid-latitudes as a response to rapid arctic warming [21][22][23] .
Here, we investigate the direct and lagged responses of vegetation productivity to seasonal compound warm-dry (CWD) and cold-dry (CCD) events across global mid-to-high latitudes over the last three decades. We used seasonal compound event indices to model concurrent climates in spring and summer. To characterize vegetation growth and productivity, we employed satellite-derived normalized difference vegetation index (NDVI) and leaf area index (LAI) in combination with observationally-constrained and state-of-the-art model-based gross primary productivity (GPP). Furthermore, highresolution observational and model-based soil moisture and runoff were used to investigate key mechanisms associated with the hydrological responses.

Results
Illustration of the compound event indices. Building on earlier studies 24,25 , we develop two univariate indices to model concurrent climate conditions, i.e., a CWD index that varies from compound cold-wet conditions to CWD conditions, and a CCD index that varies from compound warm-wet conditions to CCD conditions (see "Methods"). The two indices incorporate the dependence between temperature and precipitation and are a measure of how warm/cold and dry a point is relative to the distribution of climate conditions at a given location. We illustrate the two indices on two grid points that have strong but opposite temperature-precipitation correlation. In the case where temperature and precipitation are strongly negatively correlated, the CWD index is well aligned with the primary axis of the bivariate distribution (Fig. 1a). In the case where temperature and precipitation are strongly positively correlated, the same holds for the CCD index (Fig. 1d). As illustrated for several concurrent hot-dry and cold-dry events that occurred around the globe, the two indices well capture these events (Supplementary Figs. 1 and 2).
Notably, in the case where precipitation and temperature are strongly positively correlated, the CWD index indicates the relative anomalies of bivariate joint distribution, and some counterintuitive situations might occur relative to the univariate marginals (Fig. 1b). For instance, points might be labeled as strong CWD events (CWD index > 1.5) even though temperature is anomalously cold (temperature anomalies < 0, red dots in lower left quadrant of Fig. 1b). The CCD index exhibits similar behavior (Fig. 1c). This indicates an interesting property of the compound indices to identify strong compound conditions relative to bivariate distribution that are not necessarily extreme from a univariate perspective 3,24,26,27 .
Widespread direct and lagged impacts of springtime compound climate conditions. To evaluate the lagged summer vegetation responses to spring compound climate conditions, we compute partial correlation between CWD (CCD) spring and subsequent summer vegetation variation by controlling for the influence of summer compound climate conditions on these correlations (see "Methods"). Results show widespread negative associations between CWD spring and subsequent summer vegetation in the mid-latitudes (<50°N; 19% of the total study area, p < 0.05), indicating that CWD spring inhibits summertime vegetation growth. In contrast, positive associations are found in the high latitudes (>50°N; 8%, Fig. 2a). When considering CCD events, CCD spring primarily exhibits negative lagged effects (18%) whereas positive lagged effects are rare (1%, Fig. 2b). Notably, in line with the fact that vegetation is often impacted by both temperature and precipitation 28 , our results indicate that considering springtime temperature and precipitation concurrently can slightly better explain summertime vegetation variability than considering temperature or precipitation in isolation. Specifically, CWD climate conditions in spring are significantly correlated (p < 0.05) with summer vegetation growth in 27% of the total study area compared to 15% for temperature only and 24% for precipitation only ( Fig. 2a and Supplementary  Fig. 3).
To investigate the direct responses, we then calculate correlations between CWD (CCD) indices and vegetation growth in spring and summer separately ( Supplementary Fig. 4). There are widespread direct positive responses to CWD climate conditions in spring (52% of the study area, p < 0.05) and summer (30%) in the Northern hemisphere, while negative responses for spring (10%) and summer (22%) are distributed in mid-latitudes ( Supplementary Fig. 4a, c). CCD climate conditions predominantly show negative correlations in spring (25%) and summer (20%, Supplementary Fig. 4b, d). We further compare the direct effects of compound climate conditions with that of individual temperature or precipitation and find that CWD climate conditions can well incorporate the direct effects of both temperature and precipitation ( Supplementary Fig. 4a, c and Fig. 5).
Process-based dynamic global vegetation models are commonly used tools to simulate the responses of vegetation growth and ecosystem productivity to climate extremes 8,29 . We thus assess the ability of the current set of state-of-the-art vegetation models (TRENDYv6) 30 to capture the lagged responses to spring compound climate conditions ( Supplementary Fig. 6). Overall, we find that the models tend to overestimate the areal extent with positive lagged effects to CWD springs, and the areas with negative lagged effects to CCD springs ( Supplementary Fig. 6a, d). In contrast, the models and observation-based products generally show good agreement on the areal extent that experience negative lagged effects to CWD springs ( Supplementary Fig. 6b, c). When considering LAI, the results are similar ( Supplementary  Fig. 6e-h).
We finally assess the ability of TRENDYv6 models to replicate direct response to compound climate conditions. Overall, the models overestimate the coverage that is significantly negatively correlated to CCD climate conditions (Supplementary Fig. 7a-d).
In contrast, the areal extents with significant correlations to CWD climate conditions are generally well captured by vegetation models (Supplementary Fig. 7e-h). The results are similar when employing LAI (Supplementary Fig. 7i-p).
Regional asymmetries in direct and legacy effects of CWD event. Our results above indicate widespread positive vegetation responses to CWD climate conditions in high latitudes (>50°N), and a negative response in mid-latitudes (23.5-50°N/S). CCD climate conditions mostly exhibit negative effects in mid-to-high latitudes. To understand the effects of CWD and CCD events in more detail and the corresponding physical mechanisms, we perform a regional analysis, focusing on the areas in Fig. 2a where summer vegetation responds positively (r ≥ 0.22) and negatively (r ≤ −0.22) to spring CWD climate conditions, and the areas in Fig. 2b where summer vegetation responds negatively (r ≤ −0.22) to spring CCD climate conditions. Across the focus areas, we investigate the composites of average anomalies in GPP, LAI, soil moisture, runoff, and terrestrial water storage (TWS) (see "Methods"). Notably, over these areas, we find that the direction of direct and lagged productivity responses between the models and observation-based products overall shows good agreement, despite the difference in response magnitude (Figs. [3][4][5].

CWD events increase vegetation productivity in high latitudes.
We first analyze the direct responses of productivity to springtime and summertime CWD events across high latitudes (>50°N, Fig. 3). Productivity increases during CWD spring and summer ( Fig. 3a-c), which is consistent with vegetation responses (Supplementary Fig. 8a-c). Despite elevated spring greenness, spring water overall shows positive anomalies during CWD spring ( Fig. 3d, f, g, i). This result indicates that spring greenness during CWD conditions is not associated with dry spring across high latitudes, which is further confirmed by similar anomalies in springtime TWS ( Supplementary Fig. 8d, f). In contrast, severe water reduction is found in CWD summer (Fig. 3e, f, h, i). This suggests that despite the beneficial effects of CWD events on productivity in summer, they are associated with summer water deficit.
Next, to analyze the lagged effects of springtime CWD events, we investigate the productivity anomalies in summer under three cases, namely CWD spring but non-CWD summer, non-CWD spring but CWD summer, and consecutive CWD spring and summer. Our results indicate that springtime CWD events have positive lagged effects on summer productivity across high latitudes (Fig. 3). Specifically, we find that during non-CWD summer (that is not favorable for summer vegetation growth) preceded by CWD spring, positive anomalies are still found in summer productivity (Fig. 3a). In contrast, during CWD summer (preceded by non-CWD spring), some models and observationbased products exhibit a reduction in summer productivity (Fig. 3b). We further find that summer productivity highly increases during consecutive events (Fig. 3c). Vegetation anomalies show similar behaviors ( Supplementary Fig. 8a-c). Regarding the lagged responses of hydrological variables, CWD springs followed by non-CWD summers do not lead to water dryness, despite increased vegetation greenness (Fig. 3d, g). The magnitude of summer water deficit is similar for both cases that include CWD summer (Fig. 3e, f, h, i) and is consistent with summer TWS anomalies ( Supplementary Fig. 8e, f). These results imply that in high latitudes, summer water reductions characterized by TWS, soil moisture, and runoff are not associated with increased spring greenness but are primarily caused by summer precipitation deficit.
The productivity responses to compound climate conditions may be stronger than that to individual events through the synergistic effects of temperature and precipitation 28 . To investigate this, we compute the average anomalies in GPP and soil moisture associated with univariate events across the focus areas, which are then compared with the effects of CWD and CCD events in high latitudes (see "Methods"). Warm events can not only directly increase productivity but also show positive lagged effects ( Supplementary Fig. 9a, b). In contrast, dry events reduce productivity ( Supplementary Fig. 9e, f). This indicates that the direct and lagged positive effects of CWD events across high latitudes are mainly dominated by the warm component, while dry conditions have negative effects. Therefore, the warminduced increase in productivity slightly exceeds that associated with CWD events (Supplementary Fig. 9b). Soil moisture under warm springs shows positive anomalies ( Supplementary Fig. 9c, d), while they slightly decline during dry spring ( Supplementary  Fig. 9g, h). This suggests that the positive anomalies in soil water during CWD spring are driven by the warm component.

CWD events reduce vegetation productivity in mid-latitudes.
Here, we first investigate the direct effects of springtime and summertime CWD events across mid-latitudes (23.5-50°N/S). Springtime productivity exhibits little changes during CWD spring (Fig. 4a, c), despite dry spring (Fig. 4d, f, g, i). When considering the direct effects of CWD events in summer, the results are similar, whereas the negative magnitude of productivity in summer is larger than that in spring (Fig. 4b, c). This difference suggests CWD conditions in summer show more adverse effects on productivity than that in spring in midlatitudes. The anomalies in vegetation and TWS are consistent ( Supplementary Fig. 10).
Next, the lagged effects of springtime CWD events in midlatitudes are assessed. In cases with CWD spring but non-CWD summer, summer productivity exhibits slight anomalies (Fig. 4a), with slightly decreased summer water (Fig. 4d, g). Summer productivity and water show much higher reductions in case with consecutive events (Fig. 4c, f, i) than for the case with only CWD summer (Fig. 4b, e, h). These results are supported by the responses of vegetation indices and TWS ( Supplementary Fig. 10), revealing that springtime CWD events in mid-latitudes have negative lagged effects on summer productivity and water availability.
The direct and lagged effects of individual events are finally compared with that of CWD events in mid-latitudes. Dry conditions in spring and summer directly decrease productivity and cause soil water dryness ( Supplementary Fig. 11a-d). Moreover, dry spring depletes soil moisture earlier, which, in turn, causes dry summer and reduction in productivity during non-dry summer ( Supplementary Fig. 11a, c). This indicates that dry springs have negative lagged effects on summer productivity. In contrast, productivity and soil water show positive anomalies during warm springs, while they show negative anomalies in summer ( Supplementary Fig. 11e-h). These results suggest that the direct and lagged negative effects of CWD springs are dominated by the dry component in mid-latitudes, while the warm component mitigates the negative effects of the dry component in spring. Accordingly, the decline in productivity due to dry conditions thus exceeds that triggered by CWD events (Supplementary Fig. 11b).
Decreased vegetation productivity due to the negative synergistic effects of CCD events. Here, we first investigate the direct effects of CCD events across mid-to-high latitudes. Productivity reductions are found during springtime and summertime CCD events (Fig. 5a-c) concurrent with water reductions (Fig. 5). Vegetation and TWS show similar behaviors during CCD spring and summer ( Supplementary Fig. 12). These results reveal that CCD events in spring and summer can impose direct adverse impacts on productivity and soil water across mid-to-high latitudes. The productivity reductions in spring and summer are similar in magnitude (Fig. 5a, b), indicating that CCD events between spring and summer can cause similar damage to productivity.
We then analyze the lagged effects of springtime CCD events. Our results indicate that springtime CCD events show negative lagged effects on summer productivity and cause summer water reductions in mid-to-high latitudes (Fig. 5). Specifically, we find that in cases with CCD spring but non-CCD summer, summer productivity and water exhibit strongly negative anomalies (Fig. 5a, d, g). Moreover, summer anomalies are higher during consecutive events (Fig. 5c, f, i) than the cases including only CCD summer (Fig. 5b, e, h). Vegetation indices and TWS show similar responses ( Supplementary Fig. 12). Our results further indicate that CCD spring has more severe negative lagged effects on productivity than CWD spring. That is, we find that in comparison to cases with preceding CWD spring and consecutive CWD events, summer productivity shows higher reduction in cases with preceding CCD spring and consecutive CCD events (Fig. 4a, c versus Fig. 5a, c). Moreover, in cases with CCD spring but non-CCD summer (Fig. 5a, d, g), summer anomalies are close to those in scenarios with non-CCD spring but CCD summer (Fig. 5b, e, h). The vegetation and TWS anomalies further confirm this situation ( Supplementary Fig. 12a, b, d, e). These results suggest that the lagged effects of CCD spring can be of similar magnitude as their direct adverse effects.
We finally compare the direct and lagged effects of individual events with that of CCD events in mid-to-high latitudes. Cold conditions in spring and summer directly reduce productivity but show weak effects on soil moisture ( Supplementary Fig. 13a-d), and cold spring shows negative lagged effects on summer productivity (Supplementary Fig. 13a). Dry events show direct and lagged negative effects on productivity and soil moisture Fig. 4 The responses of vegetation productivity and hydrological variables to CWD events in mid-latitudes (23.5-50°N/S). a-c The average standardized anomalies (z-score) of GPP during CWD spring but subsequent non-CWD summer (a), non-CWD spring but subsequent CWD summer (b), and consecutive CWD spring and summer (c) for areas in Fig. 2a where summer vegetation responds negatively (r ≤ −0.22) to spring CWD climate conditions. d-f The same as a-c, but for soil moisture. g-i The same as a-c, but for runoff. The bar plots with dash lines (without dash line) indicate the average anomalies of multiple observation-based (model-based) products, and the circles indicate the average anomalies of each product. For details on data see Fig. 3. ( Supplementary Fig. 13e-h). These results imply that the negative lagged effects of CCD springs are dominated by both cold and dry components. The effects of CCD events on productivity mostly exceeds the individual dry or cold impacts ( Supplementary  Fig. 13a, b, e, f).

Discussion and conclusions
Our combined use of observation-and model-based datasets suggests that in high latitudes (>50°N), CWD events can directly and indirectly enhance vegetation productivity. Vegetation in high latitudes is primarily controlled by temperature, though lack of precipitation can suppress plant growth to some degree 2 . Hence, CWD events can directly increase productivity, despite the inhibitory effects of dry conditions. CWD events in spring can promote summer productivity by enhancing spring productivity via increased leaf area which can help to sustain higher photosynthesis levels over summer and thus increase summer productivity 29 . Notably, we find that spring greenness under CWD conditions is not associated with water dryness. Despite decreased precipitation during CWD spring, increased temperature across high latitudes may thaw frozen water in soils, offsetting the lack of water input 31,32 . In contrast, summer water availability is reduced by CWD summer conditions. This suggests that summer (soil) dryness is primarily related to summer precipitation deficit in high latitudes. We further find that the warming-induced increase in productivity exceeds that caused by CWD events. Accordingly, the direct and lagged positive effects of CWD events are dominated by the warm component whereas the dry component can constrain productivity to some degree, partially offsetting the warm-induced positive effects.
Across mid-latitudes (23.5-50°N/S), we find that springtime CWD events can cause and amplify summer drought, leading to summer productivity reduction. Water variation in mid-latitudes strongly responds to precipitation 33 . The lack of precipitation combined with increased ET during CWD spring results in depleted spring soil moisture earlier 34 . Such preconditions ultimately contribute to summer water limitations and reduce summer productivity 14,35,36 . Notably, although extremely high temperatures usually negatively impact productivity in mid-latitudes, this does not hold during spring, where the negative effects of dry events are greater than that of CWD events. Anomalously warm spring temperatures can increase vegetation growth, partially offsetting the drought-related negative effects. Accordingly, increased vegetation greenness during the warmer spring can largely alleviate the negative impacts of dry events on vegetation photosynthesis during the subsequent summer 34,37 . Generally, elevated but non-extreme temperatures seem to compensate some of the effects of decreased precipitation in temperate and boreal areas 38,39 .
In mid-to-high latitudes, we find directly and lagged adverse synergistic effects of CCD events on productivity. Temperate and boreal vegetation are most vulnerable to cold extremes 19 , which often negatively affect productivity and result in extensive reductions in productivity 20 . Moreover, our results suggest that the lagged effects of CCD spring can be of similar magnitude as direct adverse effects. Spring cold events are regarded as climatic extremes with high ecological and evolutionary importance. Destruction of foliage, flowers and unripe fruits caused by spring cold is a sustainable damage for plants because it negatively impacts on nutrient storage, growth, reproduction, leaf development, and even ultimately survival in subsequent years 19 . The coupling of cold with dry condition thus can cause long-term negative effects on vegetation due to vegetation carryover effects.
Terrestrial biosphere models tend to overestimate the area that experiences positive lagged effects to CWD spring, negative lagged effects associated with CCD spring, as well as direct negative effects to CCD conditions. This could affect the model's ability to predict future carbon dynamics under climate changes. However, it is still unclear why these models do not well replicate the observed response to environmental conditions. The bias between observational-based and model products could be due to imperfect model parametrizations and lack of important physical mechanisms. One major factor could be that the seasonal vegetation growth is not well represented in the models. Our results also show a large discrepancy between the areal proportions of negative responses of LAI to CCD condition in observational and modeled products. The overestimation of LAI in the models in their negative response to CCD conditions could be a sign of this bias. Sensitivity analyses and validation of model parameters by Fig. 5 The effects of CCD events on vegetation productivity and hydrological variables in mid-to-high latitudes. a-c The average standardized anomalies (z-score) of GPP during CCD spring but subsequent non-CCD summer (a), non-CCD spring but subsequent CCD summer (b), and consecutive CCD spring and summer (c) for areas in Fig. 2b where summer vegetation responds negatively (r ≤ −0.22) to spring CCD climate conditions. d-f The same as a-c, but for soil moisture. g-i The same as a-c, but for runoff. The bar plots with dash lines (without dash line) indicate the average anomalies of multiple observation-based (model-based) products, and the circles indicate the average anomalies of each product. For details on data see Fig. 3.
satellite and ground measurements could help to narrow down the discrepancies.
Overall, our analyses provide comprehensive insights into how vegetation productivity directly and indirectly responds to compound climate conditions. The results can serve as a basis for more reliable predictions of the summer productivity dynamics based on information from preceding springtime weather conditions.

Methods
Data sources. Long-term and multiple indicators were used as proxies for vegetation growth and photosynthetic activity, namely normalized difference vegetation index (NDVI), leaf area index (LAI), and gross primary productivity (GPP). NDVI dataset was obtained from the Global Inventory Monitoring and Modeling Studies (GIMMS) third-generation NDVI product (NDVI3g) based on retrievals from sensors on the Advanced Very High Resolution Radiometer (AVHRR) 40 . The NDVI3g dataset from Jul. 1981 to Dec. 2015 is available at an 8-km spatial resolution and a biweekly temporal resolution. The NDVI dataset has been corrected to minimize various deleterious effects, such as sensor degradation, orbital drift, and volcanic eruptions 40 .
Three types of global LAI datasets, including GIMMS, Global LAnd Surface Satellite (GLASS), and Global Mapping (GLOBMAP) LAI, were collected from 1982 to 2016. The biweekly GIMMS LAI with a spatial resolution of 8 × 8 km 2 were derived from the GIMMS-NDVI3g data using a set of neural networks 41,42 . The GIMMS LAI has been proved to be excellently accurate by validating against field measurements and other satellite-derived LAI products 42 . Using general regression neural networks, the GLASS LAI from 1981 to 1999 was generated from AVHRR reflectance data, while during 2000-2016 was derived from Moderate Resolution Imaging Spectroradiometer (MODIS) reflectance data, which are available at a 0.05°spatial resolution and 8-day temporal resolution 43 . The high quality and accurate GLASS LAI product is characterized by its spatial and temporal continuity (no gaps or missing values), temporal stability, and representation of vegetation phenology 44 . GLOBMAP LAI dataset with 0.08°spatial resolution was obtained by quantitative fusion of MODIS and AVHRR data; it has 15-day and 8-day temporal resolutions from 1981-2000 and 2001-2016, respectively. The GLOBMAP LAI dataset has been validated with field measurements and other products, demonstrating low temporal noise and high quality 45 .
Five types of GPP data from 1982-2016 were used and previous studies demonstrate that they reflect the impact of climate condition anomalies on vegetation productivity 8,29 . The first monthly GPP data with 8 km spatial resolution was based on the Monteith light use efficiency (LUE) equation but was improved with optimized spatially and temporally explicit LUE values derived from selected FLUXNET tower site data 46 . The GLASS GPP (0.05°spatial resolution, 8-day temporal resolution) was derived from the latest version of Eddy Covariance-LUE model 47 . Satellite-based near-infrared reflectance (NIRv) GPP was generated from AVHRR data and hundreds of flux stations around the world, which is available at a 0.05°spatial resolution and monthly scale 48 . In addition, the ensemble estimates of monthly FLUXCOM GPP (0.5°spatial resolution) from machine learning methods together with two flux partitioning methods driven by ERA5 or CRU meteorological forcing (here abbreviated as Flux-ERA5 and Flux-CRU GPP) are used 49 .
Soil moisture from the Global Land Evaporation Amsterdam Model (GLEAM) was used to investigate key mechanisms associated with the hydrological responses, which has a 0.25°spatial resolution with monthly time step from 1982-2016 50 . The GLEAM data are strongly constrained by observations through assimilating multisource satellite-observed soil moisture, vegetation optical depth (as a proxy for water stress), and snow water equivalents from different satellite sensors. The GLEAM data have been widely employed to investigate hydrological variation due to climate extremes or vegetation dynamics 35 .
Two types of global reconstructed terrestrial water storage (TWS) anomalies and observation-based global gridded runoff with 0.5°spatial resolution and monthly temporal resolution from 1982-2016 were also used to investigate major mechanisms related to soil water resources responses. The first TWS anomalies dataset was reconstructed by a statistical model trained with GRACE observations and forced with meteorological datasets 51 . The second TWS anomalies dataset that can reproduce the strong El Niño signal were reconstructed by combining machine learning with time series decomposition and statistical decomposition techniques 52 . The global runoff dataset (here abbreviated as GRUN) was developed by machine learning algorithms that were trained with observation streamflow, which was assessed with cross-validation against in-situ streamflow, indicating high accuracy and quality 53 .
Gridded monthly precipitation and temporal data were obtained from the Climatic Research Unit (CRU TS4.03) at 0.5°spatial resolution 54 . In addition, landcover data used in this study were based on the MODIS MCD12Q1 Version 6 data product, and land cover classification relies on the annual international geosphere-biosphere program. We aggregated the original classes into forest, shrubland, grassland, wetland, cropland, urban and others (e.g., barren and water body), with 0.5°spatial resolution. We only selected grid cells for which the dominant vegetation type occupied >60% of the grid area. We masked ecosystems dominated by cropland, urban and others, as vegetation growth in these areas is highly affected by human activity.
The fine-scale and/or coarse-scale spatial dataset used in this study were all resample to 0.5°spatial resolution using bilinear interpolation.
Terrestrial ecosystem models. We also used the monthly GPP, LAI, runoff, and soil moisture simulations for 1982-2016 from ten process-based terrestrial biosphere models participating in the TRENDY (trends in net land-atmosphere carbon exchange) v6 project (TRENDYv6) 30 . The models used in this study are CABLE, CLM4.5, DLEM, ISAM, JSBACH, JULES, LPJGUESS, LPXBern, ORCHIDEEMICT, and VISIT. All models were run by the CRU-NCEP V8 climatic datasets, which are based on a convergent product of CRU observational climate data. Global atmospheric CO 2 concentrations were collected from a combination of ice-core records and the NOAA monitoring observations 30 . In TRENDYv6, the same set of factorial simulations was derived by following a standard experimental protocol. TRENDY simulation S2 that was forced by varying both atmospheric CO 2 and climate was used in this study. A detailed introduction on these models is shown in Supplementary Table 1.
Land data assimilation system. We also used the monthly soil moisture and runoff for 1982-2016 from the Global Land Data Assimilation System (GLDAS) 55 and Famine Early Warning Systems Network (FEWS NET) Land Data Assimilation System (FLDAS) 56 . GLDAS is a global offline (uncoupled to the atmosphere) terrestrial modeling system, developed by the Hydrological Sciences Laboratory, NASA. GLDAS uses data assimilation to incorporate satellite and observation in advanced land surface models, including Noah, Catchment, the Variable Infiltration Capacity model (VIC) land surface model, to generate optimum land surface states and fluxes. The GLDAS-2.0 used the Princeton meteorological dataset for the forcing data, and the forcing of the GLDAS-2.1 was NOAA-Global Data Assimilation System atmospheric fields 55 . GLDAS provides a monthly global hydrological components dataset at 1°and 0.25°resolution from 1948 to the near present.
The FLDAS is a variant of the NASA Land Information for optimizing agricultural drought assimilation. The meteorological forcing is derived from Modern-Era Retrospective analysis for Research and Applications version 2. FLDAS produces monthly multi-model and multi-forcing estimates of hydroclimate conditions from 1982 to the near present in a 0.1°resolution 56 . In this study, the root zone soil moisture from GLDAS and top 1 m soil moisture from FLDAS were used.
Compound event indices. Though directly counting individual climate extremes is a straightforward approach to investigate simultaneous extreme events 57 , it results in very few samples, requires very long time series, and can not measure the degree to which compound conditions deviate from the average conditions 58 . The function of probability distribution is a common and effective tool to model the hydrometeorological climate anomalies 59 . For example, the standardized precipitation index, one of the most commonly used indices worldwide for detecting and characterizing meteorological droughts, is developed using univariate distribution, such as gamma probability distribution 60 . Recently, bivariate and multivariate probability distributions, such as copula theory 61 , have been used to construct hydrometeorological indices involving multiple factors to better model more complex individual hazards and compound events 10,24,25,58 . Using copula theory 61 to model the dependence structure of temperature and precipitation is an effective way to capture compound events 10,24,25,58 . Here, we compute two univariate indices to identify compound warm-dry (CWD) and compound cold-dry (CCD) events, i.e. a CWD index that tailored to the variation along compound cold-wet conditions to CWD conditions, and a CCD index that tailored to the variation along compound warm-wet conditions to CCD conditions. The two indices do not only identify concurrent extremes, but characterize the degree to which compound conditions deviate from average conditions of bivariate distribution 24,25,62 . Both indices are thus univariate quantities following standard normal distribution 24,25 . When talking about CWD and CDD events we use the following definition: CWD index ≥ 0.8 indicates CWD climate conditions, CWD index ≤ −0.8 indicates compound cold-wet climate conditions, CCD index ≥ 0.8 indicates CCD climate conditions, and CCD index ≤ −0.8 indicates compound warm-wet climate conditions. The introduction on the two indices in detail is described in the Supplementary Material. Statistical analysis. Our study areas only cover the global mid-high latitudes (i.e. 23.5°to 90°N, and 23.5°to 90°S) because of their distinctive vegetation seasonality 35 . Vegetation productivity usually shows strong seasonality, and high productivity mainly occurs in summer season. Moreover, summer productivity is associated with spring climate conditions to some degree 8 . Therefore, we focus on the impacts of spring and summer compound climate conditions on summer productivity. To evaluate direct and lagged responses of vegetation productivity, we firstly defined the spring in the Northern (Southern) Hemisphere as the period from March to May (September to November), and the summer in the Northern (Southern) Hemisphere as the period from June to August (December to February).
To investigate the lagged responses of vegetation productivity to compound climate conditions, we compute the partial correlation coefficient (r) between the summer response variables of interest (NDVI, LAI or GPP) and spring CWD index and CCD index, separately. Partial correlation analysis can robustly measure the degree of association between two random variables, with the effect of controlling random variables removed 8,35 . All variables were detrended before performing partial correlation analysis, and in these partial correlations, the covarying influences of summer CWD index (or CCD index) on the correlations between spring CWD index (or CCD index) and subsequent summer response variables have been removed. Similar partial correlation analyses are performed for spring temperature (or precipitation) to compare the lagged effects of individual climate conditions with that of spring compound climate conditions, in which the covarying influences of summer temperature and precipitation on the correlations between spring temperature (or precipitation) and subsequent summer NDVI have been removed. In addition, the areal coverage corresponding to different correlations values is calculated as a fraction of total study area that was calculated by a function of latitude.
To understand direct and lagged responses of vegetation productivity to CWD and CCD events in detail and investigate key mechanisms associated with hydrological responses, we performed a regional analysis, focusing on the composites of anomalies in GPP, LAI, soil moisture, runoff, and TWS during specific spring and subsequent summer across the areas with positive (negative) partial correlations under three scenarios from 1982 to 2016. All variables were detrended and then were transformed into standardized anomalies. The standardized anomaly (z-score) uses the formula of Z i ¼ ð X i ÀX σ Þ, where Z i denotes spring (or summer) anomaly of response variable in year i, X denotes the spring (or summer) average value of the response variable for the entire study period; σ stands for the spring (or summer) standard deviation of the response variable for the entire study period. The three scenarios are defined as: Scenario I) springtime CWD index ≥ 0.8 and summertime CWD index < 0.8, Scenario II) springtime CWD index < 0.8 and summertime CWD index ≥ 0.8, and Scenario III) both spring and summer CWD index ≥ 0.8. The three scenarios of specific springs and subsequent summers are also employed to the CCD index. The three scenarios are summarized in Supplementary Table 2. The constraints of compound climate conditions can affect vegetation productivity despite not very strong association (i.e. p < 0.2), the regional analysis is thus performed across the areas in Fig. 2a with positive (r ≥ 0.22) and negative (r ≤ −0.22) lagged responses to spring CWD climate conditions, and the areas in Fig. 2b showing negative (r ≤ −0.22) lagged responses to spring CCD climate conditions. Because few areas show positive correlations between spring CCD climate conditions and summer NDVI (Fig. 2b), these areas are not included in the regional analysis. The mean anomalies of response variables were finally computed across positive (negative) correlated areas by averaging all cases meeting the requirement for each scenario from 1982 to 2016.
In comparison to the direct and lagged effects of CWD and CCD events, the mean anomalies of response variables (GPP and soil moisture) associated with individual events (i.e., warm, cold, and dry events) were calculated across the focus areas under scenarios I and III. The springtime and summertime dry events were identified by the standardized precipitation index ≤ −0.8; the springtime and summertime warm events were characterized by the standardized temperature index (STI) ≥ 0.8; the springtime and summertime cold events were identified by STI ≤ −0.8.