Historical land-use-induced evapotranspiration changes estimated from present-day observations and reconstructed land-cover maps

. Recent results from the LUCID model intercom-parison project have revealed large discrepancies in the simulated evapotranspiration (ET) response to the historical land-use change. Distinct land-surface parameterizations are behind those discrepancies, but understanding those differences rely on evaluations using still very limited measurements. Model benchmarking studies with observed ET are required in order to reduce the current uncertainties in the impacts of land use in terrestrial water ﬂows. Here we present new estimates of historical land-use-induced ET changes based on three observation-driven products of ET. These products are used to derive empirical models of ET as a function of land-cover properties and environmental variables. An ensemble of reconstructions of past ET changes are derived with the same set of land-cover maps used in LUCID, with which we obtain an average decrease in global terrestrial ET of 1260 ± 850 km 3 yr − 1 between the preindustrial period and the present-day. This estimate is larger in magnitude than the mean ET change simulated within LUCID with process-based models, and substantially weaker than other estimates based on observations. Although decreases in annual ET dominate in deforested regions, large summertime increases in ET are diagnosed over areas of large cropland expansion. The multiple ET reconstructions carried out here show a large spread that we attribute principally to the different land-cover maps adopted and to the crops’ ET rates deduced from the various products assessed. We therefore conclude that the current uncertainties of past ET changes could be reduced efﬁciently with improved historical land-cover reconstructions and better estimates of cropland ET.


Introduction
Land-use-induced land-cover change (LULCC) has been one of the major environmental changes driven by human activities.During the last 300 years, the large-scale deforestation that occurred in the northern temperate regions has significantly contributed to the rise in concentration of atmospheric carbon dioxide and to the underlying global temperature increase (e.g., Pongratz and Caldeira, 2012;Ciais et al., 2013).
In addition to the biogeochemical impact of LULCC on climate, more direct and regionally important perturbations are driven by changes in the physical properties of the surface (biogeophysical effects).These changes are often difficult to characterize because of the multiple mechanisms involved (Davin and de Noblet-Ducoudré, 2010).Local cooling driven by an increase in surface albedo, and warming, due to reduced evaporative cooling, are two possible effects of deforestation of opposite sign (Bonan, 2008).The latter effect may dominate in the tropics, as several modelling (e.g., Nobre et al., 1991;Costa and Foley, 2000;Sampaio et al., 2007;Brovkin et al., 2009;Davin and de Noblet-Ducoudré, 2010) and observational (Gash and Nobre, 1997;Von Randow et al., 2004;Da Rocha et al., 2009;Loarie et al., 2011) studies have shown.In contrast, because of the strong snow masking effect exerted by the forest canopy, the radiative (albedo) impact of forest clearing has likely led to surface cooling at high latitudes (e.g., Betts, 2001;Govindasamy et al., 2001;Bounoua et al., 2002;Brovkin et al., 2009).
Published by Copernicus Publications on behalf of the European Geosciences Union.

Boisier et al.: Land-use-induced evapotranspiration changes estimated from observations
The climate impact of past LULCC is particularly uncertain in temperate regions, in part because of the unknown net effect of the above-mentioned radiative and non-radiative effects of deforestation, but also because the direction of change in evapotranspiration is clearly not one-sided (Sterling et al., 2013).
Evapotranspiration (ET) is a key variable of the climate system, as it affects both the energy and the water balance of the surface.Changes in ET due to LULCC or other landuse practices, such as irrigation, have received special attention over the past decade because of their potential effects on climate and water resources.Although most studies suggest that historical LULCC has led to a decrease in global ET (e.g., Gordon et al., 2005;Findell et al., 2007;Sterling et al., 2013) and consequent increase in runoff (Piao et al., 2007), the large-scale changes in ET remain quite uncertain, as do the geographical and seasonal variations of such changes (Pitman et al., 2009;Pielke et al., 2011).
Studies at regional scales have shown that ET increases resulting from widespread irrigation have induced surface cooling and other impacts on climate in India (Douglas et al., 2006;Roy et al., 2007;Guimberteau et al., 2012), in the Middle East and Asia (Lee et al., 2011), in the North American Great Plains (Adegoke et al., 2007;Mahmood et al., 2006) and in California (Lobell and Bonfils, 2008), among others.Puma and Cook (2010) have suggested that large-scale ET increases and cooling induced by irrigation may have been as large in magnitude as the opposite effect driven by deforestation, a finding consistent with the results of Gordon et al. (2005) and Haddeland et al. (2007).Observation-based studies have also shown that well-watered cropland can evaporate more than temperate forest (Baldocchi et al., 1997;Teuling et al., 2010;Sterling et al., 2013).
The LUCID project ("Land Use and Climate: Identification of Robust Impacts") has compared outputs from different climate models, each forced with the same historical change in crop and pasture area.Recent results have revealed very large uncertainties in the simulated ET responses to LULCC (Pitman et al., 2009;Boisier et al., 2012;de Noblet-Ducoudré et al., 2012).The simulated changes in ET were found to vary in both magnitude and sign across the various models, and from season to season, despite the fact that none of the models included irrigation or detailed cropland management.About one-third of this inter-model dispersion was explained by differences in the land-cover maps prescribed in each model, and the remaining two-thirds by their distinct land surface parameterizations and resultant model sensitivities to LULCC (Boisier et al., 2012).
Most impacts of historical LULCC on continental water budgets reported in the literature have been addressed through modelling studies.Few studies have estimated the historical large-scale ET changes based on observations.Gordon et al. (2005) calculated global ET change induced by deforestation and irrigation separately.They estimated a moderate decrease in total land ET of 400 km 3 yr −1 , resulting from the opposing effects of deforestation (−3000 km 3 yr −1 ) and irrigation (+2600 km 3 yr −1 ).Sterling et al. (2013) estimated the changes in terrestrial ET between potential and actual land cover based on a large and diverse record of ET measurements, including values from irrigated crops.They calculate an annual mean global ET decrease of 3500 km 3 yr −1 , i.e., a value substantially larger in amplitude than the Gordon et al. (2005) estimate.
The short period covered by ET observations and the limited number of measuring sites, explain the scarcity of studies addressing large-scale changes in ET based on observations.Nevertheless, several global gridded ET products have been recently produced, some of them being outputs of land surface models (LSMs) forced with atmospheric observations or reanalysis.Other products use simpler (semi-empirical) models to diagnose ET from surface and satellite observations of key drivers of ET.A number of them also use available ET measurements for calibration (e.g., Jung et al., 2010;Zhang et al., 2010).Many of these gridded products of ET have been evaluated in the context of the LandFlux-EVAL initiative (Jiménez et al., 2011;Mueller et al., 2011Mueller et al., , 2013)).
The aim of the present study is to provide new data-driven estimates of the past ET changes caused by LULCC (i.e., other factors such as climate and CO 2 changes being constant), and compare them with results from the LUCID climate models.We also explore the seasonal and geographical distribution of the inferred ET changes, as well as quantifying the uncertainties related both to the nature of the observation-based ET data set used and to the land-cover maps adopted.To address these objectives we start with three global products of present-day ET, each of them derived with a different approach.A multivariate regression technique is used to construct empirical models of ET as a function of key environmental drivers and land-cover properties (vegetation partitioning and leaf area index -LAI).These models, along with a set of land-cover maps constructed within the LUCID project for 1870 and 1992, are used to diagnose multiple present-day and preindustrial ET climatologies.

