# Rising water-use efficiency in European grasslands is driven by increased primary production

Christian Poppe Terán, Bibi S. Naz, Alexander Graf, Yuquan Qu, Harrie Jan Hendricks Franssen, Roland Baatz, Philippe Ciais, Harry Vereecken

<figure><img src="image_1.png" /></figure>

## To cite this version:

Christian Poppe Terán, Bibi S. Naz, Alexander Graf, Yuquan Qu, Harrie Jan Hendricks Franssen, et al.. Rising water-use efficiency in European grasslands is driven by increased primary production. Communications Earth & Environment, 2023, 4 (1), pp.95. 10.1038/s43247-023-00757-x. hal-04069263

HAL Id: hal-04069263

https://hal.science/hal-04069263

Submitted on 6 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Rising water-use efficiency in European grasslands is driven by increased primary production

Christian Poppe Terán<fnref n="1" />✉, Bibi S. Naz<fnref n="1" />, Alexander Graf<fnref n="1" />, Yuquan Qu<fnref n="1" />, Harrie-Jan Hendricks Franssen<fnref n="1" />, Roland Baatz², Phillipe Ciais³ & Harry Vereecken<fnref n="1" />

Water-use efficiency is the amount of carbon assimilated per water used by an ecosystem and a key indicator of ecosystem functioning, but its variability in response to climate change and droughts is not thoroughly understood. Here, we investigated trends, drought response and drivers of three water-use efficiency indices from 1995–2018 in Europe with remote sensing data that considered long-term environmental effects. We show that inherent water-use efficiency decreased by –4.2% in Central Europe, exhibiting threatened ecosystem functioning. In European grasslands it increased by +24.2%, by regulated transpiration and increased carbon assimilation. Further, we highlight modulation of water-use efficiency drought response by hydro-climate and the importance of adaptive canopy conductance on ecosystem function. Our results imply that decoupling carbon assimilation from canopy conductance and efficient water management strategies could make the difference between threatened and well-coping ecosystems with ongoing climate change, and provide important insights for land surface model development.

<footnote n="1">¹Institute of Bio- and Geosciences: Agrosphere (IBG-3), Research Centre Jülich, Jülich, Germany. ²Research Platform Data Analysis & Simulation, Leibniz Centre for Agricultural Landscape Research (ZALF), Müncheberg, Germany. ³Laboratoire des Sciences du Climat et de l'Environnement LSCE, CEA CNRS UVSQ, Centre d'Etudes Orme de Merisiers, Gif-sur-Yvette, France. ✉email: c.poppe@fz-juelich.de</footnote>
The exchange rate of carbon between the land surface and the atmosphere determines the productivity of natural and cultivated ecosystems and influences the concentration of CO₂ in the atmosphere¹⁻⁵. Through photosynthesis, vegetation drives the largest CO₂ flux from the atmosphere to the earth surface²,⁴, which is coupled to a simultaneous vapour outflow though open stomata⁶,⁷. Combining those fluxes by relating carbon assimilation flux to water loss, the water-use efficiency⁶,⁸,⁹ (WUE) is a commonly used indicator to examine changes of ecosystem function³,¹⁰⁻¹⁸ and ecosystem performance and fitness¹⁹. A recent study suggested that WUE constitutes one of three major axes of ecosystem function variability²⁰.

WUE increased due to stomatal closure¹¹,¹²,¹⁴,²¹,²² in response to increasing atmospheric CO₂ since the beginning of the twentieth century. This effect of CO₂ on WUE is still not fully quantified and there is no consensus on its magnitude²¹,²³. Moreover, changes of water supply and demand such as precipitation, soil moisture and vapour pressure deficit (VPD) have been identified as key factors influencing ecosystem functioning and WUE variability³,¹⁸,²³⁻²⁷. However, past large-scale spatiotemporal WUE analyses from remote sensing (RS) derived data rarely included long-term environmental effects²³, such as rising atmospheric CO₂ and atmospheric dryness, and differ even in the sign of detected WUE trends¹⁰,²⁴,²⁸. In addition, RS WUE studies investigating effects of a deficit of water resources, i.e. droughts, on WUE did not yet account for the plant physiological response to increased VPD during droughts by use of an adequate WUE index²⁹⁻³¹. Hence, the spatiotemporal heterogeneity of WUE trends and effects of droughts remains poorly understood²⁵,³⁰⁻³². Their quantification from most recent RS data is important to enhance the understanding of impacts and feedbacks of climate change on ecosystem functioning across data sources and scales²³,³³, as well as to provide insights to adequate functional discretizations and related parameter in land surface models³⁴,³⁵.

Especially in Europe, where droughts were projected to increase in frequency and severity³⁶, a detailed assessment of the WUE ecosystem function variability is required. Ranging from semi-arid to very humid hydro-climates and from croplands and grasslands to evergreen forests, the pan-European domain also enabled analyses regarding the modulations of hydro-climate and land cover to the trends and drought response of WUE that have been found in past studies¹⁸,²³,²⁷,³¹,³⁷,³⁸. Nevertheless, a clear picture of the direction, magnitude and drivers of WUE trends and drought response in European hydro-climates and land cover has not yet been drawn. On top of that, the increasing availability of in-situ data on the continent provides a unique opportunity to assess the RS data quality and on the other hand the geographical, land-use and hydro-climatic coverage of those in-situ station networks.

In this study, we defined three WUE indices. Two by dividing carbon assimilation by ecosystem water-use, i.e. evapotranspiration (ET) and transpiration (T~r~), resulting in EWUE and TWUE, respectively. Essentially, VPD has distinct effects on the carbon assimilation and the water-use of an ecosystem and is expected to increase during droughts. Therefore, the third WUE index relates carbon assimilation to canopy conductance (G~c~) (IWUE) with a linear³,¹⁴,³⁹ implementation of VPD and captures the physiological response to drought by accounting for the regulation of transpiration to atmospheric water demand through adaptive stomatal opening. A more detailed description and their calculation is given in the Method section. Here we compared the trends and drought response of these WUE indices across hydroclimates and land cover to identify the adequate WUE index for estimations of ecosystem health and performance, as well as land surface model development.

The Global Land Surface Satellite (GLASS) photosynthesis (here rather: gross primary production, GPP) RS data⁴⁰,⁴¹,

includes the long-term effect of CO₂, VPD and temperature by use of the revised eddy covariance light-use efficiency model³⁹. So, we used the GLASS GPP and ET data (0.05° resolution, see Methods for more information on the data) and reanalysis data from the Global Land Evaporation Amsterdam Model⁴² (GLEAM, T~r~, 25 km resolution), ERA5-Land⁴³ (T~r~, soil moisture, leaf area index, 9 km resolution) and COSMO-REA6⁴⁴ (meteorology, 6 km resolution) from 1995 until 2018 over Europe. We aggregated the original data to 3 km spatial resolution (to match future land surface model outputs) and 8-day temporal resolution (the coarsest resolution from the original data from GLASS RS). Relationships between ecosystem processes might vary with the seasons, e.g. different effect of droughts on WUE between during vegetation dormancy in winter and in summer, so we masked the data to the meteorological summer (June, July, August) to capture the most active period of the vegetation and to have a consistent signal of WUE drought response.

Finally, we (1) identified trends of summer WUE and differences between WUE indices and land cover types, (2) detected the response of WUE to precipitation and soil moisture in summer drought conditions for hydro-climatic zones and (3) analyzed drivers of summer WUE variability using a discovery method for causal networks of ecosystem processes (Peter & Clark Momenitory Conditional Independence⁴⁵, PCMCI + ⁴⁶). Consequently, this work quantified regional variability of the summer WUE ecosystem function along recent decades and during droughts and identified distinct impacts of environment and ecosystem processes for different hydro-climate zones and land cover types. This work is important for the ongoing development of land surface models that incorporate functional variability to environmental factors and the assessment of ecosystem health and performance in the light of threatening changes in climate and environment.

## Results

**All WUE indices decrease in Central Europe.** Several aspects of ongoing global change (e.g. anthropogenically driven increased CO₂ and resulting changes and temperature and VPD) were found to affect the WUE³,¹⁰,¹⁴,²⁶,⁴⁷,⁴⁸. To identify resulting long-term changes we conducted pixel-wise seasonal Mann-Kendall trend analyses that show spatially extensive, significant (*p* < 0.05) negative summer EWUE trends in Eastern and Central Europe (Fig. 1a).

