Terrestrial evaporation response to modes of climate variability

Large-scale modes of climate variability (or teleconnection patterns), such as the El Niño Southern Oscillation and the North Atlantic Oscillation, affect local weather worldwide. However, the response of terrestrial water and energy fluxes to these modes of variability is still poorly understood. Here, we analyse the response of evaporation to 16 teleconnection patterns, using a simple supervised learning framework and global observation-based datasets of evaporation and its key climatic drivers. Our results show that the month-to-month variability in terrestrial evaporation is strongly affected by (coupled) oscillations in sea-surface temperature and air pressure: in specific hotspot regions, up to 40% of the evaporation dynamics can be explained by climate indices describing the fundamental modes of climate variability. While the El Niño Southern Oscillation affects the dynamics in land evaporation worldwide, other phenomena such as the East Pacific–North Pacific teleconnection pattern are more dominant at regional scales. Most modes of climate variability affect terrestrial evaporation by inducing changes in the atmospheric demand for water. However, anomalies in precipitation associated to particular teleconnections are crucial for the evaporation in water-limited regimes, as well as in forested regions where interception loss forms a substantial fraction of total evaporation. Our results highlight the need to consider the concurrent impact of these teleconnections to accurately predict the fate of the terrestrial branch of the hydrological cycle, and provide observational evidence to help improve the representation of surface fluxes in Earth system models. Climate cycles cast a far-reaching shadow across evaporation from global land surfaces, but their influence varies across space and time. A team led by Brecht Martens from Ghent University tested the statistical relationship between sixteen climate cycles and land evaporation. They revealed that the El Niño Southern Oscillation (ENSO) has the strongest linkage, or teleconnection, to evaporation. But depending on season and location, other climate cycles influenced evaporation more strongly, in some cases explaining up to 40% of the month-to-month variation in evaporation. Climate cycles influence evaporation by modifying the atmospheric demand for water, rather than altering water supply itself. Yet teleconnections are not universal, and although ENSO may be the most important teleconnection in terms of percentage of land where evaporation is affected, it is not a dominant control of land evaporation. Ultimately, the research highlights the need for a nuanced, multivariate understanding of interactions between climate cycles and land evaporation.