Observation-based data sets
Several gridded data sets are used in this paper (summarized in Table 1).Three products of ET are used in a multivariate analysis to derive, respectively, three different empirical models of ET (described in Sect.2.3).
The "Global Land surface Evaporation: the Amsterdam Methodology" ET (hereafter GLEAM) uses a semi-empirical approach to derive total ET from both the soil-vegetation ET and the evaporation of rainfall intercepted by the canopy (Miralles et al., 2011).A modified Priestley-Taylor ET model is nurtured with meteorological data and a computed soil moisture stress factor.The input data set used to derive gridded ET includes information of land-cover (short versus tall canopy, vegetation optical depth), soil moisture and other meteorological variables.
The Max Planck Institute ET product (MPI) is based on the analysis of ET data from 253 eddy-covariance measurement sites (Jung et al., 2010).In situ observations of ET and of a large set of explanatory variables, including land surface properties (vegetation type and optical properties) and climatic variables, are analysed with a model tree ensemble (MTE) approach.The resultant MTE models are used to create gridded maps of monthly ET from 1982 to 2008, based on global data sets of the associated predictors (surface analyses and remote sensing data).
As with GLEAM, the Numerical Terradynamic Simulation Group ET product (NTSG) is based on a semi-empirical model of ET (Zhang et al., 2010).In this case, a modified Penman-Monteith approach is adopted to explicitly calcu-late the soil evaporation and the canopy ET components.A Priestley-Taylor model is also used for evaporation over water bodies.Global fields of ET from 1983 to 2006 are derived using remotely sensed land-cover data, NDVI and radiation, as well as meteorological data from reanalysis.Local measurements of ET from 34 FLUXNET sites are also used in NTSG to calibrate canopy conductance as a biome-specific parameter.
A short evaluation of the three ET products is presented in Appendix A. Despite the different procedures used to derive these data sets, there is a reasonably large agreement on the climatological mean distribution of ET.However, significant differences from one product to another are found in the amplitude of the interannual ET variability (Appendix Fig. A1).
Gridded data sets of surface incoming solar (S D ) and longwave (L D ) radiation, precipitation (P ) and snow water equivalent (SWE) are used as environmental constraints in our empirical models of ET.The satellite products of the Surface Radiation Budget (SRB) project and of the National Snow and Ice Data Center (NSIDC) are used for radiation fluxes and SWE, respectively.The observation-based gridded product of the Global Precipitation Climatology Centre (GPCC) is adopted for P .
Other variables, such as the net radiation or the nearsurface temperature, although suitable predictors for ET, were deliberately omitted because of their large dependencies on the type of land cover.Including them could thus have produced misleading results when evaluating the changes in ET due to LULCC.In our approach, the surface albedo and the resultant net radiation are implicitly accounted for through the land-cover partitioning.Hence, the radiative www.hydrol-earth-syst-sci.net/18/3571/2014/ Hydrol.Earth Syst.Sci., 18, 3571-3590, 2014 effect of LULCC in the diagnosed ET is included since the albedo-induced change in net radiation will follow the prescribed change in land cover.Land-cover for the purposes of this study has been simplified and grouped within five classes: evergreen trees, deciduous trees, grasses, crops and bare soil.Our choice to use only those main groups is intended to simplify the analysis and give consistent land-cover partitioning in different data sets, notably within the various plant functional types (PFTs) used in LUCID LSMs (Sect.2.2).This simplified land-cover partitioning allows capturing the spatial ET variability induced by differences in plant properties, such as canopy conductance, root depth, surface roughness or albedo.
In addition to the grid areal fraction occupied by the five classes of land cover (F v ), the spatial and seasonal LAI distribution for each land-cover class is also used to characterize the global distribution of vegetation.LAI is not used independently from F v .Instead, we make use of the Beer's law to combine the two variables and define an effective land-cover fraction per vegetation class (F * v ): where L v is the specific LAI of the vegetation group v.The light extinction coefficient k is set to the commonly used value of 0.5.The effective area occupied by the vegetation in a grid cell is then defined as the total fraction derived from the four classes of vegetation, and the effective bare soil area as (2) We use the MODIS land-cover product MCD12Q1 (Friedl et al., 2010) and the reprocessed MODIS LAI of Yuan et al. (2011) to derive monthly mean (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) maps of F v and L v , respectively.Based on the PFT classification of the MODIS land cover, we aggregated the areal fractions of the different PFTs into the main classes used here.The maps of L v were derived from both the LAI and land-cover data from MODIS, following a methodology adapted from Boisier et al. (2013).The application of this method to create biomedependent LAI maps, as well as its evaluation, is described in Appendix B.