There was a radial gradient from these continental regions outward towards significant increases in North and South Europe. The median relative change of the whole domain was −5.2%. The total relative changes here and later refer to the ratio of median of slopes of corresponding pixels to the median of intercepts of the same over the whole 24-year period 1995 – 2018. Therefore, the continental WUE change was dominated by the extensive negative trends in Eastern and Central Europe (Supplementary Fig. S1 for the extents of the regions) which changed by −17.2% and −8.8%. The negative median decreasing trend agrees with MODIS products analyses by Tang et al. 2015²⁸ who found decreasing global EWUE, although no resolved signal for Europe was shown. Other studies using MODIS data²⁹,⁴⁹ outlined decreases of a much smaller extent in Eastern Europe until 2014 than we found here, but showed positive trends over most remaining parts of Europe. Importantly, the inclusion here of (1) the long-term influences of CO₂ and VPD to GPP and (2) recent years, 2014 until 2018, might be reasons for differences to mentioned studies. The decreasing trends in EWUE are also evident in station scale and spatially upscaled EWUE estimates from eddy-covariance (EC) measurements which show mostly negative trends in Europe, and a small increase in the southern,
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Albeit different directionality of trends elsewhere, all indices agree on decreasing summer WUE in Central Europe. In the summer seasons (June, July and August) over the years 1995–2018, EWUE (a), TWUE (b) and IWUE (c) trends were calculated over monthly time series for each grid cell with the seasonal Mann-Kendall analysis. The cross hatches mark areas where the trends are significant (*p* < 0.05).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2** Contrasting drought response between WUE indices (rows) are larger to soil moisture (SSI) droughts (right column) as compared to precipitation (SPI) droughts (left column). We show the median WUE drought response of the European domain as median EWUE (a, b), TWUE (c, d) and IWUE (e, f) anomaly during precipitation (a, c, e) and soil moisture (b, d, f) droughts.</figcaption></figure>

more arid areas of the continent¹⁰,⁵⁰ as well. However, site network data was not able to confirm the decreases in Eastern Europe due to the lack of stations in that area (Supplementary Fig. S2 for the distribution of EC stations in used data in Europe).

In Fig. 1b, the summer TWUE trends show the same negative sign and a similar magnitude in Central and Eastern Europe as EWUE trends. The different trend slopes of Central and East as compared to Northern and Southern Europe was evident again. However, the extent of significant trends was not as coherent as with EWUE trends: Patchiness of the sign of TWUE trend implies higher spatiotemporal heterogeneity of T, than of ET. Europe had a median TWUE change of -3.0%. Regionally, Eastern and Central Europe TWUE changed by -12.9% and -8.9% respectively. Scandinavia experienced a TWUE change of +10.7% during the whole study period. These TWUE trends agree better than EWUE trends with previous modeling studies⁵¹.

In contrast, the median summer IWUE change across the region was positive (+4.6%). Further, in contrast to EWUE and TWUE, there were significant positive summer IWUE trends in Eastern Europe (Fig. 1c). There, IWUE changed by a median of +5.2% and on the Iberian Peninsula by +10.0%. However, as for EWUE and TWUE, also IWUE shows a negative change in Central Europe by -4.2% over the 24-year study period. We also aggregated the trend slopes over areas of selected land cover types (Supplementary Fig. S3 for a map of land cover types). We found that IWUE of evergreen needleleaf and deciduous broadleaf forests both changed by +3.0%. Further, in croplands IWUE increased by +5.6% and grassland IWUE changed by a higher margin of +24.2%.

An increase of forest IWUE was recorded in previous studies¹¹,¹²,¹⁴,²¹,²³ and mostly explained by long-term elevated CO₂ concentration in the atmosphere. Nevertheless, the magnitude of the forest IWUE increase from tree-ring studies is still debated⁵². Our results agree with the directionality of trends from eddy covariance towers of Keenan et al.¹⁴ and Mastrotheodoros et al.⁵³, but here we had a more extensive time period including more recent years and state a substantially lower magnitude of IWUE trends (0.12% year⁻¹) in forests. Frank et al.¹² also found IWUE increases in broadleaf and needleleaf forests solely driven by CO₂ based on tree ring isotope measurements. Here we did not find a distinct magnitude of increasing trends between needleleaf and evergreen forests, in contrast to tree-scale studies²²,²³.
**Contrasting drought response between WUE indices.** We further analysed the Europe-wide impacts of precipitation and soil moisture summer droughts on WUE anomalies. We first calculated the standardized precipitation index (SPI) and the standardized soil moisture index (SSI) over the whole study period for each grid cell on a monthly time scale (see Methods for more details). Then, we determined monthly anomalies of each WUE index when the drought index indicates drier than normal conditions, i.e. when it is lower than -1. Finally, the summer WUE drought response (Fig. 2) at each location is the median of its drought anomalies over all summer seasons.

Generally, we found a larger response for SSI (Fig. 2b, d, f) than for SPI droughts (Fig. 2a, c, e) with all WUE indices. EWUE and TWUE both respond predominantly negative to drought, both on or over 80% of Europe for SPI, and on over 78.9% and 66.3% for SSI droughts, respectively. For each grid cell we estimated the relative drought response by relating the median drought response to the overall median WUE over the study period for that same grid cell. The median relative SPI drought response of EWUE and TWUE over the whole domain was -8.6% and -8.5%, respectively, and the response to SSI droughts were slightly larger for EWUE (-9.4%) and lower for TWUE (-7.6%). Notably, largest regional response of both, EWUE and TWUE, happened in Central Europe during SSI droughts, -15.3% and -14.4%, respectively. The Mediterranean, however, was the only region with a median positive TWUE response, which amounted to +2.6% during SSI droughts. In contrast, IWUE responded positively to drought virtually ubiquitous (over 90% of Europe for both drought types), so that the median response to SPI droughts was +17.3% and to SSI droughts was +23.6%. Same as with the other WUE indices, the response of IWUE was larger during SSI droughts. The highest regional response was in the Mediterranean, accounting for +31.7% and +23.1% during SPI and SSI droughts. Note a gradient from Northern humid regions with lower drought response toward higher response in drier areas in the South.

We found agreement on a generally negative EWUE drought response with other past studies based on site scale FLUXNET and spatial MODIS data¹⁷,²⁹. In particular the study of Peters et al.¹⁵ based on atmospheric carbon isotope data also found increased IWUE during severe drought events for the Northern hemisphere. Using remote sensing data of ET and $T_r$ and site-scale GPP, Gu et al.²⁵ discovered discrepant EWUE (positive) and TWUE (negative) drought response. However, our results showing mostly negative EWUE drought response of forests disagree with Huang et al.²⁹, who used MODIS data and found mostly increasing EWUE in Europe. The usage of data of different scales (grid cell remote sensing water-use to site scale GPP) as well as analysis exclusively at EC site locations in their analysis could explain the different EWUE drought response between their analysis and here.

**Gradient of WUE drought response across hydro-climates.** We then analysed the WUE response to different summer drought categories in hydro-climates (Fig. 3) to identify their relationship as well as conditions with an optimal WUE³,¹⁷,⁵⁴. For that, we determined zones of hydro-climates based on annual precipitation⁵⁵ (Supplementary Fig. S1 and Table S1) and calculated median WUE anomalies (defined here as drought response) within the hydro-climate during occurrences of seven categories of summer SPI and SSI indices ranging between $\le -2$ and $> 2$ (Supplementary Table S2) on a monthly period.

In Fig. 3a we uncovered two gradients of EWUE response to SPI categories: (1, along the x-axis) from negative response during dry conditions to positive response during wet conditions with

larger response during more extreme conditions and (2, along the y-axis) increasing negative drought response in dry conditions across all hydro-climates. Here for EWUE, gradient (2) is exhibited by low response in very humid and arid and the largest response in semi-humid hydro-climates. Both these gradients were more pronouncedly evident for SSI categories (Fig. 3b). Zhao et al.¹⁷ also showed increasing negative EWUE drought response with drought severity from FLUXNET data, but in contrast to our study, they showed negative EWUE response during wet conditions.

Similarly, TWUE response also followed both gradients (1) and (2) in very humid to semi-humid regions along SPI (Fig. 3c) and SSI (Fig. 3d) categories. But in semi-arid and arid hydro-climates, we found generally smaller response of EWUE and TWUE along SPI categories and a reversed gradient (1) for TWUE along SSI categories, i.e. positive response during dry and negative during wet conditions.

This reversed gradient (1) was again observed in IWUE response but across all hydro-climates for both drought types (Fig. 3e and f). Gradient (2) in the IWUE response to SPI droughts was directed towards the largest positive response in semi-humid and humid hydro-climates, while during SSI dry and wet conditions it was towards semi-arid regions.

**Drivers of WUE change and WUE drought response.** We conducted a causal network discovery analysis to understand the relationships between meteorological, hydrological and biogeochemical variables and how they interacted to cause WUE variability for all summer months (Fig. 4) and how these interactions differed during droughts (Supplementary Fig. S8). In particular, we applied the PCMCI + method⁴⁵,⁴⁶ to create a causal network (see Methods and Supplementary Fig. S11 and S12) out of the input time series of considered variables for each grid cell. We limited this analysis to the link strengths, i.e. partial correlations, where PCMCI + yielded directed, significant ($p < 0.01$) causal links. It is clear that by definition, both, GPP and ET cause WUE variability. However, their individual contribution to WUE changes might differ depending on geography and environment, so that either GPP or ET, or both similarly, are responsible for those WUE changes. Furthermore, it is important to understand that the resulting network does not directly yield the drivers of long-term changes of a selected variable. Rather, the network depicts the drivers of small-scale changes, identifying the contribution of a certain variable at the link source (arrow beginning) to the overall variability link destination (arrow head, see examples Supplementary Fig. S11 and S12). The input time series to the PCMCI + analysis were detrended and anomalous.

