Land use change impacts on European heat and drought: remote land-atmosphere feedbacks mitigated locally by shallow groundwater

Heat and drought are projected to increase globally but may be mitigated or exacerbated by land use/land cover (LULC) change. Here, we show that remote land-atmosphere feedbacks caused by historical European LULC change led to widespread changes in the energy and water balances, drought, and heat. Using a continental-scale bedrock-to-atmosphere model, we find that LULC change following the Soviet Union collapse and European Union formation may have substantially increased cloud cover and decreased incoming shortwave radiation in western Europe, even in locations where LULC did not change. These changes to the water and energy balances had spatially heterogeneous impacts on drought and heat, including drying in the Mediterranean and Eastern Europe regions. The response of the water and energy balances to remote feedbacks was lessened in areas with shallow groundwater, indicating that local- and continental-scale responses to LULC change are influenced by the coupling between the subsurface, land surface, and atmosphere.


Introduction
Land use/land cover (LULC) change is a key anthropogenic driver of global change. LULC change, in particular shifts from natural vegetation to agriculture, can have significant impacts on local, regional, and global climate by changing the water and energy balance at the land surface (Lejeune et al 2016, Han et al 2017, Perugini et al 2017, Mohan et al 2018. For example, recent work has demonstrated decreased diurnal temperature range (Hua and Chen 2013), reduced evapotranspiration (Sterling et al 2013), and changes to regional teleconnections (Wang-Erlandsson et al 2018) caused by LULC change. However, the impacts of LULC change are diverse due to variability in transition pathways and differences in local management practices (Luyssaert et al 2014, Mahmood et al 2014.
Feedbacks between LULC and atmospheric circulation can directly impact human health and livelihoods by altering the frequency, severity, or duration of extreme weather, e.g. heat or drought. One recent example of heat and drought occurred across Europe in 2003. The Europe 2003 event was characterized by a warm early growing season combined with a pronounced precipitation deficit, leading to low soil moisture throughout Europe and a catastrophic heat wave in early August which caused more than 40 000 deaths and 13.1 billion Euros in damages (Fink et al 2004, García-Herrera et al 2010. Previous research found that the Europe 2003 event was exacerbated by a strong coupling between the subsurface (soil moisture), land surface (evapotranspiration), and atmosphere (boundary layer processes and regional circulation), and that such events are projected to increase under future warming scenarios (Samaniego et al 2018). However, the effects of LULC change are less clear. Stefanon et al (2014b) suggest that agriculture may have exacerbated late summer heat due to drought-induced crop failure. There is some evidence that historical LULC change within Europe has led to a warmer and drier summer climate (Heck et al 2001) and changes in the frequency of extreme heat (Anav et al 2010). However, other research indicates that European LULC change is contributing to cooler summer temperatures (Zampieri and Lionello 2011), or that effects of LULC change are only important in the spring (Ma et al 2016).
While previous work has shown that groundwater can strongly influence land surface processes (Zipper et al 2015, Maxwell and Condon 2016, Chang et al 2018, little is known about how groundwater influences the response of land surface and atmospheric processes to LULC change at continental scales. During the Europe 2003 event, Keune et al (2016) found that groundwater representation had a large impact on subsurface processes, an intermediate effect on land surface processes, and a small effect on atmospheric processes. However, the role of groundwater in influencing the response of ecosystems to LULC change remains poorly understood and confined to relatively local case studies (Condon and Maxwell 2014, Giménez et al 2016, Zipper et al 2017.
Thus, the impacts of LULC change on heat and drought remains highly uncertain, both during the Europe 2003 event and more generally. In this study, we use a fully-coupled 'bedrock-to-atmosphere' model to explore the sensitivity of heat and drought to LULC change in a physically-and spatially-consistent manner over all of Europe, and additionally explore the role of groundwater in mitigating land surface and atmospheric response to LULC change. We develop two new LULC datasets describing historically observed LULC change between 1990 and 2010 associated with two major geopolitical changes: the formation of the European Union and the collapse of the Soviet Union (Prishchepov et al 2012, Kuemmerle et al 2016. Specifically, we ask: (1) How sensitive are the terrestrial water and energy balances to historical European LULC change?
(2) What impact does LULC change have on heat and drought in Europe, and how do impacts vary across spatial scales (grid cell, regional, continental)?
(3) To what degree does shallow groundwater mitigate or enhance the response of the land surface and atmosphere to LULC change?

Methods
We explore these research questions using an ensemble of continental-scale simulations with the Terrestrial Systems Modeling Platform (TerrSysMP), a process-based modeling platform simulating the complete water and energy cycles (Gasper et al 2014. In the following sections we briefly describe TerrSysMP and the numerical experiments conducted in this study.

Modeling platform
TerrSysMP couples three models to close the terrestrial water and energy cycles from groundwater across the land surface into the atmosphere: (1) ParFlow (v3.1), a surface-subsurface flow model (Ashby andFalgout 1996, Kollet andMaxwell 2006) TerrSysMP and its component models are described in the above-referenced publications, therefore we provide only a brief description. Coupling between models is centered around CLM. From below, CLM is coupled to ParFlow, which simulates three-dimensional variably-saturated subsurface flow and free-surface overland flow. ParFlow and CLM are coupled via the exchange of moisture occurring at the land surface and in the root zone. From above, CLM is coupled to COSMO, which simulates atmospheric circulation and meteorological conditions, via the exchange of moisture and energy at the top of the canopy. In this study, both ParFlow and CLM have a timestep of 180 s, COSMO has a timestep of 60 s, and the models are coupled at a frequency of 180 s to accurately simulate non-linear feedbacks in the exchange processes between the subsurface, surface, and atmosphere.

Model domain, discretization, and parameterization
We applied the model setup and domain described in Keune et al (2016Keune et al ( , 2018. Spatially, this corresponds to the EUROpean COordinated Regional Downscaling EXperiment (EURO-CORDEX) domain at 0.11°resolution (Giorgi et al 2009). We used ERA-Interim reanalysis data (Dee et al 2011) as lateral boundary conditions for the COSMO circulation model. To keep the large-scale circulation of the LULC simulations consistent with ERA-Interim, we performed spectral nudging (von Storch et al 2000) on horizontal winds at wavelengths >200 km above the planetary boundary layer using a nudging coefficient of 0.05. Given the relatively coarse grid applied in the continental-scale simulations, spectral nudging provides many of the benefits of data assimilation at a fraction of the computational cost (Schubert-Frisius et al 2016). In COSMO, convection is parameterized using the Tiedtke mass flux scheme (Tiedtke 1989), which is widely used in large-scale models which do not explicitly resolve convection processes (Lejeune et al 2015, Keuler et al 2016, Prein andGobiet 2017). Hydrostratigraphic properties for the top 3 m were defined based on the UN Food and Agriculture database (UN FAO 1988), and below 3 based on the GLHyMPS database (Gleeson et al 2014). LULC data are described below. Within EURO-CORDEX, we analyzed results from a smaller focus region which excludes areas unrelated to our research questions (e.g. North Africa) and reduces potential edge effects. Simulating the entire EURO-CORDEX domain allows us to evaluate local-and regional-scale response to LULC change within a consistent large-scale integrated modeling framework.
Temporally, the simulations span 1 January-18 October 2003. To ensure that the model was in dynamic equilibrium, we first spun up CLM for 50 years with repeated 2002 meteorological forcing data, followed by an additional 65 month spin-up of the coupled ParFlow-CLM model. Results from the spinup simulations were used as initial conditions for the fully coupled TerrSysMP runs. Spin-ups were performed separately for the 1990 and 2010 LULC input datasets.

Numerical experiments
Our numerical experiments consisted of a six-simulation ensemble, allowing us to average internal model variability and identify mean responses to LULC. The ensemble comprises two land cover datasets corresponding to 1990 and 2010, with three different vegetation parameterizations accompanying each land cover dataset.

Creating LULC input datasets
To evaluate the response of the terrestrial water/ energy balances, heat, and drought to historically realistic LULC change over the European continent, we developed new land surface input datasets corresponding to the years 1990 and 2010 (figure 1). Due to the mismatch between meteorological forcing (2003) and LULC scenarios (1990 and 2010), our simulations were not intended to reproduce real-world conditions, but rather to quantify the sensitivity of heat and drought to historically observed LULC change. We selected 1990 and 2010 as LULC scenarios because they represent a before/after condition for major geopolitical events which caused significant LULC change: the breakup of the Soviet Union in 1991 and the formation of the European Union in 1992.
We elected to create new datasets because existing continental-scale LULC datasets struggle to capture some LULC change trajectories relevant to our study, specifically agricultural abandonment (Kuemmerle et al 2016). Creation of the new LULC input datasets is described in detail in the supplementary information (available online at stacks.iop.org/ERL/14/044012/ mmedia). Our new LULC input datasets use a hybrid of MODIS and HILDA land cover datasets, which are described in Friedl et al (2010) and Fuchs et al (2013Fuchs et al ( , 2015, respectively. This approach takes advantage of the high spatial resolution and plant functional type (PFT) discrimination of MODIS and the countrylevel ancillary LULC statistics used to constrain HILDA, which has fewer PFTs and does not cover the entire EURO-CORDEX domain. In our datasets, LULC differs between 1990 and 2010 for the 28 countries in the European Union+Switzerland domain covered by HILDA, but not for former Soviet regions within our domain such as Ukraine, Belarus, Western Russia, and the Balkans, in addition to Norway. Given that local agricultural abandonment rates ranging from 13% to 56% have been estimated for these regions ( In order to estimate uncertainty associated with internal model variability in the groundwater-toatmosphere simulations, we created an ensemble set of leaf area index (LAI) and stem area index (SAI) scenarios for each of the 1990 and 2010 land cover datasets based on the GLASS LAI product (Xiao et al 2014). For each CLM grid cell, we calculated the mean and standard deviation (stdev) of all LAI and SAI values for each PFT by overlaying our land cover map onto the maps of LAI and SAI. We then used the mean, the mean+stdev, and the mean-stdev as our ensemble set of 3 land cover scenarios. This procedure was repeated separately for 1990 and 2010 to create a total of 6 continental-scale LULC datasets.

Model output and analysis
TerrSysMP provides hourly output for a suite of variables spanning the bedrock-to-atmosphere domain. The primary variables related to our research questions, examples of which are shown in figure S3, quantify the water and energy budgets and include incoming shortwave radiation (SWin), latent heat flux (LE), peak hourly 2 m air temperature (T peak ), precipitation deficit (PrecDef, calculated as evapotranspiration [ET]-precipitation [P]), soil water storage (SWS) and water table depth (WTD).
We analyzed model output at multiple spatial scales. To investigate continental-scale responses to LULC change, we developed maps showing the difference between simulated results with 1990 and 2010 LULC. To investigate different geographic regions, we summarized changes from all grid cells within PRU-DENCE regions (figure 1). To investigate specific types of LULC change, we subsampled our domain to focus on seven common LULC trajectories: 1990 Crop to 2010 Grass (C-G), Grass to Crop (G-C), Crop to Forest (C-F), Forest to Crop (F-C), Crop to Crop (C-C), Grass to Grass (G-G), and Forest to Forest (F-F). The trajectories with constant LULC between 1990 and 2010 (C-C, G-G, and F-F) are intended to provide a control for comparison with the four trajectories in which LULC changes. To minimize potential confounding effects of different locations between these seven trajectories, we first identified all points in the C-G, G-C, and C-F categories, then identified the C-C cell closest to each C-G point; and the G-G cell closest to each G-C point; and the F-F cell closest to each C-F point. Except where otherwise noted, all results shown are the mean of our ensemble (3 LAI scenarios for each LULC dataset). We did not conduct a validation of our model due to our experimental design, which combined LULC from years 1990 and 2010 (to investigate changes in LULC associated with important sociopolitical events) with atmospheric forcing from the year 2003 (to investigate extreme climatic conditions). However, Terr-SysMP was previously evaluated over the same domain (EURO-CORDEX) with the same discretization (0.11°) and the same atmospheric forcing data (2003 ERA-Interim) (Keune et al 2018). To briefly summarize, Keune et al (2018) compared simulated precipitation to measurements from 1246 stations and simulated evapotranspiration with measurements from 28 eddy covariance towers. TerrSysMP had a slight wet bias (~0.2-2 mm d −1 ) for annual precipitation, which is smallest in the France, Mid-Europe, and Eastern Europe regions. During the summer months, the model bias was not consistent in sign across regions but was also consistently smaller (±1 mm) in Eastern Europe, Mid-Europe, and France. The mean evapotranspiration bias, which is representative of the model's skill at flux partitioning at the land surface, was within ±1 mm d −1 across 5 of 6 land cover types tested including multiple types of forest, grasslands, and cropland at summer and annual scales. These biases are comparable in magnitude to other regional climate simulations (Kotlarski et al 2014, Katragkou et al 2015). More broadly, a recent comparison of 9 regional climate models over the EURO-CORDEX domain found that coupled COSMO-CLM had among the smallest biases for simulated precipitation, temperature, and radiation compared to observations and reanalysis data (Davin et al 2016). As a state-of-the-art modeling platform, additional model evaluation is ongoing (Keune et al in review).

Results and discussion
Effects of LULC change on the energy balance Our new LULC datasets document substantial LULC change between 1990 and 2010, described in the supplementary information. Net LULC change is characterized by a reduction in cropland (agricultural abandonment) which is balanced by an increase in other LULCs, primarily grassland (figures 1, and S1). Agricultural abandonment was concentrated in the Eastern Europe (EA), Mediterranean (MD), Iberian Peninsula (IP), and Scandinavian (SC) regions (figures 1 and S2). At a regional scale, there are significant differences in the partitioning of incoming energy into latent and sensible heat between LULC scenarios, which follow a latitudinal gradient (figures 2 and S3). Specifically, changes in LULC led to statistically significant increases in mean summer latent heat fluxes and decreases in sensible heat fluxes in the southern portions of the domain, most strongly the Mediterranean, Eastern Europe, and the Iberian Peninsula regions. The latent and sensible heat changes in the Mediterranean and Eastern Europe regions are dominated by changes to the north and west of the Black Sea ( figure 1(a)).
We also observe continental-scale shifts in summer incoming shortwave and longwave radiation (figure 2) and cloud cover (figure 3) associated with LULC change. Over the western portion of the domain, incoming shortwave radiation decreases while cloud cover and incoming longwave radiation both increase. Mean changes in incoming longwave radiation at the regional scale are approximately −45% to +57% of the simulated changes in incoming shortwave radiation. The Mediterranean is the only region which has reductions in both incoming shortwave and longwave radiation; combined, these changes in incoming radiation are approximately counterbalanced by changes in the turbulent fluxes (latent and sensible heat). At the grid-cell scale, cloud cover is strongly correlated with differences in incoming longwave radiation (R 2 =0.61) and incoming shortwave radiation (R 2 =0.48).
Our results indicate that LULC-driven changes in cloud cover and the summer energy balance may be widespread over most of Europe and occur in regions of both agricultural expansion and agricultural abandonment. The observed changes in cloud cover are temporally extensive, occurring through the entire summer (figure S5), and robust across our ensemble of LAI scenarios (figure S6). In contrast to previous work in Australia (Ray et al 2003, Nair et al 2011, our simulations indicate that decreases in cloud cover are more prevalent in areas with substantial agricultural abandonment such as Eastern Europe (figures 1(e) and 3(a)), which may also contradict observations showing enhanced cloud cover over forests relative to agricultural land in France (Teuling et al 2017). Our simulations differed from these prior studies in that abandoned agricultural land was predominantly converted to grassland rather than forest (figure 1). However, the simulated changes in cloud cover in our experiments are likely driven in part by both local and regional-scale changes in energy partitioning, and may be related to convection parameterization (Hohenegger et al 2009, Prein et al 2016. For instance, some areas (e.g. much of Mid-Europe; figures 2(c) and (f)) experience significant changes in energy partitioning, with an increase in latent heat fluxes and decreases in sensible heat fluxes, but an increase in cloud cover ( figure 3(a)), which may contradict recent case studies with a convection-resolving model (Bosman et al 2019) and points to the importance of finer-resolution future simulations to resolve potential issues associated with convection parameterizations (Prein et al 2016. The importance of remote land-atmosphere feedbacks is further supported by comparison among LULC trajectories for various components of the energy balance. For the turbulent fluxes (latent and sensible heat) we found that all LULC trajectories resulted in increases in latent heat and incoming longwave radiation and decreases in sensible heat and incoming shortwave radiation, with few significant differences between LULC trajectories (figures 2(a), (d), (h), (k)). The lack of differentiation between LULC trajectories indicates that, while there are changes in water and energy partitioning at the grid-cell level caused by LULC change, they cannot be distinguished from larger atmospheric-scale changes in incoming shortwave and longwave radiation associated with local and remote land-atmosphere feedbacks and driven primarily by cloud cover.

Effects of LULC change on heat and drought
There is also substantial spatial heterogeneity in the response of heat to LULC change (figures 4(a)-(c)). , the horizontal line shows the median value for that LULC trajectory. Distributions with the same letter do not significantly differ from each other (p>0.05 in Tukey Honest Significant Differences test), and the caret symbol (^) indicates that the mean of that sample does not differ significantly different from 0 (p>0.05 using one-sample t-test). In maps, colorbars reverse for some variables, such that red always corresponds to hotter/drier conditions. Over the entire domain, differences in peak 2 m air temperature between LULC scenarios range between −1.2°C and +1.0°C at the grid cell scale (1st and 99th percentiles) with significant decreases at the grid cell scale over much of the southern half of Europe ( figure 4(c)). While the mean overall difference between LULC scenarios is small (−0.09°C), the mean absolute change is larger (0.31°C) and changes exceed ±0.5°C in 19.9% of terrestrial grid cells. There are few significant differences in heat among LULC trajectories ( figure 4(a)), but significant regional-scale decreases in peak temperature in 7 of 9 PRUDENCE regions, with the largest mean changes in the France (−0.30°C), Mediterranean (−0.26°C), and Eastern Europe (−0.20°C) regions; peak temperature increases by a small but statistically significant amount in the British Isles (+0.02°C) and Scandinavia (+0.12°C).
These changes in peak temperature are driven by LULC-induced remote land-atmosphere feedbacks leading to differences in available energy. Changes in regional mean sensible heat between LULC scenarios are strongly correlated with the observed variability in peak temperature differences (R 2 =0.64). Our results suggest that increased cloud cover led to decreased available energy and lower sensible heat fluxes under 2010 LULC, which in turn reduced regionally-averaged peak summer temperatures. In turn, changes in regional mean sensible heat are strongly associated with differences in SWS and incoming shortwave radiation (multiple linear R 2 =0.90), suggesting that changes in both energy availability and water availability associated with LULC change drive shifts in the distribution and magnitude of heat.
Changes in total summer precipitation (figures 4(d)-(f)) exhibit less spatial coherence than the observed changes in heat. At the grid-cell scale, differences in precipitation between LULC scenarios can exceed ±160 mm but are rarely statistically significant ( figure 4(f)). At regional scales positive and negative changes cancel each other out to a large degree, with mean changes smaller than ±3.5 mm due to close spatial proximity between increases and decreases in precipitation ( figure 4(f)). This is most likely a result of small-scale shifts in precipitation in the simulations (e.g. a storm track shifting by a few kilometers would lead to neighboring increases and decreases in precipitation) and is consistent with high spatial variability in simulated precipitation changes under potential future land cover trajectories (Cherubini et al 2018). Similarly, summer precipitation deficit (evapotranspiration-precipitation), which integrates water consumption and supply as an indicator of drought, varies substantially at the grid-cell level primarily associated with changes in precipitation (R 2 =0.79) and to a lesser degree changes in ET (R 2 =0.11). Changes in summer SWS, an indicator of plant available water, are also strongly associated with the observed differences in precipitation (figure 4(e)) but do not differ significantly among our paired groups of LULC trajectories (figures 3(k) and 4(h)).
From a temporal perspective, there are several distinct pathways by which LULC change either exacerbates or reduces drought through changes of plant water use and supply (figure 5). In some regions, changes in precipitation and ET partially or wholly counteract each other. For example, the Mediterranean region has an increase in both precipitation and ET with 2010 LULC relative to 1990 LULC, and because the increase in ET is greater than the increase in precipitation, the net change in precipitation deficit and SWS are both more drier with 2010 LULC. In contrast, in the France region, increases in ET are approximately balanced by increases in precipitation of a similar magnitude, while in the Scandinavian region both precipitation and ET decrease by a comparable magnitude. In other regions, changes in ET and precipitation both exacerbate drought conditions. For example, the increases in ET in the Mediterranean are approximately double those observed in Eastern Europe, but coincident precipitation decreases in Eastern Europe leading to a net change in precipitation deficit at the end of the growing season which is comparable to the Mediterranean regions.
Synthesizing local and remote sensitivity to LULC change The results discussed above indicate that LULC change can lead to significant regional-to continental-scale changes in available energy (e.g. incoming radiation,  (k), the horizontal line shows the median value for that LULC trajectory. Distributions with the same letter do not significantly differ from each other (p>0.05 in Tukey Honest Significant Differences test), and the caret symbol (^) indicates that the mean of that sample does not differ significantly different from 0 (p > 0.05 using one-sample t-test). In maps, colorbars reverse for some variables such that red always corresponds to hotter/drier conditions.
figure 2), primarily by changing atmospheric moisture transport, small scale circulation, and cloud formation (figure 3). These atmospheric, 'top-down' changes in available energy due to LULC change overwhelm potential land surface, 'bottom-up' grid-cell-level differences in water and energy partitioning (e.g. latent versus sensible heat, figure 2) between PFTs. Thus, the effects of LULC change on heat and drought is not confined to locations where LULC is changing but propagate to surrounding areas via remote land-atmosphere feedbacks largely driven by changing cloud cover.
We also found that widespread and significant reductions in extreme heat associated with historic LULC change which may have been driven by the changes in energy partitioning and cloud cover discussed above. Peak summer temperatures were significantly lower with 2010 LULC than 1990 LULC in 7 of 9 regions, with the largest reduction in France-the region which suffered most severely during the Europe 2003 event (Zaitchik et al 2006). Previous studies of European LULC change impacts have primarily compared modern LULC to preindustrial vegetation as a baseline (Heck et al 2001, Findell et al 2007, Christidis et al 2013, Stefanon et al 2014b, Lejeune et al 2016. While useful for understanding the sensitivity of landatmosphere feedbacks to LULC change, comparisons to natural vegetation cannot quantify how modern LULC change may impact heat and drought. Our results, which are likely to be conservative estimates of LULC change, indicate that ongoing regional LULC changes occurring over decadal timescales have the potential to significantly influence drought and heat at continental scales. This is further evidence that LULC may be part of a portfolio of strategies to combat negative effects of climate change at regional scales (Cherubini et al 2018, Seneviratne et al 2018.
While our modeling approach is physically-based, it is nevertheless subject to uncertainty, particularly since the focus of our study is extreme conditions which are challenging to simulate (Lhotka et al 2018). Model uncertainties are associated with (i) model input data, e.g. the accuracy of the input LULC data (Fuchs et al 2015), LULC parameterization , and hydrostratigraphic data (Gleeson et al 2014); (ii) the structure of the model, e.g. the simulated timestep and grid resolution entailed parameterization of some exchange processes (Hohenegger et al 2009; and (iii) the mathematical representation of processes within the model, e.g. alternate parameterizations for land surface-atmosphere coupling or atmospheric convection (Davin et al 2016. Options for reducing these uncertainties, such as using observational data to constrain simulations (Stegehuis et al 2013), data assimilation (Zhang et al 2018, Naz et al 2019, and finer model discretization  are becoming more feasible as computational power continues to increase and may provide a path towards reduced uncertainty in future work.

Influence of shallow groundwater
Soil moisture-atmosphere coupling has been shown to influence the onset and exacerbation of heat and drought in Europe, including the 2003 event (Ferranti and Viterbo 2006, Fennessy and Kinter 2011, Stefanon et al 2014a. Our results extend this coupling further into the subsurface by identifying groundwater as a potential stabilizing mechanism for the local response of the water/energy budgets, heat, and drought to LULC change. Across our focus domain, the water table is within 4 m of the land surface over~79% of our domain ( figure S3). Comparing among proximal locations with different LULC trajectories, we find that changes in peak temperature, precipitation deficit, and SWS are not significantly affected by LULC trajectory when the WTD is very shallow (WTD<1 m) but change significantly when WTD>1 m (figure S8). These results indicate that shallow groundwater may reduce the sensitivity of the water and energy balance, as well as associated heat and drought indicators, to both local and remote LULC change.
The reduced sensitivity to LULC change in areas with shallow groundwater is consistent with recent studies demonstrating that groundwater can significantly impact water and energy transport across the critical zone (Miguez-Macho and Fan 2012, Keune et al 2016, Maxwell and Condon 2016, Zeng et al 2018. In areas where the water table is shallow, SWS is high (figure 6(a)) and not affected by LULC change (figure S8) regardless of potential changes to the energy balance resulting from local LULC or remote feedbacks (figures 2 and 4). This may stabilize the portion of incoming energy going into the latent heat flux (evaporative fraction), which is positively correlated with SWS at the grid cell scale (figures 6(b) and S7). The evaporative fraction shows a complex relationship with heat ( figure 6(c)). In water-limited settings, reductions in latent heat (lower evaporative fractions) are associated with increases in peak temperature because there is less evaporative cooling and a higher sensible heat flux. In contrast, in energy-limited environments where virtually all incoming shortwave radiation is going to latent heat, warmer periods are associated with higher available energy and greater latent heat rates.
These results suggest that groundwater is both locally and regionally important as a driver of the land surface and atmospheric response to LULC change which has been neglected by previous continentalscale research. For instance, inaccessible groundwater may be an important factor contributing to the observed increase in the frequency and duration of heat waves in southern Europe (Fischer and Schär 2010), which is the region of Europe over which the water table is deepest ( figure S3). Given the important role of groundwater in controlling land surface processes (Kollet and Maxwell 2008, Lowry and Loheide 2010, Hain et al 2015, Booth et al 2016, particularly during dry periods (Zipper et al 2015, Fang et al 2017, Taufik et al 2017, this points to an important need to explicitly include groundwater in LULC change simulations. While groundwater impacts on land surface processes are often thought of as a local process occurring at relatively small spatial scales (e.g. in zones of groundwater convergence), we demonstrate that explicit representation of groundwater in a process-based manner can significantly impact simulated land surface and atmospheric response to LULC change even at large spatial scales (Fan 2015).

Conclusions
In this study, we evaluated the sensitivity of the energy/water balances, heat, and drought to historically observed LULC change over Europe. We harmonized multiple LULC datasets to document spatially heterogeneous LULC change between 1990 and 2010, with extensive agricultural abandonment occurring in Eastern Europe, the Mediterranean, and Iberian Peninsula regions. Using a fully coupled bedrock-toatmosphere model we found that: (1) The impacts of LULC change on the water and energy balances were not confined only to areas undergoing LULC change, but may also lead to regionally significant changes in drought and heat; (2) LULC change between 1990 and 2010 decreased simulated peak air temperature over much of the domain, including the region hardest hit by the Europe 2003 event (e.g. France), which was associated with decreases in available energy caused by increased cloud cover; and, (3) groundwater near the land surface (WTD<4 m, especially WTD<1 m) plays a potentially important role in reducing sensitivity to both local LULC change and remote land-atmosphere feedbacks by maintaining the soil water supply.
Overall, we conclude that historical European LULC change which took place over a relatively short duration (20 years) likely had an impact on the water and energy balances at regional to continental scales.

Acknowledgments
SCZ thanks the Bundesministerium für Bildung und Forschung (BMBF) 'Green Talents' program for funding his research stay in Germany. JK acknowledges support from the European Research Council (ERC) under grant agreement no. 715254 (DRY-2-DRY). TerrSysMP simulations were run on the European Center for Medium-range Weather Forecasts (ECMWF) High Performance Computing Facility. We