LUCID simulations
Global climate simulations performed in the context of the LUCID project (Pitman et al., 2009) are intended for two purposes.We want first to estimate the past LULCC-induced ET changes together with the associated uncertainties arising from historical land-cover reconstructions.Results from LU-CID have indeed demonstrated that current climate models implement the historical crop and pasture data sets in different ways, increasing the spread in the simulated impacts of LULCC, notably on ET (de Noblet-Ducoudré et al., 2012).
To account for this uncertainty, we use the same vegetation Reconstructed ΔET"

estimates
Data-driven ET LC (6)   1870 & 1992" From 6 LSMs (LUCID) maps of 1870 and 1992 prescribed in each LUCID land surface model (LSM) to force the observation-based empirical models of ET derived here (Sect.2.3) and, with them, to perform several reconstructions of past and present-day ET (Fig. 1).
Secondly, we use the data simulated within LUCID to contrast the modelled historical changes in ET with the ones diagnosed here.
The set of LUCID simulations used here were carried out by six atmospheric global circulation models (AGCMs) coupled to LSMs (Table 1).The modelling experiment includes two types of simulations, each with the same prescribed present-day sea-surface temperature/sea-ice coverage (from 1970 to 1999) and atmospheric CO 2 concentration (set to 375 ppm).The types of simulation differ only in the land cover prescribed in the corresponding LSMs, respectively representing the vegetation distribution of 1870 and 1992.The land-cover distribution results from the combination of the "background" map of vegetation of the host model (potential or observed, depending on the LSM) with the historical cropland and pasture distribution of the SAGE (Ramankutty and Foley, 1999) and HYDE (Goldewijk, 2001) data sets.Hence, the historical LULCC that is finally prescribed varies from model to model due to differences in the Hydrol.Earth Syst.Sci., 18, 3571-3590, 2014 www.hydrol-earth-syst-sci.net/18/3571/2014/ background land-cover maps and in the strategies adopted to include the agricultural data.Further details of the LUCID simulations are given by de Noblet-Ducoudré et al. (2012).All the LUCID data were interpolated and analysed in a common rectangular grid of 2.0 • .