First, we scrutinized the causal networks resulting from all summer months (drought as well as nondrought months included). GPP was significantly linked to EWUE variability on a greater share of Europe than ET (Fig. 4d) but both correlated on average and in absolute values almost equally strong with EWUE ($+0.66 \pm 0.21$ and $-0.6 \pm 0.36$, Fig. 4e, uncertainty range refers to the standard deviation). On the Iberian Peninsula, GPP variability was the predominant cause of EWUE changes on 92.9% of its area (where the correlations of the link GPP → EWUE were above +0.66). Around France GPP and ET did both control EWUE. However, in Central and Eastern Europe 94.7% and 86.4% of the regions had rather ET caused EWUE changes, and a much smaller area with strong links between GPP and EWUE.

Generally, PCMCI + determined less grid cells with directed links from GPP towards TWUE and IWUE than to EWUE (compare totals of Fig. 4d, h, l). Further, the mean link strength between GPP and TWUE as well as IWUE was lower than to EWUE (compare Fig. 4e, i, m). While the zones where GPP and
<figure><img src="image_6.png" /><figcaption>**Fig. 3** The magnitude of WUE response to droughts depends on the drought severity and the hydro-climate. Shown here are EWUE (a, b), TWUE (c, d) and IWUE (e, f) median anomalies (color map) to indicate the drought response during instances of standardized precipitation (SPI) (a, c, e) and soil moisture (SSI) (b, d, f) anomaly categories (x-axes) and hydro-climates (y-axes, Supplementary Fig. S1 and Tables S1 and S2).</figcaption></figure>

water-use influence were both strong around France and north of the Black Sea remained when investigating TWUE (Fig. 4b), the predominance of GPP influence in the Iberian Peninsula and the Mediterranean on TWUE was not present. Rather, TWUE was mostly transpiration controlled (Tr→TWUE < -0.66 on 76.8% of Europe). Only in Scandinavia along the coast and on Iceland GPP prevailed primal driver of TWUE. We found again less occasions where GPP caused IWUE (Fig. 4l), and on at least 90% of the total area Gc dominated over GPP in the control of IWUE variability.

Then, we aggregated mean partial correlations from the directed links of the PCOCI + analysis towards WUE (Fig. 4f, g, j, k, n, o) and towards GPP, Tr and Gc (Supplementary Fig. S6 and S7) over hydro-classes and land cover. In grasslands, the GPP → EWUE link was especially strong for all hydro-climates (Fig. 4f), in exchange for a weak link between ET and EWUE (Fig. 4g). Concurrently, also in grasslands, there was a weak connection between temperature as well as soil moisture to GPP (Supplementary Fig. S6a and b) as well as strong Tr dependent Gc (Supplementary Fig. S7b). We found that besides the predominance of the Gc control on IWUE in other ecosystems, in

grasslands and evergreen needleleaf forests GPP still accounted for a fair share of IWUE variability.

Just across hydro-climates, we discovered stronger links from GPP to EWUE (Fig. 4f) and weaker links from ET to EWUE (Fig. 4g) in semi-arid than in more humid areas. An influence of hydro-climate on the GPP to TWUE link was only clearly visible for deciduous broadleaf forests and grasslands (Fig. 4j, both towards stronger GPP → TWUE links in very humid areas), the strongest Tr → TWUE links, however, were in semi-humid and humid hydro-climates (Fig. 4k). Aside from the low average influence of GPP on IWUE across the continent that we described before, Fig. 4n also shows that this link still played a role for evergreen needleleaf forests and grasslands and it was modulated by hydro-climate, so that the weakest GPP → IWUE links were in semi-arid and the strongest in very humid areas. The Gc influence on IWUE shows lower correlations generally in semi-arid zones throughout all land cover and in grasslands across all hydro-climates.

The input for a second PCOCI + analysis were only those values in the time series of the considered variables that concurred with a soil moisture drought (SSI < -1). We did the same processing of the data to create the European interaction
<figure><img src="image_7.png" /><figcaption>**Fig. 4** EWUE change is caused by gross primary production (GPP) on the Iberian Peninsula and evapotranspiration (ET) variability in Central and Eastern Europe, while we determined less directed causal links to TWUE and IWUE. Here we show partial correlations of GPP versus ET to EWUE (a), GPP versus transpiration (T<sub>r</sub>) to TWUE (b) and GPP versus canopy conductance (G<sub>c</sub>) to IWUE (c) from the PCMCI + analysis. Green colors show a stronger link between GPP and WUE variability, and pink colors a stronger link between water-use and WUE, while gray shows low correlations and dark purple strong correlations of both. On the right side, we show for each map the relative area to the total land surface where PCMCI + yielded directed links from GPP and water-use WUE (d, h, l) and the respective mean partial correlations (circle) across the continent with the standard deviation (error range marker) (e, i, m). Lastly, we aggregate the partial correlations of those links (f, g, j, k, n, o) over land cover (x-axes, ENF Evergreen Needleleaf Forest, DBF Deciduous Broadleaf Forest, MF Mixed Forest, GR Grasslands, C Croplands) and hydro-climates (y-axes).</figcaption></figure>

maps of GPP and ET links in Supplementary Fig. S8a, b and c as we did in Fig. 4a, b, and c. For the bars of relative area covered by significant links (Supplementary Fig. S8d, h and l), the mean link strengths (Supplementary Fig. S8e, i and m) and the aggregations of mean link strengths across hydro-climate and land cover (Supplementary Fig. S8f, g, j, k, n and o), we rather show the

absolute difference of the corresponding values between the drought and the non-drought PCMCI + analyses, to highlight the aggregated response of the links to droughts (see Methods for a more detailed description).

The number of grid cells and therefore the relative area with significant links between GPP, water-use or G<sub>c</sub> have generally
decreased (Supplementary Fig. S8d, h and l), but by a rather small magnitude (always less than 5%). Further, the influence of GPP on EWUE increased and the influence of ET on EWUE decreased during droughts (Supplementary Fig. S8e). Geographically, this meant more extensive sole GPP control over EWUE in Southern Europe, and more co-control of EWUE by GPP in Central and Eastern Europe. This was noticeable through more green colors in the South and shifts from pink and purple to dark purple in Central and Eastern Europe between Fig. 4a and Supplementary Fig. S8a. Similarly, GPP influence on drought TWUE increased and $T_r$ decreased (Supplementary Fig. S8i), mostly noticeable in Central and Eastern Europe (compare Fig. 4b and Supplementary Fig. 8b). However, droughts did not change the interactions of IWUE drivers much. The already small GPP influence on IWUE in general decreased during droughts (Supplementary Fig. S8m), but the ubiquitous $G_c$ influence persisted.

Furthermore, the aggregated drought response of GPP and water-use links to WUE along hydro-climate and land cover were more nuanced. We found a contrasting GPP → EWUE drought response between semi-arid (negative) and more humid hydro-climates (positive), particularly in grasslands and croplands (Supplementary Fig. S8f). The response of the ET → EWUE link to drought was also reversed (positive) in semi-arid grasslands than in other ecosystems (negative). In deciduous broadleaf and mixed forests, as well as in croplands, the drought-related link strength reductions of ET to EWUE and $T_r$ to TWUE, were most dominant, while evergreen needleleaf and mixed forests experienced the strongest increases in GPP influence on both EWUE and TWUE. The semi-arid regions showed a contrasting response of the $T_r$ → TWUE link in deciduous broadleaf forests, where $T_r$ influence increased. Again inverted drought response of GPP → IWUE link in semi-arid hydro-climates response was most dominant in grasslands and croplands, and for $G_c$ → IWUE in deciduous broadleaf and evergreen needleleaf forests. The reduction of GPP influence on IWUE during droughts was most accentuated in very humid evergreen needleleaf forests.

## Discussion

To comprehensively describe the variability of the important WUE ecosystem function along recent climate change, we conducted analyses on the summer-time trends, drought response and the causes of variability across the European continent using novel RS and reanalysis data. Here we did not include investigations on particular effects of the newly included influences of CO₂ and VPD in the RS GPP data on our WUE analysis. But these developments are welcome because the missing influences on GPP and resulting uncertainties on long-term WUE trends were a main point of criticism on RS-based WUE studies²³. Nevertheless, the RS GPP model implementations and underlying data uncertainty must still be scrutinized with respect to WUE. Future works might therefore use control data without these inclusions and test their incorporation to assess the effect on WUE variability.

A Mann-Kendall seasonal trend analysis showed EWUE and TWUE decreased in Central and Eastern Europe and increased in Northern Europe and the Iberian Peninsula. IWUE on the other hand mostly increased but aligns with the other indices on decreases in Central Europe. The dominant increase of IWUE was in accordance with previous studies from different data sources and scales¹⁰⁻¹²,¹⁴,²¹⁻²³,⁴⁵,⁴⁶. But here we also showed highly resolved, spatially extensive, significant decreases of all three WUE indices in Central Europe and of EWUE and TWUE in Eastern Europe which could not yet be confirmed in in-situ data yet, due to the lack of eddy covariance measurements there.