INTRODUCTION
Intra-annual and decadal variability in Earth's climate is largely driven by large-scale and periodic changes in the (coupled) state of ocean and atmosphere, such as the El Niño Southern Oscillation (ENSO). These modes of climate variability do not only affect meteorological conditions in regions where they occur, but can also influence weather patterns in remote areas, being commonly referred to as teleconnections. 1 As such, these modes of climate variability are expected to affect the general dynamics of ecosystems and biogeochemical cycles inland. 2,3 Terrestrial evaporation (or evapotranspiration) is a critical variable in the water cycle that can serve as a diagnostic of ecosystem activity, and that is expected to be influenced by these teleconnections. 4 Terrestrial evaporation consists of the vaporization of water through vegetation (transpiration), directly from soil (bare soil evaporation), or from wet canopies (interception loss) and snowcovered surfaces (sublimation). Transpiration accounts for the vast majority of the flux at the continental scale, 5 and is closely related to the process of photosynthesis-the assimilation of carbon dioxide by plants-hence, acting as a nexus between the water and carbon cycles. 6,7 Due to the high energy requirements to evaporate water, terrestrial evaporation also controls the partitioning of radiation at the ecosystem and the cycling of energy in the atmosphere, thereby affecting air temperature, humidity and cloud formation. [8][9][10][11] The skill of Earth system models to predict future changes in the global water cycle relies on our ability to provide insights on the response of land ecosystems to internal climate variability. 4,[12][13][14][15] As such, quantifying the sensitivity of terrestrial evaporation to the dominant modes of climate variability-and understanding the mechanisms through which these modes affect the flux-is critical to disentangle the effects of climate change on hydrology. 4,[12][13][14][15] Some studies have investigated the effect of specific teleconnections on hydroclimatological variables, such as precipitation, 16,17 temperature, 18 river discharge 19,20 and soil moisture. 21 Others have also highlighted the link between teleconnection patterns and general ecosystem dynamics, at the global-scale. 14,15 However, the response of terrestrial evaporation at global scales remains unknown to a large degree. Yet, a few studies have linked intra-annual variability in terrestrial evaporation to climate oscillations, but only at the regional-scale and focussing on specific teleconnetions relevant for the area of interest (e.g., refs. 22,23 ).
Only recently, the emergence of global-scale, long-term datasets of terrestrial evaporation-derived from satellite and in situ data-has enabled a global analysis of the intra-annual variability in terrestrial evaporation and its link to large-scale modes of climate variability. 4,16,24 The study by Miralles et al. 4 showed a strong control of ENSO over terrestrial evaporation in Australia, southern Africa and eastern South America. Anomalously low volumes of evaporation in these regions during El Niño were related to rainfall deficits in these areas, resulting in higher evaporative stress and thus lower evaporation rates. In fact, ENSO is often considered to be the primary mode of internal climate variability for many regions on Earth. 25,26 However, the overall effect of other fundamental modes on global terrestrial evaporation has not yet been explored. In addition, the dominant effect of ENSO is presumably concentrated in the tropics and the Southern Hemisphere, 4 while in the Northern Hemisphere, weather patterns and ecosystem functioning are likely to be influenced by other teleconnections. 1,15 Therefore, a direct analysis of the effect of teleconnections on terrestrial evaporation is still missing, which hampers the ability to analyze the effects of internal climate variability on the global hydrological cycle. 27 Here, we investigate the impact of leading modes of internal climate variability on global terrestrial evaporation. We use an error-based merger of observational datasets based on the Global Land Evaporation Amsterdam Model (GLEAM), 28,29 the Penman-Monteith-Leuning (PML) model, 24 and the Model Tree Ensemble (MTE). 30,31 A simple supervised learning algorithm, called the Least Absolute Shrinkage and Selection Operator (LASSO), 32 is used to relate the temporal variability in the evaporative flux, and that of its local climatic drivers, to climate indices (CIs) specifically designed to diagnose the modes of climate variability. This approach allows us to disentangle the influence of inter-dependent and cross-correlated processes. 13,33 Ultimately, our objectives are to identify the dominant modes of climate variability controlling land evaporation dynamics across the globe, and to find the local climatic drivers through which these modes act upon evaporation.