Multivariate analysis and evapotranspiration reconstructions
The goal of this study is to calculate LULCC-induced changes in ET.This implies that we develop a method that calculates ET given a specific land-cover map, at the global scale.The methodology we developed to achieve such calculations is presented in this section and summarized in Fig. 1.
We use an implementation of Multivariate Adaptive Regression Splines (MARS) for Python (py-earth) to construct empirical models of ET based on the present-day data described in Sect.2.1 .This technique builds an additive model of piecewise linear (hinge) functions of the predictor basis, capturing non-linear dependencies of a given variable to multiple explanatory variables (Friedman, 1991).The coefficients of the basis are calculated in a two-step iterative computation that uses a least squares method to minimize the error between partial predictions and observations.The first step (forward phase) adds a large number of basis functions, while the second step (pruning phase) selects the best basis sample, leading to a simpler model with a high predictive performance.
Some advantages of MARS, compared to other multivariate methods that could have been used here (e.g., regression trees or neural networks), are the high predictive capabilities, transparency and continuous character of the model constructed (Hastie et al., 2009).
The global ET products presented in Sect.2.1 are used to derive three different MARS-based ET models as a function of a unique set of explanatory variables (Fig. 1).In the computation of these models we account for both the spatial and temporal variability of the input variables.The effective fraction of the five land-cover classes used (F * v ) are not included as independent predictors in the ET models, but are used to weigh the environmental predictors.Further details on the MARS analysis and an evaluation of the predictive skill of the resulted ET models are described in Appendix C.
Preindustrial (PI) and present-day (PD) reconstructions of ET are then computed by forcing each MARS-based ET model with the land-cover maps of 1870 and 1992 used by six LSMs in the LUCID simulations (Sect.2.2).Hence, we derive 18 different pairs of reconstructed ET that only differ in the choice of the F v values used in Eq. ( 1) and the ET information deduced from each product.
It should be noted that all reconstructions are carried out with the same environmental data used in the MARS computations (Sect.2.1) and the MODIS-based LAI of each land-cover class (L v ).Therefore, the diagnosed ET changes between PI and PD represent instantaneous responses to LULCC, and do not include any biogeophysical feedback involving a perturbation in downward radiation, precipitation or snow cover, nor in a biome-dependent LAI.(Ramankutty and Foley, 1999).In most LSMs, such change mainly occurs at the expense of forest.Given that pastures are considered as natural grasses in most LU-CID LSMs, the prescribed change in grass area from 1870 to 1992 results from the balance between the positive change when natural biomes are converted to managed pastureland, and the negative change when natural grasses are converted to croplands.The rules defining how forest and natural grassland are decreased when crops and pasture expand are thus crucial in drawing up the final figure for deforestation, as discussed by de Noblet-Ducoudré et al. (2012).For the LSMs assessed here, the global deforested area between 1870 and 1992 varies strongly, ranging from ∼ 4 to 10 million km 2 for JS-BACH and TESSEL, respectively.These two examples illustrate the discrepancy in the LULCC reconstructions that result from different protocols.When the agricultural units (crops and pastures in this case) are allocated over natural grasslands as priority, the resulting change in forest area is comparatively low (as in the JSBACH case).The opposite case (favoured deforestation) occurred in the case of TES-SEL.
Foliage density is a key parameter in both the simulated and diagnosed ET.As described in Sect.2.1, the present-day LAI as well as the changes in LAI from 1870 to 1992 are implicitly accounted for in the diagnosed ET through Eq. ( 1).Yet, the change in LAI as diagnosed with our MODIS-based reconstructions is in itself an interesting result.Such an estimate is also important if we are to correctly interpret the diagnosed LULCC-induced changes in ET.Hence, based on the LUCID land-cover maps and the seasonally and geographically varying LAI values we derived from MODIS data for each of the four classes of vegetation (L v , see Appendix B), we computed the monthly distribution of LAI for 1870 and 1992 as (3)   Since deforestation is the dominating perturbation between 1870 and 1992, and given that forest LAI is usually larger than that of short vegetation (see Appendix Fig. B1), most areas of the globe show decreases in LAI, and the spatial pattern of such changes is similar to that of deforestation.Increases in LAI between 1870 and 1992 are also diagnosed, notably in summer and over regions where LULCC is dominated by transitions from natural grasses to crops with a higher LAI (the occurrence and amplitude of such transitions depends on the land-cover data set chosen; not shown).
Looking at Fig. 2d in detail, the distribution of L highlights some regions that could be particularly sensitive to LULCC in terms of ET, such as the southern part of North America or southeastern Amazonia.These regions show comparatively large decreases in LAI, despite moderate deforestation between 1870 and 1992.As described next, the impact of LULCC on the mean annual LAI in temperate and boreal regions results mainly from the changes in summer, when canopies are fully developed.
Figure 3 shows the monthly mean LAI reconstructed for 1992 and L averaged over four regions we are particularly interested in because of their strong historical LULCC.These regions, as indicated in Fig. 2d, correspond to the land areas within defined domains in Eurasia (hereafter EA), in North America (NA), in South America (SA) and in South-east Asia (SEA).The reconstructed mean LAI depicted in Fig. 3a (solid lines) is plotted along with the ones simulated (prescribed in some cases; de Noblet-Ducoudré et al., 2012) by the LUCID LSMs (dashed lines).In both cases, the mean obtained from the six LSMs ±1.0 inter-model mean absolute deviation (MD) are shown.The MODIS-based reconstructed LAI shows a clear seasonal cycle in the northern temperate regions (EA and NA), with a winter minimum of near null LAI and a maximum in summer of around 2.0 m 2 m −2 .In these two regions there are large differences between the present-day reconstructed and simulated LAI.The latter overestimates the former by nearly 1.0 m 2 m −2 throughout the year.The mean LAI in SA and in SEA are larger than those observed in EA and NA both in the models and reconstructions, and without clear seasonal cycles, indicating a large contribution of tropical and subtropical evergreen forest in these regions.
A large spread between the LUCID models LAI compared to the reconstructions is also noticeable in Fig. 3a.Given that the reconstructed LAI accounts for the various LUCID LSMs' maps of land cover, the spread in the simulated LAI highlights intrinsic model dependencies (parameterizations).Although remarkable, this spread is not surprising, given the very different treatments of the vegetation phenology included in the LUCID LSMs assessed here.For instance, three of them compute LAI (JSBACH, LPJmL and ORCHIDEE) while the other three models prescribe it.Further, those models prescribing LAI have used different data sets as reference,  These results echo recent model inter-comparison studies reporting large uncertainties and a systematic overestimation in modelled LAI when compared to satellite data (e.g., Anav et al., 2013).Besides the known uncertainties in the observation-based LAI data sets (e.g., Garrigues et al., 2008;Fang et al., 2013), Anav et al. (2013) suggest that the model LAI biases could in part be attributed to missing land-surface parameterizations accounting for nutrient (Nitrogen) limitation or ozono effects on gross primary production; both of these processes would increase the simulated plant carbon allocation and LAI.
The reconstructions show clear decreases (∼ 10 %) in LAI from 1870 to 1992 in the four regions assessed as having intense LULCC (solid lines in Fig. 3b).In EA and NA, the amplitude of L is maximized in the corresponding summer, while in SA and SEA there is a substantial year-round decrease.The effect of using different LULCC data is particularly important in the reconstructed L over NA and EA.The MD between the individual results (around 0.15 m 2 m −2 in summer; grey shaded area in Fig. 3b) is as large as the mean L.
Also noticeable is how different the reconstructed L are compared to the simulated ones (dashed lines in Fig. 3b).This is particularly clear in EA and NA, where the simulated mean L shows an opposite seasonal pattern with respect to the one reconstructed.In the other two regions the models underestimate L in all seasons when compared to the reconstructions.
The contrasting diagnosed versus modelled LAI responses to LULCC are in part explained by an overestimation of the model cropland LAI during the growing season compared to the values derived from MODIS -a feature that is particularly marked in those LSMs that simulate LAI (de Noblet-Ducoudré et al., 2012).However, it should also be noted that the simulated LAI responds to changes in both land cover and climate (accounting for interactions/feedbacks between land surfaces and the atmosphere), while the latter is not accounted for in the MODIS-based reconstructions derived in this study.