As indicated by Kühn et al.¹⁹, increasing WUE is a trait that is related to enhanced plant performance and fitness irrespective of the WUE index at plant scale. We assumed the relationship of WUE on ecosystem performance and fitness to persist on grid cell scale, but found inconsistency of WUE trend direction between indices (negative for EWUE and TWUE and positive for IWUE). Therefore, we concluded that the IWUE index, which considers the eco-physiological adaptation of $G_c$ and hence compares better with plant scale WUE, is best for inference of ecosystem performance and fitness. Consequently, declining IWUE from our trend analysis indicated ecosystems that potentially did not adapt well to environmental changes and have weakened functionality and performance: the Central European region, in the Central Alps, Balkan and western Turkey as well as Southern Finland (Fig. 1c). Irrespective of the hydro-climate, we observed increased IWUE indicated enhanced ecosystem performance, particularly in grasslands (Supplementary Fig. S4c), and decreased performance in mixed forests. Other land cover have contrasting directionality of IWUE trends and performance across hydro-climates. We note that increased plant IWUE is also a response to stress, and increased performance is not directly associated with increased ecosystem health. Nevertheless, plants and ecosystems that can efficiently adapt IWUE positively to changes in environmental conditions have an advantage with ongoing climate change. Future studies might verify the scaling and from plant to ecosystem and grid cell IWUE trends and the inference of ecosystem performance, but for that more extensive in-situ ecosystem observation data is needed. Because of the contradicting trends, drought-response and drivers between WUE indices, analyses on ecosystem performance from large-scale and spatiotemporal continuous data such as remote sensing, reanalysis and land surface model outputs might focus specifically on the IWUE index. These results create concrete functional references for land surface model developers as well as insights into relationships between ecosystem functioning and performance.

We did not find distinct IWUE trends between needleleaf and broadleaf forest land cover types, exhibiting the known discrepancy between leaf-scale and ecosystem scale WUE²³. We hypothesized that the potential bias from representing spatial homogeneity of ecosystems by aggregation of grid cells by dominant land cover is not significant, because EC measurements (which are based on single land cover ecosystems) also did not show significantly different WUE trends between land cover either²³. The potential spatial aggregation bias also averaged off, by use of an adequately high number of grid cells with corresponding dominant land cover and might yield statistically more sound conceptual conclusions on discrete ecosystems than EC measurements from only a few dozen sites with short time series and underrepresented geographical regions and hydro-climates. Still, and importantly, the data used in our study was subject to several uncertainties. Those include lack of representation of vegetation and atmosphere interactions, e.g. the influence of WUE changes on above-canopy VPD that was not present in the meteorological data from COSMO-REA6 used here. Further, the remote sensing and reanalysis models are forced with different input data, so that drought response of leaf area index (and further GPP, ET and WUE) might occasionally be shifted in time. Though, we assessed this effect over an 8-daily time resolution to be negligible. Furthermore, the vegetation response from long-term environmental effects from rising CO₂ and VPD in the used GLASS GPP data might not be reflected in e.g. ERA5-Land and GLEAM vegetation and soil moisture data, and may have caused bias in the TWUE and IWUE trends. We found that the GPP and ET data used in our study compared well with in-situ measurements by validating time-series from corresponding grid cells with time series from eddy covariance sites (See Methods and
Supplementary Tables S3 and S4). The $T_r$ data was compared to plant sap flow measurements upscaled to forest stand (see Methods and Supplementary Table S5) but did not yield satisfactory results. This could be due to uncertainties in the scaling process and the different ecosystem representation of the grid cell data and upscaled sap flux data. However, large-scale spatial and in-situ $T_r$ data alternatives remain scarce. Further development of in-situ and spatial $T_r$ estimates will be appreciated, and in future those should be evaluated on WUE ecosystem function estimations. WUE calculations, trends and drought response then depend on the combination of uncertainties for each variable, amplifying the uncertainty range for the WUE indices. A future study will compare WUE estimations from eddy covariance sites with corresponding grid cells in our data and land surface model output, to evaluate this uncertainty in more detail.

Further we found negative drought response of EWUE and TWUE to both precipitation and soil moisture droughts, while IWUE responded positively (Fig. 2). All WUE indices showed larger response during more severe drought conditions, and a gradient across hydro-climates (i.e. towards highest response at semi-humid hydro-climates for EWUE, humid for TWUE and semi-arid for IWUE). Despite a clear relationship between precipitation and soil moisture, we discovered more pronounced WUE response gradients during soil moisture droughts. We expected this higher sensitivity of WUE to soil moisture than to precipitation due to its direct connection to the vegetation dynamics (vegetation uptake of water through roots in the soil). The intriguing contrasting TWUE drought response between arid and humid hydro-climates, for example, was much more pronounced during soil moisture drought (Fig. 3d) and less distinctive during precipitation droughts (Fig. 3c). It is noteworthy that the differing root depths of vegetation land cover (especially shallow versus deep rooted) might have influenced the drought response signal. We investigated the difference in WUE response between soil moisture droughts in different soil depth layers (see Methods) and we did not discover apparent differences. So, we focused the further discussion only on WUE drought response to soil moisture droughts in the surface layer (0–7 cm depth).

To conclude on these contrasting TWUE drought response between arid and humid hydro-climates, it is helpful to first discuss the gradient of IWUE drought response across hydro-climates. Those were exclusively positive with the largest response in semi-arid hydro-climates (Fig. 3f). We concluded two points from here: Firstly, the physiological drought response of $G_c$ reduction was not accompanied by a proportional decrease of GPP, resulting in the ubiquitous positive response of IWUE. This was in accordance to the current understanding of plant bio-chemistry drought response⁵⁴. Secondly, humid ecosystems responded much less sensitive to droughts than arid. This was evident in weaker responses in humid hydro-climates than in arid for all drought severity categories, and exhibits the two opposing drought water management strategies of plants: The isohydric strategy common in arid ecosystems, where the stomatal conductance is quickly regulated to limit transpiration during droughts, and the anisohydric strategy common in humid ecosystems, where the plants do not regulate $G_c$ until close to dehydration to maximize carbon assimilation³⁵,⁵⁶. Importantly, our analyses were conducted at monthly time scale, where VPD and soil moisture were found to be coupled⁵⁷, and under such conditions, the isohydric – anisohydric framework is a good predictor for $G_c$ variability⁵⁸. The results show that there was still some adaptation of the $G_c$ in humid hydro-climates even during moderately dry conditions, expressed through slightly increased IWUE. This can be the result of upscaling heterogeneity of different plant dehydration thresholds and occurrence of

anisohydric species in humid climates. Additionally, the drought definitions by standardized indices did not differentiate between specific water needs of ecosystem types. Further, the different net water availability during same drought category between locations of different mean soil water magnitudes might also have confounded these gradients. However, the gradient towards more dynamically adapted stomata and $G_c$ and increased IWUE in dry hydro-climates is clear, and there, it even resulted in increased TWUE (Fig. 3d), which could be due to the fact that decreasing $G_c$ let $T_r$ decrease stronger than GPP, leading to a net water saving per carbon assimilated during droughts. That behaviour remained exclusive for arid and semi-arid ecosystems, suggesting that humid ecosystem functioning with their ineffective and risky water strategy behaviour that evolved with water resource abundance³⁶,⁵⁶, will possibly be inhibited by dehydration (and/or shift to more drought tolerant vegetation) along with the expected increase in drought occurrence and severity across the whole European continent. Upcoming studies should project ecosystem functioning and performance with future climate scenarios and assess the state of humid ecosystems by accounting for ecological instead of soil moisture drought categories.

We concluded with analyses of causal networks from PCMCI + to determine differences between the detected drivers of WUE variability across regions, hydro-climates and land cover types in Europe. A striking result was the GPP caused EWUE variability in Southern Europe which is not evident for TWUE variability accordingly (Fig. 4a, b). Therefore, the weaker link between ET and EWUE in Southern Europe comes from the inclusion of evaporation from soil. Especially in semi-arid, less vegetated regions during dry conditions, soil evaporation will increase as long as there is water in the soil but $T_r$ will decrease due to the closure of stomata (as seen in Fig. 3d). Hence, the resulting ET variability is lower than $T_r$ variability in semi-arid regions and resulted in less significantly explained EWUE variability by ET and the predominance of GPP influence on EWUE (Fig. 4f). However, TWUE variability showed the differences between isohydric, more $T_r$-controlled, arid ecosystems and anisohydric, less $T_r$ controlled, humid ecosystems well (Fig. 4b, j and k). In addition, the influence of GPP that was present on EWUE and TWUE, but not obvious for IWUE suggested that IWUE variability was predominantly controlled by $G_c$ through dependence of GPP on $G_c$ (Fig. 4c, m, n and o). Hence, IWUE depended mostly on the water management strategy and the physiological adaptability to environmental changes. Importantly though, in the most humid regions (near the coast of Norway, in Iceland and the Alps), there is high influence of GPP on IWUE, potentially pointing at a decoupling of GPP and $G_c$.