RESULTS
Key climatic drivers of terrestrial evaporation, such as precipitation, total incoming radiation and near-surface air temperature, are known to be affected by different teleconnection patterns. Therefore, this effect is likely to propagate to the dynamics in land evaporation. 4 Figure 1a shows that monthly evaporation anomalies are significantly (p < 0.05) driven by the main modes of climate variability in about half of the land surface. Clear hotspots, where up to 40% of the variance in land evaporation can be explained by the CIs can be identified for different seasons (Fig. 1b). In these seasonal hotspots, terrestrial evaporation is typically sensitive to a limited number of modes, suggesting both the control of these modes over the local meteorology and the control of the latter over evaporation. Figure 1 shows that throughout the year, variability in terrestrial evaporation in Amazonia is influenced by multiple teleconnections: although ENSO is identified as an important driver of variability in December-February in the east, and in March-May in the south, evaporation dynamics in the rainforest are also sensitive to the Indian Ocean Dipole (IOD) pattern, the Atlantic Multidecadal Oscillation (AMO), the Pacific Decadal Oscillation (PDO), and the Tropical North Atlantic (TNA) pattern. During June-August, the TNA has also a distinct effect in the arc of deforestation (i.e., south of the Brazilian rainforest). Given the high volumes of precipitation in these regions, dynamics in terrestrial evaporation are primarily driven by variations in energy supply, as shown in Fig. 2b (note that while we used total incoming radiation as a proxy for the energy supply, results using surface net radiation show similar patterns as shown in Supplementary Fig. 13). Rainfall interception loss also forms a substantial portion of the total evaporative flux in the rainforest, 29,34 making Amazonian evaporation sensitive to changes in rainfall dynamics as well. The sensitivity of rainfall interception loss to ENSO due to its influence on precipitation variability was already documented by Miralles et al. 4 Figure 2a shows that this mechanism mostly explains the impact of ENSO in the northeast of the rainforest during December-February (see also Supplementary Fig. 5), but not for the arc of deforestation where rainfall dynamics are not affected. In the latter, the impact of ENSO and TNA on evaporation comes from their influence on the energy supply ( Fig. 2a and Supplementary Figs. 10 and 11), mostly due to changes in air temperature and cloud cover. Note that ENSO is also known to induce substantial variability in tropical Atlantic seasurface temperatures, driving the TNA and Tropical South Atlantic (TSA) patterns that in turn alter the meteorology in the northeast of Brazil 35 and 12). Finally, the multidecadal AMO typically modulates the intensity of the ENSO cycle, 37,38 hence the effect of AMO on evaporation in Amazonia is indirect.
In the southwestern Pacific, ENSO, IOD, AMO and the Southern Annular Mode (SAM) influence land evaporation during June-August and September-November (Fig. 1b), especially in the east of Australia, a semiarid region driven by water availability (Fig. 2b). Results in Fig. 2a show that in June-August, precipitation in central east Australia is affected, mainly by SAM and ENSO (see also Supplementary Fig. 7). On the other hand, the hotspot in September-November is due to a combined effect on the water and energy supply induced by ENSO and SAM (Supplementary In contrast to the Southern Hemisphere, the effect of individual modes of climate variability in the Northern Hemisphere is generally concentrated to specific regions. This confined effect is in agreement with the findings by Zhu et al. 15 on the impact of teleconnections on global carbon fluxes. In northwestern Europe, winter time terrestrial evaporation is shown to be sensitive to modes of variability known to modulate the European winter climate: the North Atlantic Oscillation (NAO), and the East Atlantic (EA) pattern. 37,[39][40][41] During December-February, evaporation dynamics in the region are not just influenced by the energy supply ( Fig. 2b) but also by precipitation, due to its influence on rainfall interception. 34 The positive phases of NAO and EA typically bring more precipitation and higher temperatures to northern Europe during winter, 41  The net of dots is presented at a 2°r esolution to aid visibility, and highlights regions with a statistically significant R 2 (p < 0.05). For comparison, the regions of interest were aligned with the ones in Fig. 1. b primary climatic driver of terrestrial evaporation. Classification is obtained by ranking the regression coefficients from the LASSO model targeting the anomalies of terrestrial evaporation based on the anomalies of precipitation and total incoming radiation. Pixels where the regression coefficient for precipitation is top-ranked are classified as water-driven, whereas pixels where the regression coefficient for radiation is top-ranked are classified as energy-driven. NS indicates non-significant pixels (Fig. 2b), and explains the small hotspot along the west coast of Portugal shown in Fig. 1b (not highlighted, but statistically significant). Figure 2a shows that more towards the centre of the European continent, the energy supply during winter is affected, mainly by both NAO and EA ( Supplementary Fig. 9), thereby the influence of these teleconnection patterns on evaporation shown in central Europe.
In eastern Europe and western Russia, the East Atlantic Western Russia (EAWR) pattern starts affecting evaporation dynamics in March-May in the north, and gradually becomes more dominant in summertime (Fig. 1b). Figure 2a shows how the available energy, driving land evaporation during that time of the year (Fig.  2b), is clearly impacted by the EAWR (Supplementary Fig. 11). These findings agree with those by Ionita et al., 42 and indicate that the positive phase of the EAWR typically results in negative temperature anomalies over eastern Europe and western Russia ( Supplementary Figs. 10 and 11), resulting in the negative anomaly in land evaporation observed in Fig. 1b. In the northeast of Russia, terrestrial evaporation is mainly sensitive to the IOD, and the West Pacific (WP) and East Pacific/North Pacific (EPNP) patterns. Both the WP and EPNP induce substantial variability in air temperatures across the region, 1,43,44 resulting in anomalies in the energy supply ( Fig. 2a and Supplementary Figs. 10 and 11), which propagate to the dynamics in terrestrial evaporation.
In North America, terrestrial evaporation responds to different modes of variability, with the EPNP, NAO, WP, the Northern Annular Mode (NAM) and the Pacific-North American (PNA) pattern being the most dominant ones. Figure 1b shows that the EPNP pattern affects evaporation dynamics in the area surrounding the Great Lakes from September to May. The impact of the EPNP on evaporation in North America is opposite to its impact in Siberia ( Supplementary Figs. 1-4). This is due to the contrasting effects of the EPNP on the energy supply in both areas ( Supplementary Figs. 9-12). Figure 1b also shows that the Labrador Peninsula is largely influenced by NAO and NAM from December to May, with negative anomalies in evaporation associated to a positive phase of NAO and NAM (Supplementary Figs. 1 and 2). This effect is driven by the impact of NAO and NAM on the supply of energy in the region, as shown in Supplementary Figs. 9 and 10. NAO is sometimes considered as a regional-scale expression of the large-scale NAM, both having analogous climatological impacts in North America. 39,41,45 The positive phases of NAO and NAM also result in higher temperatures southwards, along with positive radiation anomalies at the East Coast of the United States, 39,45 as illustrated in Supplementary Fig.  9. Despite the strong energy control over terrestrial evaporation in this region (Fig. 2b), the impact of NAO and NAM on evaporation appears lower and is not highlighted in Fig. 1b. The PNA is mainly affecting evaporation in the Canadian Prairies: higher evaporation is associated to the positive phase of the PNA during the winter months (Fig. 1b), due to the impact of PNA on the energy supply in the region (Fig. 2a and Supplementary Fig. 9).
Finally, the hotspot in the southwest of the United States, particularly Mexico, is primarily affected by ENSO, TSA and SAM (Fig. 1b). Given the low precipitation volumes in the area, evaporation dynamics are mainly water-driven (see Fig. 2b). However, only ENSO leads to substantial changes in precipitation, preferentially in winter time ( Supplementary Fig. 5), while the remaining modes act upon evaporation by affecting the atmospheric demand for water ( Supplementary Figs. 9 and 10). During the warm phase of ENSO (i.e., El Niño), the area typically receives more precipitation 37 (Supplementary Fig. 5), resulting in less evaporative stress and positive anomalies in terrestrial evaporation. 4 Figure 2a shows that while precipitation patterns are sensitive to the fundamental modes of climate variability in rather confined areas, total incoming radiation is impacted worldwide. This suggests that most modes affect terrestrial evaporation through altering the atmospheric demand for water. This particularly applies to areas and seasons in which water availability is not limiting the evaporative flux, like northern Europe, Russia, east of North America, and the equatorial belt ( Fig. 2b and Supplementary Fig. 13). Nonetheless, changes in precipitation might still affect the evaporation dynamics in these wet regions through their impact on interception loss, which represents a substantial fraction of the total evaporative flux in forested regions such as Amazonia. 29,34 It should be noted that many of the Northern Hemisphere teleconnections-generally affecting the energy supply-have their primary impact during winter, especially in western Europe, which might in fact make them less important in absolute terms given the lower volumes of evaporation in winter time. In arid and semiarid regions, such as Mexico, the horn of Africa and the east of Australia, the impact of the teleconnection patterns on the energy supply is generally not affecting evaporation dynamics. When existent, the impact of the modes on terrestrial evaporation in these regions is induced by changes in water supply (Fig. 2a). Results presented in Fig. 1b also show that the response of terrestrial evaporation to most modes of climate variability has a short latency, with average lag times ranging between 0 and 2 months.

DISCUSSION
We find that ENSO affects evaporation dynamics in the largest portion of land (Fig. 3, and Supplementary Figs. 14 and 15). This result confirms the notion that ENSO is a major mode of internal climate variability, affecting a wide range of land-surface processes. 4,14,19,25,26,46 However, as shown in Fig. 3, several teleconnections affect terrestrial evaporation over similarly extensive areas. In addition, when several modes influence evaporation concurrently, ENSO is often not the most dominant one. When considering only the dominant mode acting upon evaporation in each pixel, AMO is by far the most important in 0% 4% 6% 4% 8% NH SH 2% 6% 2% Fig. 3 Land area subject to the influence of each teleconnection on terrestrial evaporation. Dark shades show the percentage of land area where each mode dominates the variability in terrestrial evaporation (i.e., the mode is top-ranked in the LASSO model), and light shades show the percentage of land area where each mode affects the variability in terrestrial evaporation (i.e., the oscillation has a regression coefficient > 0 in the LASSO model). Percentages are calculated relative to the total land area, excluding pixels with missing data (e.g., in high latitudes during winter). Importance ranking is based on the magnitude of the regression coefficients from the resulting LASSO model targeting the anomalies of terrestrial evaporation based on the CIs (Supplementary Figs. 1-4 both hemispheres. Nonetheless, given the multidecadal periodicity of AMO, its effect on terrestrial evaporation comes mostly through its influence on other modes of climate variability, such as ENSO. 37,38 Also the TSA and TNA patterns dominate variability in terrestrial evaporation in relatively large regions of both hemispheres, mainly along the equatorial Atlantic. 35 Except for the IOD, which also affects evaporation in both hemispheres, most other teleconnections have their dominant impact in confined areas of the Northern Hemisphere (Fig. 1).
Overall, our results highlight the control of the main modes of climate variability upon the return-flow of water from land to atmosphere. External climate forcing impacts the intensity and frequency of these modes, 12,13,47,48 which will have regionaldependent impacts on local water resources. More efforts to provide observational evidence on the sensitivity of global hydrology to internal climate variability-and to separate its effect from externally induced changes-are a necessary step to benchmark Earth system models, and increase their skill to project the dynamics of our hydrosphere.

Evaporation and climate datasets
Global datasets of terrestrial evaporation, precipitation, net radiation and total incoming radiation are constructed for the period 1982-2012. All datasets are linearly resampled from their native spatial resolution to a common 0.5°latitude-longitude grid, and aggregated to a monthly temporal scale. To avoid spurious results due to similar seasonal patterns, the raw data are converted into monthly, standardised and de-trended anomalies by (a) subtracting the multi-year monthly average for each month of the year, (b) normalising by the standard deviation and (c) subtracting the linear trend at the pixel level. To obtain a robust dataset of terrestrial evaporation in this study, de-trended and standardised anomalies of the GLEAM v3.1a, 28,29 PML 24 and FLUXNET-MTE 30,31 global datasets of terrestrial evaporation are merged based on their relative errors (Supplementary Methods and Supplementary Fig. 16). Total incoming radiation and net radiation are sourced from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim re-analysis dataset, 49 and precipitation from the Multi-Source Weighted-Ensemble Precipitation (MSWEP) 50 dataset. Results are thus affected by the uncertainties of the ERA-Interim re-analysis and the MSWEP datasets. Nevertheless, both have been validated against in situ data, and compared against similar datasets in numerous studies, 50,51 and it has been concluded that both datasets are robust and reliable in terms of intraannual variability.
The sixteen selected modes of climate variability (Table 1) are diagnosed by CIs, generally obtained by temporal analysis of sea-surface temperature or geopotential height fields in specific regions or points on Earth. These indices provide a simple statistical summary of complex and inter-related processes resulting in the modes of internal climate variability. 39 CIs are obtained from various centres from the National Oceanic and Atmospheric Administration (NOAA) and the Climate Explorer (ClimEx) of the Royal Meteorological Institute of the Netherlands (KNMI).

Approach
First, the impact of sixteen modes of climate variability-diagnosed by their CIs-on land evaporation is analysed by targeting the dataset of terrestrial evaporation, using the CIs as predictive features in a LASSO framework. To account for the unknown response time of terrestrial evaporation to variations in the climate, time lags ranging from 0 (i.e., no lag) to 8 months are introduced to each of the CIs, resulting in 144 predictive features in total (i.e., 16 CIs x 9 lags/CI). Acknowledging the dependency among CIs ( Supplementary Fig. 17), resulting from interactions between different climate phenomena impacting each of the CIs, 13,33 the LASSO is selected here as a regression method, as it has proven to perform well when features are cross-correlated. 32 The LASSO has the unique property to perform both feature selection and regularization at the same time, typically resulting in a predictive model with only a subset of the selected features, which simplifies the interpretation of the final model. The method shrinks the regression coefficients towards zero by aiming at finding coefficients that minimize a penalized residual sum of squares: 32 , whereβ is the p-dimensional vector with the estimated regression coefficients, n the number of training samples in the dataset, y i the value of the target variable in sample i, p the number of features, x ij the value of feature j in sample i and λ a hyper-parameter controlling the amount of shrinkage. The latter is tuned by minimizing the prediction error of the model, estimated using a fivefold cross-validation. Second, the relative control over terrestrial evaporation of the two local climatic drivers (i.e., precipitation and total incoming radiation) is quantified by fitting a LASSO model for each pixel targeting the time series of terrestrial evaporation, and using precipitation and radiation as predictive features. Note that the total incoming radiation is used here as a proxy for the energy supply, as it is to a large extent not affected by the land surface, and thus almost independent from terrestrial evaporation itself (except through its indirect effect on temperature, humidity and clouds 11 ). This is in contrast to other studies which typically use net radiation-or potential evaporation-to define evaporation regimes. Supplementary Fig. 13, however, shows that the evaporation regimes derived using net radiation are similar as the ones discussed in Fig. 2b, and are comparable with maps published in literature. 8 Third, the LASSO models are applied to target precipitation and radiation based on the same 144 predictive features used in the first step (i.e., the lagged CIs). While the first step identifies the key modes of climate variability that impact evaporation for a certain pixel, the second step allows to hypothesise through which meteorological driver these modes may be acting upon. The third step confirms whether the mode of variability actually impacts the expected climatic driver. Altogether, we identify not just the dominant teleconnections controlling land evaporation, but also the local meteorological drivers through which these modes act upon evaporation.

Model evaluation
The true prediction error of each fitted model (including the hyperparameter tuning, hence feature selection procedure) is estimated using an external fivefold cross-validation on top of the cross-validation used for hyper-parameter tuning. This results in a nested cross-validation that avoids biases in the estimates of the true prediction error due to the internal feature selection procedure, 52 and avoids the use of an independent test set for model evaluation. The explained variance R 2 of the final LASSO model is then calculated as: , where RSS is the residual sum of squares calculated using the nested fivefold cross-validation, and TSS is the total sum of squares. R 2 values are tested for statistical significance (p < 0.05) by comparing against a null hypothesis of no predictive power. The null distribution of the R 2 of a LASSO model without any predictive power is obtained using a nonparametric permutation test. 53 Therefore, all LASSO models are fitted on a random training dataset, obtained by shuffling the original predictive features in time. The null distribution is then nonparametrically described by the obtained R 2 values, and compared with the anologue of the original experiment. To control the false discovery rate at a level of 5% in this multiple hypothesis testing problem, the obtained p-values were adjusted according to the Benjamini-Hochberg procedure. 54 Using this method, the critical R 2 value to declare significance in case of the experiment-targeting evaporation using the CIs, ranges between 4 and 6%, depending on the season. Similar values are obtained for the other experiments.

Feature ranking
Pixel-based ranking of features (i.e., either the climatic oscillations or the drivers of evaporation) is based on the magnitude of the regression coefficients in the fitted LASSO model targeting the variable of interest (i.e., either evaporation, precipitation, or total incoming radiation). To account for the potential effect of climate modes on the variable of interest over different time lags, regression coefficients were first summed overall lags per CI to obtain Supplementary Figs. 1-12. An area-based rank of the climatic oscillations in Fig. 1 is then obtained by calculating the average pixel-based rank per CI, over the area of interest.

Code availability
The computer codes have been implemented in open source software, and are available through the GitHub channel of the Laboratory of Hydrology and Water Management -Ghent University: https://github.com/ lhwm.