Diagnosed changes in evapotranspiration
Mean annual changes in ET ( ET) between the preindustrial period (PI) and present-day (PD) were calculated based on each of the three ET products, and using the six pairs of land-cover maps from the LUCID LSMs (Fig. 1).Those 18 reconstructed ET climatologies were averaged to produce the mean annual ET displayed in Fig. 4.This multi-product mean highlights the main patterns of LULCC-induced ET.A decrease in annual ET is observed in most areas with landcover perturbations; this is true for South America, Africa, India and Oceania.However, in the northern temperate regions the signal is not as strong, except in the southernmost regions of change in North America.The annual mean ET is moderate in most cases, and only the regions mentioned above show relative anomalies above 5 %.The spread obtained from the various estimates is quite large and of an order of magnitude similar to the diagnosed mean annual ET  (see contour line in Fig. 4).As described next, the nature of the LULCC data adopted represents the main source of uncertainty in our estimates.Table 2 summarizes the global changes in annual ET in terms of volume of terrestrial water vapour.We diagnosed from the reconstructions an average decrease in global ET of around 1250 km 3 yr −1 between PI and PD.Relative to this value, we obtain a large uncertainty between individual reconstructions (MD of 850 km 3 yr −1 ), arising both from the ET product used (MD of 470 km 3 yr −1 ) and from the LULCC data set adopted (MD of 780 km 3 yr −1 ).
For comparison, Table 2 also indicates the global ET simulated by the LUCID AGCM/LSMs, as well as the ones diagnosed in previous studies based on observations.LUCID simulations show a model-mean ET of −760 km 3 yr −1 , with a large inter-model dispersion (MD of 720 km 3 yr −1 ).Both the diagnosed and simulated global ET are substantially weaker in amplitude than the recent estimate of Sterling et al. (2013).
Although irrigation is not explicitly considered in our diagnosed ET, it should be partially accounted for in an implicit way given that two ET products (MPI and NTSG) have used in situ ET observations for calibration, including irrigated crops.Considering irrigation, our diagnosed ET fall within the estimates of Gordon et al. (2005, −400 km 3 yr −1 ) and Sterling et al. (2013, −3500 km 3 yr −1 ).However, in view of the large uncertainties associated with the LULCC data sets, comparisons between independent estimates might be misleading.Further, the two references mentioned above have computed past ET climatologies with potential (preagricultural) vegetation maps and, in consequence, their prescribed land-cover perturbations between the past and the present are likely to be larger in amplitude than the ones used in this study.
As can be seen in Fig. 4 Figure 5 illustrates the monthly mean ET diagnosed with the three ET products on average over the four regions of study.Clear differences are observed between the various estimates.There are, however, some patterns that distinguish the results from MPI and NTSG on the one hand, and that from GLEAM on the other.The reconstructed ET based on MPI and NSTG, although with biases between them, show similar seasonal patterns in EA, NA and SA.In EA, ET is characterized by positive values during the early boreal summer and a minimum (a decrease in ET in the case of MPI) during the autumn.In NA, the seasonal pattern of ET is even clearer: the results from these two products show a decrease and an increase in ET in spring and in summer, respectively.In SA, the MPI and NSTG-based ET also show a clear seasonal pattern, characterized by negative values during most of the year, particularly large ET decreases in the austral winter, and increases in late summer.
The ET reconstructions based on GLEAM show yearround decreases of ET between PI and PD in all four regions assessed.Although ET derived with this product is clearly different from those based on MPI and NTSG, the seasonality shows some similarities with the other products in EA and NA.By contrast, the result from GLEAM in SA shows a year-round decrease in ET that clearly differentiates it from the other two products.
For comparison, the multi-product mean ET is also plotted in Fig. 5, as well as mean ET resulting from the LU-CID climate model simulations.The uncertainties are large and of the same order of magnitude in both cases.Nevertheless, some consistent signals of ET can be determined, such as the year-round decreases in SA and SEA.The distinct seasonal pattern of ET in EA is suggested in both the reconstructed and the simulated results.By contrast, the simulated mean ET does not show a clear pattern in NA, as compared to the one shown by the reconstructions.
The regional means ET depicted in Fig. 5 suggest some robust seasonal features, but they still mask some contrasting ET responses to LULCC within the regions assessed, notably across the areas affected by strong LULCC in northern temperate latitudes.This can be observed in Fig. 6, which presents the spatial distribution of the seasonal (2-monthly) mean ET in North America and Eurasia.The ET differences between PI and PD derived from the three assessed ET products show positive and negative anomalies of large amplitudes during the late spring and summer.The North American Great Plains are particularly affected by large ET values in all the three reconstructions, with anomalies exceeding 10 % in many areas (the relative changes are not shown).Clear differences can be observed between the various estimates.Those based on GLEAM (Fig. 6a) and NTSG (Fig. 6c) are at both extremes, respectively showing the strongest negative and positive ET values.Meanwhile, the reconstruction based on MPI (Fig. 6b) shows clear simi-larities in the seasonal and spatial patterns of ET with that based on NTSG, in accordance with the regional mean ET shown in Fig. 5.Both cases present a clear late summer maximum in ET over the North American area of large cropland expansion (Fig. 2).
The distribution of ET shown in Fig. 6 reveals spatially coherent signals compared to specific land-use transitions between 1870 and 1992.As mentioned above, increases in ET are obtained in the MPI and NTSG reconstructions in regions where crops were partially allocated in place of grass, such as in northern North America (Fig. 2).On the other hand, regions with an increase in grass and decrease in forest, such as in southern North America or in the mid-Eurasian area north of the Caspian Sea (areas of large expansion in pasture), show systematic decreases in ET in all the three reconstructions.To better understand the different ET sensitivities to LULCC estimated from the three data sets assessed, in the section that follows we try to quantify the ET response associated with three different types of land-cover transitions.