We investigated the large increasing trend of grassland IWUE (Supplementary Fig. S4 and S5) with respect to its drivers. For that, we depicted the linear regressions of relevant link strengths (as independent variables) as a function of the grassland IWUE trend (dependent variable) slope in Fig. 5. Interestingly, the larger the slope of IWUE trend, the stronger was the GPP → IWUE and VPD → $G_c$ and the weaker was the $T_r$ → $G_c$ link. This pattern matched to the regulation of $G_c$ and $T_r$ to dry conditions that we described above, through down regulated $T_r$ variability and more GPP variability. Comparing the trends of relevant variables confirmed the emergent picture: high positive IWUE trends corresponded with positive GPP trends and negative $G_c$ trends. Coincidentally, soil moisture droughts became less likely because of increasing soil moisture trends, but atmospheric water demand, as VPD, rose, causing the described physiological responses of $G_c$ and $T_r$ in the long term. Combined with the GPP increase, which was presumably caused though the enhanced light-use efficiency with rising atmospheric $CO_2$⁵⁹, among other environmental factors, these resulted in the large IWUE increases
<figure><img src="image_8.png" /></figure>

**Fig. 5** Strongly increasing IWUE in European grasslands are the result of regulated T<sub>r</sub> and increasing GPP. We show in (a) the link strengths (y-axis, color lines) as linear regression function of the IWUE trend (x-axis). The 95% confidence interval was plotted but is minimal for all. The slope and p values of the regression are given in the figure legend. In (b) the linear regression of trend slopes of other relevant variables (y-axis, color lines) are plotted against the respective IWUE trend (x-axis). Again, the slope and p values are given in the figure legend. The stronger the positive IWUE trend, we note stronger GPP influence in (a, blue line) and the stronger the GPP trend slope, too (b, red line). At the same time, although vapour pressure deficit (VPD) is rising (b, light blue line), transpiration (T<sub>r</sub>, b, purple line) was regulated by decreasing canopy conductance (G<sub>c</sub>, b, dark blue line), resulting in the strong IWUE increases.

in European grasslands. Future studies should look in more detail on the factors driving the large positive grassland GPP trends and its decoupling from G<sub>c</sub> variability from in-situ data on smaller scales and whether this signal is represented in current land surface models.

We saw the overarching pattern of semi-arid ecosystems reacting distinctively than humid, croplands and forested ecosystems in the drought response of causal networks (Supplementary Fig. S8). Semi-arid ecosystems effectively control T<sub>r</sub> without inhibiting carbon assimilation, which was evident in unchanged GPP link strengths to EWUE and TWUE as shown in Supplementary Fig. S8f and j, but partly strengthening the T<sub>r</sub> → TWUE link (Supplementary Fig. S8k) and generally strengthening the G<sub>c</sub> → IWUE link (Supplementary Fig. S8o). This underlines the conclusions on the opposing TWUE drought response between isohydric ecosystems in arid hydro-climates and anisohydric ecosystems in humid ecosystems through different water management strategies.

We also emphasize the response of causal networks in grasslands to drought that corresponded with the response to long-term increasing VPD described above. Grasslands showed a rather small change in the T<sub>r</sub> → TWUE link (Supplementary Fig. S8k) and no changes to the G<sub>c</sub> variability as well (Supplementary Fig. S8o). Rather, GPP influence on IWUE increased especially in grasslands, again pointing at a potential decoupling of G<sub>c</sub> and GPP through other environmental factors like increased atmospheric CO<sub>2</sub> that would also explain the large increasing grassland IWUE trend along with continuous drying of the atmosphere (see increasing VPD trend).

Long-term environmental changes and resulting in trends of ecosystem variables over the 24-year study period might have had an impact on influences between those variables, thus changing the causal network over time. Those long term effects on causal networks were not considered here. Prospective work will deal with differences in e.g. decadal causal networks around WUE to determine the changes in the drivers over time, that were not considered here. We included a wide range of meteorological, hydrological and ecosystem processes time series in the PCMCI + analysis to comply with the causal sufficiency assumption and have as many possible confounding factors included as possible. But scarce and fragmentary data kept us from including factors that might have played a role on ecosystem processes change: local atmospheric CO<sub>2</sub> concentrations and human forest and crop management, inter alia.

All in all, the study advanced the understanding of the WUE ecosystem function, using new remote sensing data and relating results to existing studies from different data sources and scales. Concretely, we highlight:

1. Decreasing summer IWUE trend as indicator for potentially inhibited ecosystem function and performance in Central Europe;

2. High increasing European grassland summer IWUE trend despite atmospheric drying was caused by maintained T<sub>r</sub> through adaptive G<sub>c</sub> and increased, decoupled GPP;

3. Contrasting TWUE drought response and IWUE drought sensitivity across hydro-climates that demonstrated the contrast between isohydric and anisohydric plant water management.

In this study, we provided references for the further development of land surface models and its discretization and parameterization and assessments of European ecosystem functioning and performance changes. Future work will scrutinize whether these effects are detected by in-situ observations and the representation of these effects by current implementation and parameterization of plant functional types in state-of-the art land surface models.

## Methods