Sensitivity of evapotranspiration to specific land-cover transitions
Although the ET products assessed show a general agreement in the global distribution of ET (Appendix Fig. A1), specific biome-dependent ET should explain the differences between the three diagnosed ET responses to LULCC.These differences arise when the long-term mean ET rates of the various products are compared over dominant types of land cover.Figure 7 illustrates this comparison for three major groups of vegetation: crops, grasses and forest.Regions with dominant land cover are defined as the grid cells (at 1.0 • of resolution) showing at least 75 % of their area covered by the corresponding class.Averaged over the northern extratropical regions (above 20 • N), the monthly mean ET of the three products shows similar rates and seasonal patterns over grasslands and forest (Fig. 7).However, a larger contrast is observed over croplands, where the ET from GLEAM clearly underestimates those from the other two products.This difference is consistent with the diagnosed ET response to LULCC in areas of large cropland expansion, larger (with positive anomalies in some cases) in the cases of MPI and NTSG with respect to that of GLEAM (Fig. 6).However, this analysis does not allow direct comparisons of ET across the various land-cover classes, because they were obtained from different regions and therefore from different climate regimes.
In order to quantify the response of ET to specific changes in land cover under equivalent environmental conditions, we have examined at each grid cell the transition between 1870 and 1992 from a given land-cover type (A) to another (B) using the following rule: where F v is the change in the areal fraction of a generic land-cover type v. Hence, if A is totally converted to B, then δ A→B equals 1.0.If A is partially converted to B, and no other transition occurs simultaneously in a given grid cell, the change in ET as a response to a total transition from A to B can be estimated as the ratio between the actual ET and δ A→B .Figure 8 illustrates how the diagnosed ET in the northern extratropics relates to the three major land-use transitions we are particularly interested in: forest-to-grass, forest-to-crops and grass-to-crops.Figure 8a shows, as an example, the local (grid cells) ET deduced from the MPI product in July, plotted as function of each type of transition.In order to avoid misleading results from mixing simultaneous transitions, the only grid cells retained are those for which the selected transition is, as a minimum, four times larger in amplitude than the fractional area change of all other (not involved) land-cover units (this factor was defined by inspection with the criteria of retaining a significant number of grid cells for the analysis).This example illustrates the dominant ET response to the specific transition selected.That is, a decrease in ET when forests are replaced with grasses, an increase in ET when cropland is allocated at the expense of grasses, and a no clear signal to the forest-tocrops transition.The normalized mean ET responses to each transition (sensitivity) may then be quantified as the slope of the linear fit between ET and δ A→B .This analysis is illustrated as dashed lines for the particular case shown in Fig. 8a, and generalized for all seasons and product-based ET reconstructions in Fig. 8b.
For the northern extratropical regions, a similar behaviour is observed for the three ET products for a transition from forest to grassland, which is characterized by decreases of ET throughout the year but maximized in the boreal summer, when the mean anomalies reach ∼ −35 mm month −1 .This preferred ET response to deforestation is consistent with what is expected from reduced LAI and the underlying decrease in transpiration and in evaporation from intercepted rainfall.By contrast, when a land-cover transition involves changes in the area under crops, the ET sensitivity shown by the three products disagrees both in sign and magnitude.The results based on NTSG and GLEAM respectively show a summertime increase and decrease in ET for a transition from forest to crops, while MPI produces a weak ET response.When grasses are converted to crops, NTSG and MPI show clear increases in ET, while GLEAM-based anomalies are weak compared with the other products.These patterns highlight then different estimates of crops' ET rates across the various products assessed, in line with the results shown in Fig. 7.

Hydrol
The results shown in Fig. 8b help to interpret seasonal patterns of ET derived with the reconstructions.The distinct ET response to LULCC derived with the MPI and the NTSG products, which is characterized by summer ET increases in areas of large expansion of cropland, agrees with their ET sensitivities to transitions from forest or grassland to cropland.Hence, the results based on NTSG show large positive ET anomalies, particularly when cropland replaces grassland (Fig. 6).The reconstruction based on MPI shows a similar ET response to these types of land-cover transitions but of lower amplitude.In turn, the ET sensitivity based on GLEAM reveals weak to negative ET anomalies for the three major land-cover transitions from 1870 to 1992, explaining the dominant ET decreases between the PI and PD diagnosed with this product.In summary, the analysis described in this section suggests that the disagreement between the changes in ET deduced with the different products is principally a result of their specific cropland ET estimates.Moreover, our study provides seasonal and spatial details of large-scale ET changes that have not been discussed in any earlier work based on observations.Our results have demonstrated that while most parts of the globe show annual mean ET decreases, extensive areas in the Northern Hemisphere extratropics have experienced ET increases, specifically during the growing season in regions of large historical cropland expansion.In those regions the impact of past LULCC on annual ET masks strong seasonally varying changes in ET and points to the necessity of having access to observation-based reconstructions at the seasonal timescale.
Previous results from the LUCID intercomparison project have revealed very large uncertainties in modelled ET response to LULCC between the preindustrial period and the present, and an important fraction of these uncertainties was attributed to the reconstructed historical scenarios of LULCC (Boisier et al., 2012).To account for this uncertainty, the past changes in ET were estimated here using six different historical scenarios of LULCC previously used in LUCID.
We have diagnosed a global land annual ET change of around −1250 km 3 yr −1 between the preindustrial period and present-day.This value is larger in amplitude than the simulated mean response to LULCC obtained from the LU-CID climate models, and is placed between the data-driven  estimates of Gordon et al. (2005)  In addition to the land-cover reconstructions, another source of uncertainty in the diagnosed ET change arises from intrinsic sensitivities to LULCC deduced from the various ET data sets adopted here.Our results show that these sensitivities are mainly related to the products' present-day crop ET estimates; we therefore highlight the necessity of revisiting how data sets treat crops.
The large differences in LAI shown by the LUCID models, and the associated uncertainties in simulated ET, reiterate previous findings and point out the need for in-depth evaluations of the vegetation phenology simulated in LSMs (e.g., Richardson et al., 2012;Anav et al., 2013).
Including irrigation is crucial to proper assessments of past land-use-induced ET changes (Gordon et al., 2005).Although not explicitly accounted for in this study, irrigation should be partially included in our ET reconstructions based on the MPI and NTSG products, given their FLUXNET data calibration.Regarding this aspect, it is noteworthy that the diagnosed ET based on these two products, despite the very different procedures used to derive them, show similar sensitivities to LULCC and, therefore, consistent spatial and seasonal patterns of ET change between the preindustrial period and the present.
Given that the preindustrial ET reconstructions were calculated using present-day data of the environmental drivers (precipitation, radiation and snow cover) and, implicitly, with current atmospheric CO 2 concentrations, the estimated ET changes do not consider any feedback that involves these drivers.This assumption could be particularly important for precipitation.Previous results from LUCID have shown that precipitation responds to LULCC synchronously with the changes in ET, hence amplifying the impacts on ET when atmospheric feedbacks are accounted for (Boisier et al., 2012).
The set of diagnosed ET presented here were derived with empirical multivariate models of ET.It is important to recall that the ET products used to derive these models were themselves obtained with empirical or semi-empirical methods.Although a number of site-level ET measurements were used to calibrate two products, these remain very limited in a global-scale context.Considering this, our results should be carefully interpreted since they involve uncertainties inherent both in the multivariate analysis carried out here and in the nature of the ET product used.Accounting for the specific product-based ET sensitivities to LULCC is therefore crucial.This is why we apply the same methodology to three ET products originally derived in quite different ways.
The increasing numbers of ground-based observations and satellite data, combined with statistical tools, allow accurate estimates of the current large-scale ET to be derived (e.g., Jung et al., 2010).Here, we have demonstrated that similar methods are also suitable for constraining uncertainties in the historical changes in ET, bringing a new class of estimates independent of global climate model simulations.Similar methods could also be applied to evaluate the historical impact on other key variables of the climate system (e.g., Boisier et al., 2013), driven by LULCC or by another climate forcing, as well as to perform future projections.

Appendix C: Multivariate ET models
As described in Sect.2.3, three empirical ET models were constructed using a multivariate regression tool based on Multivariate Adaptive Regression Splines (MARS).Each of them resulted from the analysis of each ET product assessed here (Table 1).To construct those models, we first considered the mean ET in a grid cell (E g ) as the linear combination of the components associated with each land-cover class (E v ).That is, where x represents an array of the environmental predictors, including the monthly mean, the long-term monthly mean and the long-term annual mean of P , L D , S D and SWE.If E v is an additive model of the basis defined by the elements of x (as the hinge functions used in MARS), we can rewrite E g as with α 0 and α v,i K , the parameters of the MARS model computed for each of the 60 basis, resulting from the 5 landcover classes [N (v)] by 12 environmental predictors [N (i)].H K represents the hinge functions obtained for a given basis.The knots position and number [N (K)] are automatically selected in the MARS routine.
The MARS analysis is performed with the monthly data of ET and the explanatory variables from 1984 to 2006, the overlapping period between all data, except for F * v , for which we have used the mean (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) monthly values derived from MODIS as present-day data (Table 1).The complete data set was previously regridded (averaged) onto a 1.0 • rectangular projection, since all data are available at equal or higher resolutions (Table 1 The full data set comprises more than 3 × 10 6 observations (number of pixel-months).In order to increase the computational efficiency, we use a random subset of the full record (about 9 %) as input during the training process.Preliminary tests showed that the predictive performance of the constructed models is not improved by using larger samples in the training data.
In order to evaluate the predictive skill of the MARS-based ET models, we reconstructed ET with each model and the complete predictor data set, and we compared them to the actual ET.Appendix Fig. C1 shows the scatter diagrams between the observed ET (product) and the reconstructed ones.Considering the full record (i.e., collecting together all the pixel-month values from 1984 to 2006), the reconstructed ET explains about 90, 91 and 95 % of the observed ET variance in the case of GLEAM, NTSG and MPI, respectively.The typical errors (MAE) are, in the same order, 8.2, 6.9 and 5.3 mm month −1 .These errors are reduced by a factor ∼ 2, when the climatological monthly values are considered, indicating that an important fraction of the errors occurs at the interannual timescale.This is somehow expected since the input data set record is much larger in space (about 11 000 pixels) than in time (276 months).The skill of MARS models to predict the spatial and seasonal distribution of ET is reasonably good despite the relative simplicity of the set of explanatory variables used and their independence with respect to the ET data sets.The environmental data adopted here are in most cases different to those originally used to derive the ET products (Table 1).Yet, the typical error of the reconstructed ET climatologies, ranging from ∼2.9 mm month −1 (MPI) to ∼4.1 mm month −1 (GLEAM), remains lower but of the same order than the diagnosed changes in ET in regions of important land-cover perturbations (see Fig. 6).
Besides the predictive skills of the MARS ET models derived here, an important fraction of the resulting errors in the ET reconstructions should be attributable to the omission of key drivers of ET originally accounted for in a given product, such as the soil moisture in the case of GLEAM.In contrast, we consider that the contribution of the land-cover partitioning to the spatial ET variability of the products is fairly well captured by the reconstructions and, therefore, the MARS models are able to estimate changes in ET driven by LULCC.

Figure 1 .
Figure 1.Methodology followed to estimate the LULCC-induced evapotranspiration (ET) changes between the preindustrial period and the present.The data set used (boxes) includes gridded maps of present-day ET, land cover (LC), leaf area index (LAI), radiation (S D , L D ), precipitation (P ) and snow cover (SWE) (see Table 1).Solid arrows indicate the input data set used to compute the MARS-based ET models ( E). Dotted arrows indicate the data used as forcing in the computation of the ET reconstructions.

3
Figure2a-cshow the change between 1870 and 1992 in the fraction of crops, grasses and forest as deduced from the LUCID land-cover maps (the mean changes deduced from the various LSM PFT maps are shown).The expansion of cropland, notably in North America and in western Eurasia, is a main feature of LULCC between 1870 and 1992.The change in crop area prescribed in LUCID LSMs reflects this pattern, in line with the historical land-use change deduced from SAGE(Ramankutty and Foley, 1999).In most LSMs, such change mainly occurs at the expense of forest.Given that pastures are considered as natural grasses in most LU-CID LSMs, the prescribed change in grass area from 1870 to 1992 results from the balance between the positive change when natural biomes are converted to managed pastureland, and the negative change when natural grasses are converted to croplands.The rules defining how forest and natural grassland are decreased when crops and pasture expand are thus crucial in drawing up the final figure for deforestation, as discussed by deNoblet-Ducoudré et al. (2012).For the LSMs assessed here, the global deforested area between 1870 and 1992 varies strongly, ranging from ∼ 4 to 10 million km 2 for JS-BACH and TESSEL, respectively.These two examples illustrate the discrepancy in the LULCC reconstructions that result from different protocols.When the agricultural units (crops and pastures in this case) are allocated over natural grasslands as priority, the resulting change in forest area is comparatively low (as in the JSBACH case).The opposite case (favoured deforestation) occurred in the case of TES-SEL.Foliage density is a key parameter in both the simulated and diagnosed ET.As described in Sect.2.1, the present-day LAI as well as the changes in LAI from 1870 to 1992 are implicitly accounted for in the diagnosed ET through Eq. (1).Yet, the change in LAI as diagnosed with our MODIS-based reconstructions is in itself an interesting result.Such an estimate is also important if we are to correctly interpret the diagnosed LULCC-induced changes in ET.Hence, based on the LUCID land-cover maps and the seasonally and geographically varying LAI values we derived from MODIS data for each of the four classes of vegetation (L v , see Appendix B), we computed the monthly distribution of LAI for 1870 and 1992 as

Figure 2 .
Figure 2. Differences between 1992 and 1870 in the fractional area (%, absolute) covered by (a) crops, (b) grass and (c) forest.Maps of land-cover change correspond to those prescribed in LUCID LSMs (model mean).(d) MODIS-based reconstructed change in annual mean leaf area index from 1870 to 1992.Contour lines indicate the four regions used later for specific analysis.

Figure
Figure 2d shows the change in annual mean LAI [ L = L(1992) − L(1870)] computed with the six LUCID LULCC data sets (the average from the six reconstructions is shown).Since deforestation is the dominating perturbation between 1870 and 1992, and given that forest LAI is usually larger than that of short vegetation (see Appendix Fig.B1), most areas of the globe show decreases in LAI, and the spatial pattern of such changes is similar to that of deforestation.Increases in LAI between 1870 and 1992 are also diagnosed, notably in summer and over regions where LULCC is dominated by transitions from natural grasses to crops with a higher LAI (the occurrence and amplitude of such transitions depends on the land-cover data set chosen; not shown).Looking at Fig.2din detail, the distribution of L highlights some regions that could be particularly sensitive to LULCC in terms of ET, such as the southern part of North America or southeastern Amazonia.These regions show comparatively large decreases in LAI, despite moderate deforestation between 1870 and 1992.As described next, the impact of LULCC on the mean annual LAI in temperate and boreal regions results mainly from the changes in summer, when canopies are fully developed.Figure3shows the monthly mean LAI reconstructed for 1992 and L averaged over four regions we are particularly interested in because of their strong historical LULCC.These regions, as indicated in Fig.2d, correspond to the land areas within defined domains in Eurasia (hereafter EA), in North America (NA), in South America (SA) and in South-

Figure 3 .
Figure 3. (a) Monthly mean LAI in 1992 and (b) LAI change from 1870 to 1992 averaged over four regions of study (indicated Fig. 2d).Reconstructed (MODIS-based) and simulated (LUCID) LAI are shown as solid-black and dashed-green lines, respectively.Mean LAI is plotted along with a range of ±1.0 mean absolute deviation, resulting from individual reconstructions (based on the different LULCC data; grey shading), and from the individual model results (green shading).

Figure 4 .a
Figure 4. Land-use-induced change (PD-PI) in annual evapotranspiration ( ET) estimated from multiple products of ET and LULCC maps (the average ET is shown).Contour line encloses the areas where the mean absolute deviation of ET, computed within the various estimates, exceeds 20 mm yr −1 .

Figure 5 .
Figure 5. Monthly mean evapotranspiration (ET) change (PD-PI) averaged over the four regions of study: (a) Eurasia, (b) North America, (c) South America and (d) Southeast Asia.Product-based estimates and the simulated changes in ET are illustrated.Shaded areas indicate the mean ±1.0 mean absolute deviation from the individual estimates or model data.

6
Discussion and conclusions This study presents novel observation-based estimates of LULCC-induced changes in evapotranspiration (ET) between the preindustrial period and present-day, together with associated error bars based on (a) uncertainties in the historical reconstruction of global land-cover distribution, and (b) uncertainties on the data-driven global ET products.

Figure 7 .Figure 8 .
Figure 7. Monthly mean ET associated with dominant classes of vegetation in the Northern Hemisphere extratropics (> 20 • N).Results from GLEAM, MPI and NTSG illustrated as solid, dashed and dotted lines, respectively.Regions with dominant vegetation (areal fraction > 75 % in 1.0 • grid) are shown in the bottom panel.
Figure B1.(a)Reconstructed global distribution of leaf area index (LAI) in July based on MODIS data.(b) Monthly mean LAI averaged over the four regions of study: Eurasia (EA), North America (NA), South America (SA) and southeast Asia (SEA).Reconstructed and observed values indicated as solid line and red dots, respectively.(c) Regional and monthly mean LAI per land-cover class: crops (black line), grasses (red), evergreen trees (green) and deciduous trees (blue).

Figure C1 .
Figure C1.Reconstructed vs. observed monthly evapotranspiration (ET) based on the (a) GLEAM, (b) MPI and (c) NTSG ET products.The full data set (pixel-months) and the long-term monthly mean values are illustrated as grey and black dots, respectively.The mean absolute error (MAE) and coefficient of determination (R 2 ) between the reconstructions and the observations are indicated.
Discussion and conclusions are presented in Sect.6.

Table 1 .
Data set used in this study.

Hydrol. Earth Syst. Sci., 18, 3571-3590, 2014 www
, areas with strong LULCC such as in North America or Eurasia do not show clear annual ET responses to LULCC.The weak annual ET over these regions results from contrasting ET from season to season, and between the various ET product-based estimates. .hydrol-earth-syst-sci.net/18/3571/2014/
and Sterling et al. (2013).However, our reconstructed past changes in ET show a very large dependency on the land-use maps adopted, as previous results from LUCID have also shown.Hence, no straight-forward comparisons can be made between independent estimates that prescribe different LULCC.Constraining the current protocols used to reconstruct maps of land cover and deriving realistic historical scenarios of land-use change Hydrol.Earth Syst.Sci., 18, 3571-3590, 2014 www.hydrol-earth-syst-sci.net/18/3571/2014/