**GPP data.** We used GLASS12B02 V40 (http://www.glass.umd.edu/introduction.html)⁴⁰ GPP data, which is in 0.05° spatial and 8 daily temporal resolution and was created with the revised EC-LUE model⁴¹ and includes long-term effects of CO₂, temperature (T) and VPD on GPP. While we left the assessment of the impact of the inclusion of those long-term environmental variables in GPP on WUE variability for future studies, we still validated this GPP product with FLUXNET2015 and ICOS Drought-2018 data (see their data descriptions below). To facilitate the data handling in this work and streamline comparisons with upcoming outputs from a land surface model, we aggregated this data from the original 0.05° resolution to the 3 km European CORDEX domain⁶⁰ by bilinear interpolation. Then we compared the site GPP time series with those of corresponding grid cells in Supplementary Tables S3 and S4 and the mean ET value of stations against the mean of corresponding grid cells in Supplementary Fig. S10. We also plotted the geographical distribution of R² values between station and GLASS GPP in Supplementary Fig. S2a. We assessed the validation as good (mean R² value of 0.67 across 60 FLUXNET2015 and 0.7 across 52 ICOS Drought2018 sites), but point out potential bias in regions where there less station density or no station data available at all (e.g. Southern and Eastern Europe).

**ET data.** The ET data used in this work is from GLASS11B02 V41 (http://www.glass.umd.edu/introduction.html), is calculated by a multi-model ensemble approach merging five process-based ET data⁴⁰ and has a 0.05° spatial and 8-daily
temporal resolution. Similar to the GPP data described above, we linearly interpolated the GLASS ET data from the original 0.05° resolution to the 3 km European CORDEX grid. Importantly, GLASS ET is given as latent heat in W m⁻², so we converted the values to mm day⁻¹ by multiplying all values by 0.035, assuming a constant enthalpy of vaporization decoupled of variable temperature, which was applied in numerous studies before because of the small effect of variable enthalpy to the overall outcome of the conversion. Additionally to the validation with FLUXNET2015 and ICOS Drought-2018 ET (calculated from latent heat flux, see below) we also compared GLASS ET with reanalysis ET data from GLEAM3.5a (https://www.gleam.eu/) and ERA5-Land (https://cds.climate.copernicus.eu/cdsapp#/dataset/reanalysis-era5-land). In Supplementary Fig. S9 we plot mean FLUXNET2015 ET from stations against mean ET of the corresponding grid cell in GLASS, GLEAM and ERA5-Land ET. We decided to do our analyses with GLASS ET based on this result: The slope of GLASS ET against station ET shows better representation across sites than ERA5-Land and GLEAM. Again, in Supplementary Tables S3 and S4 we present the R² between GLASS ET and station FLUXNET2015 and ICOS Drought-2018 ET time series. There is good agreement in the ET variability of GLASS and station ET, yielding mean R² across stations of 0.74 and 0.71 for FLUXNET2015 and ICOS Drought-2018 respectively. The geographical distribution of R² between stations and GLASS ET is also shown in Supplementary Fig. S2b.

**T<sub>r</sub> data.** Here we used GLEAM3.5a (https://www.gleam.eu/) T<sub>r</sub> data<sup>42</sup> in all our analysis. The data is created in a land evaporation model and uses data assimilation of soil moisture from satellite observations. It is available in 0.25° spatial and daily resolution. We again bilinearly interpolated in space to the 3 km European CORDEX grid, and aggregated the daily data to 8-daily means to make direct comparisons with the GLASS RS data. We compared GLEAM and ERA5-Land (https://cds.climate.copernicus.eu/cdsapp#/dataset/reanalysis-era5-land) T<sub>r</sub> with sap flow measurements from trees that we upscaled to forest stand T<sub>r</sub> from sapfluxnet<sup>61</sup> (http://sapfluxnet.creaft.cat/). The results are shown in Supplementary Table S5. We opted for GLEAM T<sub>r</sub> because of the higher mean R² across stations (0.42) than ERA5-Land (0.39). We add that the comparisons were biased because the upscaling from tree sap flow to tree stand transpiration did only represent one species (although this species dominated the forest) and did not represent understory transpiration at all while GLEAM transpiration represents whole ecosystem transpiration. Still, we were not aware of currently existing, more adequate homogeneous in-situ T<sub>r</sub> data spanning multiple sites to verify transpiration directly and leave this out for future studies.

**Meteorology data.** Meteorological variables we used in this study were used to determine droughts (precipitation), hydro-climates (precipitation, annual sums<sup>55</sup>, see Supplementary Fig. S1 and Table S1) and/or environmental drivers in the PCMCI + analysis (precipitation, temperature, shortwave incoming radiation, wind speed, relative humidity). All those values originated from the COSMO-REA6 (https://reanalysis.meteo.uni-bonn.de/?COSMO-REA6, a meteorological, high-resolution regional reanalysis over the European CORDEX domain in 6 km spatial and hourly time resolution with boundary conditions from ERA-Interim<sup>44</sup>. Again, we bilinearly interpolated the data to 3 km and aggregated to 8-daily means. This meteorological data also forces a land surface model, which will be used to elaborate on the current implementation of the WUE ecosystem function in Europe in an upcoming study. The VPD was calculated from the temperature and relative humidity after Allen et al.<sup>62</sup>:

$$E_S = 0.6108 \times e^{\frac{(T_{r} - T_{0})}{0.373}} \quad (1)$$

$$E_A = \frac{E_S \times RH}{100} \quad (2)$$

$$VPD = E_S - E_A [\text{kPa}] \quad (3)$$

Where E<sub>S</sub> is the saturated vapour pressure and E<sub>A</sub> the actual vapour pressure, T the temperature and RH the relative humidity.

**Hydrology data and drought indices.** We took soil moisture data from ERA5-Land (https://cds.climate.copernicus.eu/cdsapp#/dataset/reanalysis-era5-land) at two different depth layers (0 – 7 cm, 7 – 28 cm). We bilinearly interpolated from the native 0.09° spatial resolution to the 3 km European CORDEX grid and aggregated the daily values to 8-daily means. After we compared the WUE drought response based on drought indices of both, we opted to just use the surface soil moisture layer (0 – 7 cm) because the WUE response to soil moisture drought in both soil layers did not differ. The drought indices of precipitation and soil moisture of both layers were calculated using the climate indices python package<sup>63</sup> by fitting the soil moisture and precipitation of the whole 1995 – 2018 time series in every grid cell into a gamma distribution. Then we calculate standardized anomalies from the transformed time series:

$$SXI = \frac{X - \bar{X}}{\sigma_X} \quad (4)$$

Where X is the time series of the respective water resource (precipitation or soil moisture from both layers), and X̄ and σ<sub>X</sub> the mean and the standard deviation of

the respective time series at the grid cell level. SXI represents the resulting standardized drought index, specific for the water resources, the standardized precipitation index (SPI) and the standardized soil moisture index (SSI) of the surface soil layer.

**Leaf area index and land cover data.** We obtained the leaf area index for the PCMCI + analysis from the ERA5-Land. We bilinearly interpolated from the 0.09° native spatial resolution to the 3 km European CORDEX grid and aggregated the daily data to 8-daily means.

The land cover data is IGBP-Modified MODIS 20-category data from MCD12Q1 (https://lpdaac.usgs.gov/products/mcd12q1v006/). We upscaled the data from 500 m native spatial resolution to our 3 km European CORDEX working grid by selecting the dominant land cover within the native grid cells whose centers lay inside the coarser 3 km grid cell. Therefore, the land cover indications refer in this paper refer to the grid cell-specific dominant land cover type by area.

**In-situ observation data and validation.** We used in-situ observation data from three different networks, that are FLUXNET2015<sup>64</sup> (60 sites, 569 site-years) and ICOS Drought-2018<sup>65</sup> (42 sites, 577 site-years) for GPP and ET (latent heat) sand SAPFLUXNET (54 sites, 180 site years) for sap flow.

The ICOS Drought-2018 data used the FLUXNET format<sup>64</sup> and methodology that is also used in FLUXNET2015, so here we describe how we used them for both. We used the GPP variable estimated by night-time partitioning from net ecosystem exchange with variable threshold on friction velocity dependence ('GPP_NT_VUT_REF' variable). We masked this variable based on estimations of data quality of net ecosystem exchange ('NEE_VUT_REF_QC' variable), which is the measured variable before GPP gets partitioned from it. Similarly, to use in-situ ET observations, we took the measured latent heat flux from those data ('LE_F_MDS' variable) that uses the marginal distribution sampling technique (MDS) for gap filling. This variable was again masked by data quality measure ('LE_F_MDS_QC' variable). Then, we converted W m⁻² to mm day⁻¹ by multiplying all values by 0.035, assuming a constant enthalpy of vaporization decoupled of variable temperature, which was applied in numerous studies before because of the small effect of variable enthalpy to the overall outcome of the conversion.

As in-situ T<sub>r</sub> data remains scarce, we opted to use a pan-European in-situ data of sap flux measurements and upscale them to the tree stand with the method seen in Nelson et al.<sup>66</sup>; We normalized daily tree sap flow values to the unit basal area of each tree and averaged the values for each species present in the stand. In all stands, species where sap flow was measured covered 90% of the total stand basal area. Then we multiplied the species-specific sap flow per basal area by the basal area of each species in the stand and summed the resulting T<sub>r</sub> of all species to obtain the stand level T<sub>r</sub> in mm day⁻¹.

For the validation, we first found the closest grid cell in the 3 km European CORDEX grid, which the spatial data where regridded to before (see above), to the coordinates of the observation station. Closest grid cell means the grid cell where the grid cell center has the shortest distance to the observation coordinates. We extract the GPP, ET and T<sub>r</sub> variables from the aggregated spatial data and calculated R² between the time series of the station and the grid cell level data of the spatial data.

**Determination of WUE drought response.** We then identified the drought response in each location or grid cell as the difference of the median WUE only in drought years to the overall median WUE in that grid cell. For this, we defined a drought where the monthly SSI or SPI is less than -1. Similarly, drought severity categories were defined over the monthly SPI and SSI bins based on thresholds from Supplementary Table S2. Respective drought response refer to the difference between the median WUE of the variable values masked to when the drought index is within the given bin and the overall median WUE at this location.

**Resulting European data for WUE analyses.** We aggregated the spatial data to one single netCDF file. This data is in the 3 km European CORDEX grid, in 8-daily time-steps from 1995 until 2018 and comprises of variables from GLASS remote sensing (GPP, ET), GLEAM reanalysis (ET, T<sub>r</sub>, surface and root layer soil moisture), ERA5-Land reanalysis (ET, T<sub>r</sub>, surface and first sub-surface layer soil moisture, leaf area index), and COSMO-REA6 meteorology (temperature, precipitation, relative humidity, incoming shortwave radiation, wind speed). The data can be made available upon request.

**Water-use efficiency indices.** The WUE indices are calculated from the selected data described above. We calculated the EWUE index based on the GPP and ET from GLASS remote sensing:

$$EWUE = \frac{GPP}{ET} [\text{gC mm}^{-1}] \quad (5)$$

The TWUE index is based on GLASS GPP and GLEAM T<sub>r</sub>:

$$TWUE = \frac{GPP}{T_r} [\text{gC mm}^{-1}] \quad (6)$$
To calculate the IWUE index, that considers the canopy conductance G_c rather than actual water-use per se, we first define G_c as:

$$G_c = \frac{T_r}{VPD} [\text{mm kPa}^{-1}] \qquad (7)$$

With T_r from GLEAM and VPD calculated from COSMO-REA6 (see above). Therefore, IWUE results as:

$$IWUE = \frac{GPP}{G_c} = \frac{GPP}{T_r} \times VPD = TWUE \times VPD [\text{gC kPa mm}^{-1}] \qquad (8)$$

**Trend analysis.** The WUE slope, intercept and p values were calculated in each grid cell using the non-parametric, seasonal Mann-Kendall analysis from the python package pycman.kendall⁶⁷ for monthly values from 1995–2018 masked to the summer months June, July and August. The slope is determined as the Theil-Sen estimator slope and the p value with a two-tailed test. We estimated median trends and intercepts of geographic regions and land cover by first determining the grid cells corresponding to a certain region or to one particular or a group of land cover and then calculating the median slope and intercept values corresponding to those grid cells, respectively. The relative median trend is then calculated by dividing the resulting median trend by the median intercept. We chose the intercept as the relative reference because it represents the regression value at the start of the time series, in our case i.e. beginning of 1995 and by multiplying with the number of years, i.e. 24, we got the total relative change of the trend regression during this study period.

**Causal network discovery analysis.** PCMCI+⁴⁵,⁴⁶ is a causal network discovery method, which can estimate contemporary and lagged dependencies between time series. Important assumptions to entitle causality in PCMCI + is causal and time series stationarity and causal sufficiency. To fulfill causal sufficiency, we included as many as 11 variables related to ecosystem process variability (i.e. GPP, ET, T_r, temperature, relative humidity, VPD, G_c, shortwave incoming radiation, wind speed, precipitation, soil moisture and the leaf area index) and selected possible links of physically plausible relationships using the selectedLinks parameter and a significance level threshold of 0.01. Although some links of the network were not analysed here, their inclusion secured causal sufficiency (consideration of all possible common influencing factors to a variable) and they will be scrutinized in future investigations. We approximated causal stationarity by masking the time series to the summer season to omit changing relationships between variables, e.g. shortwave incoming radiation→GPP in deciduous forests is given in summer but not in winter⁶⁸. We further acknowledged changes in the causal relationships between the variables during special conditions, such as extreme high temperatures⁶⁹, and in the long-term over the 24-year period but we argue that they are negligible in the 8-day timely resolution in our data. Time series stationarity is satisfied by detrending and using seasonal anomalies of the time series of each variable that have already been masked to the summer months, before inputting into the PCMCI + algorithm of the python package tigrstone (https://github.com/jakobrunge/tigrstone). We conducted three analyses whereby we added one respective WUE index to the input variables and adapted the selectedLinks parameter to match the possible links between the different components and their respective WUE index and no lag of causal relationships. Directed causal links between variables are iteratively found with conditional independence tests. We remark that some relationships among the selected links are of nonlinear nature. Nonlinear conditional independence tests in the PCMCI + algorithm exist, however the detection power of linear links will be inhibited. We chose to assume a linear causal network as we are predominantly interested in how the strength of particular links vary across networks aggregations of regions, climates and land cover. Therefore we opted for the linear partial correlation conditional independence test for the cost of potential false positives and false negatives, but is still able to detect small nonlinearities. We minimized the impact of this uncertainty by limiting the network to physically possible and plausible links with the selectedLinks parameter as described above. For the full methodological description of PCMCI + we refer to Runge et al.⁴⁶. The network including all possible links through our selectedLinks control is shown in Supplementary Fig. S11. Importantly, the total possible network may not make sense as such but enables PCMCI + to unravel the different causal networks of hydro-climates and land cover. An exemplary causal graph output from PCMCI + for a grid cell corresponding to the Wüstebach ICOS station (DE-RuW) for the EWUE index is shown in Supplementary Fig. S12.

To identify shifts in the causal network resulting from a particular environmental condition, two causal networks have to be determined: One with time series including only time periods with this particular effect and one with time periods excluding this effect. The differences in the detected links and link strengths between those two networks are the effect of the environmental condition on the causal network. We investigated the effects of drought on the causal network by splitting the grid cell time series of all variables into two: one only containing values corresponding with SSI values lower than −1, and another non-drought time series, containing all other values. Then we subtracted the absolute link strengths of the non-drought from the drought network to have the difference

network representing all link changes during droughts.

$$Link_{response} = Link_{drought} - Link_{non-drought} \qquad (9)$$

Where Link is a chosen link of the causal network, Linkdrought the strength, or partial correlation, of the chosen link in the causal network during droughts and Linknon-drought the strength of the chosen link in the non-drought causal network. Therefore, we applied PCMCI + to each land surface cell in the study grid for three WUE indices and for non-drought periods, drought periods and overall, resulting in a high dimensional output with six causal networks from mentioned variables for each cell.

## Data availability

We made use of only publicly available data without restricted access. In addition to the more detailed information in the method section, in this declaration we state the sites where the data are accessible publicly based on the used variables. The GLASS remote sensing GPP and ET and leaf area index data are available at http://www.glass.umd.edu/introduction.html. The ERA5-Land reanalysis transpiration, ET, soil moisture and leaf area index data are stored at https://cds.climate.copernicus.eu/cdsapp#/dataset/reanalysis-era5-land. The GLEAM reanalysis ET and transpiration can be found at https://www.gleam.eu/. The meteorological variables, temperature, precipitation, relative humidity, wind speed and shortwave incoming radiation are from the COSMO-REA6 reanalysis at https://reanalysis.meteo.uni-bonn.de/?COSMO-REA6. The station data from FLUXNET2015 are accessible at https://fluxnet.eu/data/fluxnet2015-dataset/, the ICOS Drought-2018 at https://www.icos-cp.eu/data-products/YVR0-4898, and SAPFLUXNET at http://sapfluxnet.creaf.cat/.

## Code availability

The python, bash and NCL code generated to aggregate the data, analyse and plot the data can be obtained by request from the author and are published at https://doi.org/10.5281/zenodo.6522085.

Received: 20 July 2022; Accepted: 13 March 2023; Published online: 27 March 2023

## References

1. Green, C. & Byrne, K. A. Biomass: Impact on Carbon Cycle and Greenhouse Gas Emissions. In Encyclopedia of Energy 223–236 (Elsevier). https://doi.org/10.1016/B0-12-176480-X/00418-6. (2004).

2. Keenan, T. F. & Williams, C. A. The Terrestrial Carbon Sink. Annu. Rev. Environ. Resour. 43, 219–243 (2018).

3. Beer, C. et al. Terrestrial Gross Carbon Dioxide Uptake: Global Distribution and Covariation with Climate. Science 329, 834–838 (2010).

4. Friedlingstein, P. et al. Global Carbon Budget 2021. Earth Syst. Sci. Data. 14, 1917–2005 (2022).

5. Reich, P. B. The Carbon Dioxide Exchange. Science 329, 774–775 (2010).

6. Cowan, I. R. & Farquhar, G. D. Stomatal function in relation to leaf metabolism and environment. Symp. Soc. Exp. Biol. 31, 471–505 (1977).

7. Wright, I. J. et al. The worldwide leaf economics spectrum. Nature 428, 821–827 (2004).

8. Hatfield, J. L. & Dold, C. Water-Use Efficiency: Advances and Challenges in a Changing Climate. Front. Plant Sci. 10, 103 (2019).

9. Introduction: Water Use in Relation to Productivity. in Stable Isotopes and Plant Carbon-water Relations 3–8 (Elsevier). https://doi.org/10.1016/B978-0-08-091801-3.50007-6. (1993).

10. Huang, M. et al. Change in terrestrial ecosystem water-use efficiency over the last three decades. Glob. Change Biol. 21, 2366–2378 (2015).

11. Adams, M. A., Buckley, T. N. & Turnbull, T. L. Diminishing CO2-driven gains in water-use efficiency of global forests. Nat. Clim. Change 10, 466–471 (2020).

12. Frank, D. C. et al. Water-use efficiency and transpiration across European forests during the Anthropocene. Nat. Clim. Change 5, 579–583 (2015).

13. Liu, J. et al. Response of global land evapotranspiration to climate change, elevated CO2, and land use change. Agric. For. Meteorol. 311, 108663 (2021).

14. Keenan, T. F. et al. Increase in forest water-use efficiency as atmospheric carbon dioxide concentrations rise. Nature 499, 324–327 (2013).

15. Peters, W. et al. Increased water-use efficiency and reduced CO2 uptake by plants during droughts at a continental scale. Nat. Geosci. 11, 744–748 (2018).

16. Peñuelas, J., Canadell, J. G. & Ogaya, R. Increased water-use efficiency during the 20th century did not translate into enhanced tree growth: Tree growth in the 20th century. Glob. Ecol. Biogeogr. 20, 597–608 (2011).

17. Zhao, J. et al. Physiological and environmental control on ecosystem water use efficiency in response to drought across the northern hemisphere. Sci. Total Environ. 758, 143599 (2021).
18. Belmecheri, S. et al. Precipitation alters the CO₂ effect on water-use efficiency of temperate forests. *Glob. Change Biol.* **27**, 1560–1571 (2021).

19. Kühn, N. et al. Globally important plant functional traits for coping with climate change. *Front. Biogeogr.* **13**, (2021).

20. Migliavacca, M. et al. The three major axes of terrestrial ecosystem function. *Nature* https://doi.org/10.1038/s41586-021-03939-9. (2021)

21. Marchand, W. et al. Strong overestimation of water-use efficiency responses to rising CO₂ in tree-ring studies. *Glob. Change Biol.* **26**, 4538–4558 (2020).

22. Wang, M., Chen, Y., Wu, X. & Bai, Y. Forest-Type-Dependent Water Use Efficiency Trends Across the Northern Hemisphere. *Geophys. Res. Lett.* **45**, 8283–8293 (2018).

23. Lavergne, A. et al. Observed and modelled historical trends in the water-use efficiency of plants and ecosystems. *Glob. Change Biol.* **25**, 2242–2257 (2019).

24. Liu, X., Feng, X. & Fu, B. Changes in global terrestrial ecosystem water use efficiency are closely related to soil moisture. *Sci. Total Environ.* **698**, 134165 (2020).

25. Gu, C. et al. Discrepant responses between evapotranspiration- and transpiration-based ecosystem water use efficiency to interannual precipitation fluctuations. *Agric. For. Meteorol.* **303**, 108385 (2021).

26. Guerrieri, R. et al. Disentangling the role of photosynthesis and stomatal conductance on rising forest water-use efficiency. *Proc. Natl. Acad. Sci.* **116**, 16909–16914 (2019).

27. Mathias, J. M. & Thomas, R. B. Global tree intrinsic water use efficiency is enhanced by increased atmospheric CO₂ and modulated by climate and plant functional types. *Proc. Natl. Acad. Sci.* **118**, e2014286118 (2021).

28. Tang, X. et al. How is water-use efficiency of terrestrial ecosystems distributed and changing on Earth? *Sci. Rep.* **4**, 7483 (2015).

29. Huang, L. et al. A global examination of the response of ecosystem water-use efficiency to drought based on MODIS data. *Sci. Total Environ.* **601–602**, 1097–1107 (2017).

30. Huang, M., Zhai, P. & Piao, S. Divergent responses of ecosystem water use efficiency to drought timing over Northern Eurasia. *Environ. Res. Lett.* **16**, 045016 (2021).

31. Yang, Y. et al. Contrasting responses of water use efficiency to drought across global terrestrial ecosystems. *Sci. Rep.* **6**, 23284 (2016).

32. Wang, M. et al. Divergent responses of ecosystem water-use efficiency to extreme seasonal droughts in Southwest China. *Sci. Total Environ.* **760**, 143427 (2021).

33. Medlyn, B. E. et al. How do leaf and ecosystem measures of water-use efficiency compare? *New Phytol* **216**, 758–770 (2017).

34. Keith, D. A. et al. A function-based typology for Earth's ecosystems. *Nature* https://doi.org/10.1038/s41586-022-05318-4. (2022)

35. Joshi, J. et al. Towards a unified theory of plant photosynthesis and hydraulics. *Nat. Plants* https://doi.org/10.1038/s41477-022-01244-5. (2022)

36. Lehner, F. et al. Projected drought risk in 1.5 °C and 2 °C warmer climates. *Geophys. Res. Lett.* **44**, 7419–7428 (2017).

37. Bai, Y. et al. Variation in ecosystem water use efficiency along a southwest-to-northeast aridity gradient in China. *Ecol. Indic.* **110**, 105932 (2020).

38. Cooley, S. S., Fisher, J. B. & Goldsmith, G. R. Convergence in water use efficiency within plant functional types across contrasting climates. *Nat. Plants* https://doi.org/10.1038/s41477-022-01131-z. (2022)

39. Umair, M., Kim, D. & Choi, M. Impact of climate, rising atmospheric carbon dioxide, and other environmental factors on water-use efficiency at multiple land cover types. *Sci. Rep.* **10**, 11644 (2020).

40. Liang, S. et al. The Global Land Surface Satellite (GLASS) Product Suite. *Bull. Am. Meteorol. Soc* **102**, E323–E337 (2021).

41. Zheng, Y. et al. Improved estimate of global gross primary production for reproducing its long-term variation, 1982–2017. *Earth Syst. Sci. Data* **12**, 2725–2746 (2020).

42. Martens, B. et al. GLEAM v3: satellite-based land evaporation and root-zone soil moisture. *Geosci. Model Dev.* **10**, 1903–1925 (2017).

43. Copernicus Climate Change Service. ERA5-Land hourly data from 2001 to present. https://doi.org/10.24381/CDS.E2161BAC. (2019)

44. Bollmeyer, C. et al. Towards a high-resolution regional reanalysis for the European CORDEX domain: High-Resolution Regional Reanalysis for the European CORDEX Domain. *Q. J. R. Meteorol. Soc.* **141**, 1–15 (2015).

45. Runge, J. et al. Inferring causation from time series in Earth system sciences. *Nat. Commun.* **10**, 2553 (2019).

46. Runge, J. Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. In *Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence* (AUAI Press, 2020).

47. Dekker, S. C., Groenendijk, M., Booth, B. B., Huntingford, C. & Cox, P. M. Spatial and temporal variations in plant water-use efficiency inferred from tree-ring, eddy covariance and atmospheric observations. *Earth Syst. Dyn.* **7**, 525–533 (2016).

48. Zhang, Q. et al. Response of ecosystem intrinsic water use efficiency and gross primary productivity to rising vapor pressure deficit. *Environ. Res. Lett.* **14**, 074023 (2019).

49. Xia, L. et al. Water use efficiency of net primary production in global terrestrial ecosystems. *J. Earth Syst. Sci.* **124**, 921–931 (2015).

50. Zhang, L., Xiao, J., Zheng, Y., Li, S. & Zhou, Y. Increased carbon uptake and water use efficiency in global semi-arid ecosystems. *Environ. Res. Lett.* **15**, 034022 (2020).

51. Huang, M. et al. Seasonal responses of terrestrial ecosystem water-use efficiency to climate change. *Glob. Change Biol.* **22**, 2165–2177 (2016).

52. Gong, X. Y. et al. Overestimated gains in water-use efficiency by global forests. *Glob. Change Biol.* **28**, 4923–4934 (2022).

53. Mastrotheodoros, T. et al. Linking plant functional trait plasticity and the large increase in forest water use efficiency. *J. Geophys. Res. Biogeosciences* **122**, 2393–2408 (2017).

54. Yang, Y. et al. Evolution of stomatal closure to optimize water-use efficiency in response to dehydration in ferns and seed plants. *New Phytol* **230**, 2001–2010 (2021).

55. Jafari, M., Tavili, A., Panahi, F., Zandi Esfahan, E. & Ghorbani, M. Introduction. in *Reclamation of Arid Lands* 1–19 (Springer International Publishing). https://doi.org/10.1007/978-3-319-54828-9_1. (2018)

56. Sade, N., Gebremedhin, A. & Moshelion, M. Risk-taking plants: Anisohydric behavior as a stress-resistance trait. *Plant Signal. Behav.* **7**, 767–770 (2012).

57. Novick, K. A. et al. The increasing importance of atmospheric demand for ecosystem water and carbon fluxes. *Nat. Clim. Change* **6**, 1023–1027 (2016).

58. Novick, K. A., Konings, A. G. & Gentine, P. Beyond soil water potential: An expanded view on isohydricity including land-atmosphere interactions and phenology. *Plant Cell Environ.* **42**, 1802–1815 (2019).

59. Cai, W. & Prentice, I. C. Recent trends in gross primary production and their drivers: analysis and modelling at flux-site and global scales. *Environ. Res. Lett.* **15**, 124050 (2020).

60. Giorgi, F., Jones, C. & Asrar, G. R. Addressing Climate Information Needs at the Regional Level: the CORDEX Framework. *World Meteorological Organization Bulletin* **58**, 175–183 (2009).

61. Poyatos, R. et al. SAPFLUXNET: A global database of sap flow measurements. https://doi.org/10.5281/ZENODO.3971689. (2020)

62. Allen, R., Pereira, L., Raes, D. & Smith, M. FAO Irrigation and drainage paper No. 56. *Rome Food Agric. Organ. U. N.* **56**, 26–40 (1998).

63. Adams, J. climate\_indices, an open source Python library providing reference implementations of commonly used climate indices. (2022).

64. Pastorello, G. et al. The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data. *Sci. Data* **7**, 225 (2020).

65. Drought 2018 Team & ICOS Ecosystem Thematic Centre. Drought-2018 ecosystem eddy covariance flux product for 52 stations in FLUXNET-Archive format. https://doi.org/10.18160/YVR0-4898. (2020)

66. Nelson, J. A. et al. Ecosystem transpiration and evaporation: Insights from three water flux partitioning methods across FLUXNET sites. *Glob. Change Biol.* **26**, 6916–6930 (2020).

67. Hussain, M. & Mahmud, I. pyMannKendall: a python package for non parametric Mann Kendall family of trend tests. *J. Open Source Softw.* **4**, 1556 (2019).

68. Krich, C. et al. Estimating causal networks in biosphere-atmosphere interaction with the PCMCI approach. *Biogeosciences* **17**, 1033–1061 (2020).

69. Krich, C. et al. Decoupling between ecosystem photosynthesis and transpiration: a last resort against overheating. *Environ. Res. Lett.* **17**, 044013 (2022).

70. Jülich Supercomputing Centre. JURECA: Data Centric and Booster Modules implementing the Modular Supercomputing Architecture at Jülich Supercomputing Centre Journal of large-scale research facilities, 7, A182. https://doi.org/10.17815/jlsrf-7-182

## Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 871128 (eLTER PLUS). The authors gratefully acknowledge the computing time granted through JARA on the supercomputer JURECA,⁷⁰ at Forschungszentrum Jülich.

## Author contributions

C.P.T. and B.S.N. conceived and designed the study. C.P.T. processed the data and performed the analyses. B.S.N., A.G., R.B., H.J.H.F., Y.Q., H.V. and P.C. made suggestions to the analyses and helped interpret the results. C.P.T. wrote the manuscript and edited the suggestions from the co-authors.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-023-00757-x.

**Correspondence** and requests for materials should be addressed to Christian Poppe Terán.

**Peer review information** *Communications Earth & Environment* thanks the anonymous reviewers for their contribution to the peer review of this work. Primary Handling Editor: Aliénor Lavergne. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023