Evapotranspiration simulations in ISIMIP2a—Evaluation of spatio-temporal characteristics with a comprehensive ensemble of independent datasets

We analyse ET from historical simulations of the ISIMIP2a project using all models from the global water, biomes and agriculture sectors that simulate ET using the model's default ET scheme. The simulations used in this analysis are based on three distinct meteorological forcing datasets: the Global Soil Wetness Project Phase 3 (GSWP3, http://hydro.iis.u-tokyo.ac.jp/GSWP3/), the Princeton Global Meteorological Forcing Dataset version 2 (PGMFD v.2, Sheffield et al 2006) and the Water and Global Change (WATCH, Weedon et al 2011) Forcing Data methodology applied to ERA-Interim data (WFDEI, Weedon et al 2014). We use naturalized runs (i.e. without human impact) only, except from a few simulations from biomes models (see table 1). Due to the large number of models involved, we cannot provide a detailed description of each model here, and rather suggest the interested reader to either consider the publications listed in table 1 or to consult the summary information listing the most relevant characteristics of each model (freely accessible via the ISIMIP website, www.isimip.org/impactmodels/).

From the pool of crop model simulations of the ISIMIP2a agriculture sector, we only select harmonized simulations (see Elliott et al 2015 for details). All crop model simulations use fixed fertilizer application rates except for LPJmL and LPJ-GUESS, as fertilizer input is irrelevant for this version of these models. As previously mentioned in table 1, ISIMIP2a crop model simulations are provided as pure crop runs (i.e. assuming a specific crop is growing all over the world), implying that the ET output from these models corresponds to the amount of water evaporated or transpired from a specific crop type under a given irrigation scenario. Among all crop model simulations available within the ISIMIP2a framework, we focus here on ET from non-irrigated maize crops (denoted as ET_maize), as this crop type is among the most dominant crop types in the regions of interest (see section 3.2 and figure A1) and its growing season matches the summertime growing season on the Northern Hemisphere (see figure A2), making the comparison of ET_maize with ET_tot more straightforward. Note that additional preprocessing steps are required prior to performing these comparisons (listed in appendix B).

To additionally include simulations used in the scope of other model inter-comparison projects than ISIMIP, we also include estimates from simulations used in the E2O Water Resources Reanalysis v.1 (Schellekens et al 2017). Please refer to the latter publication for a short description of each of the contributing (land surface, global hydrological and simple water balance) models. Note that SWBM is also part of E2O but has been excluded from this ensemble, as its version is identical to the one of SWBM provided in ISIMIP2a. PCR-GLOBWB and Water GAP are also part of both inter-comparison projects, but kept in both ensembles due to differences in the model versions.

Besides the aforementioned simulations, we also analyse ET output from two major land reanalysis products: ERA-Interim/Land and MERRA-2. ERA-Interim/Land has been selected, as it is commonly used as a reference for quantifying land surface conditions (89 citations as of September 14, 2017, according to Web of Science). MERRA-2 is the most recent reanalysis advancement from NASA that uses refined precipitation corrections Reichle et al (2016). Although ET from MERRA-2 is known to have an anomalously high share of bare soil evaporation Schwingshackl et al (2017), the global average of total ET is well within the range of ET from other reanalyses Bosilovich et al (2016), making it a sufficiently good candidate for our analysis.

We use an ensemble of recent diagnostic datasets suitable for identifying potential biases in model-based estimates of ET (see table 1). This ensemble consists of some well established ET estimates from the recent past (a subset of those diagnostic datasets has also been used to generate the ensemble of diagnostic ET in LFE, see section 2.3), but also includes more recent datasets such as GLEAM v3.1 Martens et al (2016). Please refer to the references listed in table 1 for further details on the individual datasets. Please note that the version of Fluxnet-MTE employed here uses a modified temperature and precipitation forcing compared to the previous version presented in Jung et al (2009).

LandFlux-EVAL (LFE) is an ensemble based land ET product that itself is based on individual diagnostic, model and reanalysis products available in the early 2010s Mueller et al (2013), which we use here as a reference dataset (96 citations as of September 14, 2017, according to Web of Science). Although we cannot argue that LFE is free of biases, we can make the conservative assumption that its ensemble statistics (in particular the provided quantile statistics) are a suitable estimate for quantifying the probable range of ET over land. Please note that some of the diagnostic and model datasets used in the LFE ensemble are also included as individual datasets in this analysis, and hence there is some degree of dependency between the diagnostic (model based) estimates and both LFE DIAGNOSTIC (LFE LSM) and LFE ALL.

In order to allow an inter-comparison of all datasets (mostly available at 0.5° resolution), we have bi-linearly interpolated all input data to a 1° × 1° regular latitude-longitude grid (which corresponds to the spatial resolution of LFE). As LFE can be interpreted as an ET reference product (see section 2.3), we favoured to proceed with this lower resolution.

Besides global land (excluding Antarctica), we focus on the following regions (see table 2): (1) the Great Plains region (spatial delineation according to the USGS physiographic divisions of the conterminous US, https://water.usgs.gov/lookup/getspatial?physio), (2) Central Europe (37°N–53°N and 1°W–20°E, representing the European land area with the highest JJA temperature anomalies in summer 2003 according to Schär et al (2004)) and (3) western Russia (50°N–60°N and 35°E–55°E, according to Dole et al (2011)). These regions have been selected, as they contain large areas of agricultural land that has been affected by severe heat waves and droughts during the time period covered by most of the assessed datasets. The spatial extent of the study areas is also indicated in figures A2 and A1. For each of the regions, area averages of ET were derived by weighting each grid cell by its land area (based on ISLSCP II Global Population of the World, Balk et al 2010).

We apply two distinct approaches to determine the overall level of similarity between the individual datasets focusing on (a) a combination of spatial and temporal variability (replicating the method described and applied in Sanderson and Knutti (2012) and Sanderson et al (2015) for the univariate case) and (b) spatial variability alone (following Mueller et al 2011). For both approaches, we first create a subset only containing data from the period of temporal overlap (1989–2005) using all datasets except from GETA 2.0, MODIS Global ET and WECANN (which are not of sufficient length). If any of the grid cells (from any dataset and time step) contain missing data, values corresponding to those grid cells are removed from all datasets for all time steps. From the m = 75 (m = 89 for ET_maize) datasets remaining, we then proceed as follows:

In approach (a), we first reform the elements of the four-dimensional input field (latitude, longitude, time and dataset) to a two-dimensional matrix X of size m by n (where n corresponds to the number of all non-missing observations) and transform this into anomalies (based on all time steps remaining) ΔX. To only preserve the dominant modes of ensemble variability, a singular value decomposition (SVD) is performed on ΔX and truncated to t = 9 modes (which is well within the range of suitable truncation values found in Sanderson et al (2015)). The corresponding loadings matrix U (size m by t) can now be used to derive the difference matrix D_svd by calculating pairwise Euclidean distances as

$$\begin{equation} {\delta _{ij}} = {\left\{ \mathop \sum \limits^t_{l = 1} {\left[ {{\bf{U}}\left( {i,l} \right) - {\bf{U}}\left( {j,l} \right)} \right]^2}\right\} ^{1 / 2}}. \end{equation}$$

In approach (b), we simply calculate the temporal average (denoted 'tavg') and reform the result to a two-dimensional matrix Y of size m by p (where p corresponds to the number of remaining grid cells, which is the same for all datasets). We then calculate pairwise Euclidean distances δ_ij in between datasets i and j, resulting in the difference matrix D_tavg (dimensioned m by m):

$$\begin{equation} {\delta _{ij}} = {\left\{ \mathop \sum \limits^p_{l = 1} {\left[ {{\bf{Y}}\left( {i,l} \right) - {\bf{Y}}\left( {j,l} \right)} \right]^2}\right\} ^{1 / 2}}. \end{equation}$$

The distance matrices D_svd and D_tavg are then subjected to hierarchical clustering by applying the hclust function in R statistics software with Ward's clustering criterion (Ward 1963, ward.D2). The results are visualized in R using the Complex Heatmap package (version 1.14.0, available at https://bioconductor.org/packages/release/bioc/html/ComplexHeatmap.html). Clusters arevisually highlighted using a fixed threshold corresponding to the 99th percentile of the distances to enhance comparability of the results. To compare these clusters against the influence of the assessed factors (i.e.the data category, forcing, number of soil layers, ET scheme and model choice) on overall ET variability, we also assess the fraction of variation explained by each of those factors. This is achieved by means of applying a permutational multivariate analysis of variance (PERMANOVA, Anderson 2001) on the Euclidean differences D_svd and D_tavg.

Figure 1 shows global patterns of time averaged total ET (i.e. excluding crop models) grouped by data category. Spatial patterns of the ensemble averages are mostly similar among data categories and delineate the average location of the governing hydroclimatic regimes, such as the Intertropical Convergence Zone. The ensemble spread (as expressed by the IQR, i.e. the difference between the 75th and 25th percentile) in tropical rainforest regions is generally largest for the ISIMIP simulations and LFE ALL, whereas the uncertainties (in absolute values) in this region are considerably less pronounced in the diagnostic and E2O ensemble. Similarly, the spatial patterns of the relative ensemble spread (shown by means of the quartile coefficient of dispersion, QCD, which is a relative measure of dispersion, Bonett 2006) in the diagnostic ensemble are very well represented by both ISIMIP sectors and by LFE ALL. Although the magnitude of the uncertainties in the diagnostic ensemble is smaller than in the ISIMIP sectors, it is still higher than the relative ensemble spread in E2O. Note, however, that these small inter-model differences are linked to the fact that the E2O models are all forced by WFDEI, while the ISIMIP ensembles consist of three different forcings.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results of the SVD-based cluster analyses are shown in figures 2 (ET_tot) and figure 3 (ET_maize). As the associated clusters explain most of the variability among the assessed datasets (i.e. incorporating variability in both their spatialand temporal domains), we treat those as our main results, but occasionally draw comparisons to the clustering results based on time-averaged ET (visualized in figures A3 and figure A4). In the following, we highlight important details of the cluster diagrams by discussing individual parameters potentially having a measurable influence on the differences among the analysed ensemble.

At first glance, we notice a few apparent outliers in the displayed difference matrices, most strikingly for WB-MTE (ET_tot), PM-MU CSIRO, MPI-HM and the crop models LPJ-GUESS and ORCHIDEE-CROP (ET_maize only). As the differences among the contributing models are affected not only by the modelling structure but also the calibration process taken, it is not surprising to see such apparent differences between the non-calibrated models LPJ-GUESS and ORCHIDEE-CROP and all other datasets (the less pronounced signal for ORCHIDEE-CROP is potentially due to the fact that this model includes Nitrogen cycling). We also see that the long-term average global land mean ET of these models is not an outlier, suggesting that the differences are mainly due to anomalous patterns in the spatial domain. This is also supported by the clustering results of time-averaged ET, where outliers in the dendrogram mostly coincide with anomalies in the long-term mean (e.g. WANG-ET for ET_tot, figure A3 or PEGASUS for ET_maize, figure A4).

The meteorological forcing has, in general, a very noticeable impact on ET estimates when considering variability across both space and time. There are a number of clusters whose members are almost exclusively based on the same forcing (e.g. all PGMFD-forced simulations form a single cluster for ET_tot; cluster 2 in figure 2), underlining the dominant role of this parameter. However, the forcing apparently only has very minor influence on the spatial variability of time-averaged ET_tot (figure A3), indicating that most of the differences in between the employed forcing datasets are due to differences in the temporal evolution of ET. For ET_maize, there are a few more models for which the forcing dominates the differences, arguably due to a higher sensitivity of these models to diurnal forcings and other parameters which play a more important role when considering crop-specific ET aggregated over growing seasons.

Similarities among the individual datasets can also be reasonably well explained by their data category. While ET simulations from the ISIMIP global water and biomes models are very similar, crop models show substantial differences to all other realizations (the majority of the crop model simulations are members of the same cluster; cluster 3 in figure 3), reflecting the missing (or only rudimentary) representation of crops in the water and biomes models. E2O simulations also share most of their space-time variability (all but one of the E2O simulations fall into the same cluster; cluster 4 in figure 2 and cluster 5 in figure 3). This could be due to stronger similarities among the models participating in this project. Diagnostic datasets mostly show similarities in the spatial variability of the time averages (for ET_tot in figure A3, all but two diagnostic datasets are associated with the same cluster). It is also noticeable that ET diagnostics from the LandFlux project (i.e. PM-MU LANDFLUX, PT-FI LANDFLUX and SEBS LANDFLUX) form their own cluster in the spatio-temporal dendrogram (most pronounced in cluster 7, figure 2). For ET_tot, we also note that the LFE products are quite similar to most of the diagnostic datasets (sub-tree of cluster 1 in figure 2), adding support that LFE can indeed be used as a reference.

The other two parameters displayed in the cluster diagrams are the (binned) number of soil layers (water and crop models only) and the ET parametrization scheme. These parameters apparently only play a minor role for the clustering. While there are still a couple of instances where the ET scheme seems to contribute to small differences among the datasets (e.g. sub-trees of clusters 1, 2 and 4 in figure 2 indicate that datasets using the Penman-Monteith formulation show particularly little differences), there is no apparent influence of the number of soil layers. It must also be noted that there are many cases where both the ET scheme and the model are the same (e.g. cluster 5 in figure 2), suggesting that (in these cases at least) the associated small differences are due to the identity of the model. These findings are in contrast to what we would have expected from the literature, which suggests that ET schemes may explain a substantial part of the overall variance in ET (see section 1), whereas our results rather suggest that the forcing dominates this variance.

Besides looking at the individual clusters, it is also important to assess the fractional variance in ET for each of the factors considered (figure 4). In all cases, the assessed factors are capable of explaining more than 90% of the variabilities in the distance matrices (D_svd and D_tavg), and the identified proportions of explained variance are all statistically significant (p < 0.01). Irrespective of the distance metric and the type of ET estimate, differences in the forcing, model choice and ET scheme together account for more than half of the variabilities. In fact, the combined effect of model choice and ET schemes alone can explain at least 48% of the variabilities (which is again in line with previous findings that stress the importance of ET schemes). When considering D_svd, the forcing plays an important (for ET_tot even the leading) role, whereas differences in D_tavg can mainly be explained by the choice of the model. In spite of its minor role for the cluster analysis, the ET scheme can explain more than a quarter of the variabilities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now have a good overview of the structural differences among the assessed ensemble of datasets. Based on this information (using the largest sub-ensembles shown in figures 2 and figure 3), we can now go ahead and investigate differences among the clusters with respect to both the temporal evolution of ET and the representation of extreme events.

Zoom In Zoom Out Reset image size

Figures 5 and figure 6 present global and regional time series of both ET_tot and ET_maize, respectively. Among the diagnostics- and reanalysis-based estimates we see a few noticeable outliers for both ET_tot and ET_maize, most notably SEBS LANDFLUX (for which the global and regional averages exceed the IQR of the entire time span). The magnitude of the spread (IQR) among those estimates is comparable to the magnitude of the spread of the displayed cluster-based sub-ensembles. In the global domain, the diagnostics- and reanalysis-based estimates of ET_tot are slightly positively biased with respect to the other sub-ensembles shown. Those differences are most pronounced for cluster 2 (which is also located well below the other cluster ensembles and is also offset with respect to the median of LFE ALL), suggesting an underestimation of global land ET_tot, which is potentially related to the dominance of PGMFD-forced model simulations in this cluster. However, the respective sub-ensemble for ET_maize (cluster 4 in figure 6) is well embedded in the overall spread for all regions considered. We also observe that crop models from the ISIMIP agriculture sector tend to disagree more on ET_maize than datasets from the other clusters, although their median is well in agreement with the medians of the other sub-ensembles (cluster 3 in figure 6; note that this ensemble consists of only one model prior to 1979). Except from this, the medians of the clusters agree reasonably well, both in terms of their absolute magnitude (which is close to the median of LFE ALL) and their temporal variability.

The prominent jump in ET_tot from 1978–1979 apparent for cluster 1 within the Central European and global domains (figure 5) is due to the switch in the forcing data from WATCH to WFDEI at this time (when considering a sub-ensemble of only WFDEI-forced simulations, this discontinuity becomes even more pronounced, not shown). This artificial change signal is caused by differences in the reanalysis product used to create these datasets (ERA-40 for WFD, ERA-Interim for WFDEI), and is mainly due to the revision in the average aerosol loadings in North Africa and Europe affecting the short-wave incoming solar radiation in these regions (Weedon et al 2014).

A particularly important aspect to assess in greater detail is the response of the different datasets and sub-ensembles to drought events in the study areas (drought years indicated in rows 2–4 in figures 5 and figure 6; in addition, monthly time series for each event year are shown in figure 7). Irrespective of the region and event considered, month-to-month changes in ET anomalies in the displayed sub-ensembles are similar to the changes in the ensemble mean of the diagnostic estimates. However, we note that there is a tendency of the model-dominated sub-ensembles to exaggerate negative ET anomalies (relative to the diagnostic estimates), which could be related to overly strong soil water depletion in some of the model simulations, to models overestimating stomatal closure, to some models not representing irrigation, or to missing non-stomatal flux components.

For each of the considered regions, there has been at least one pronounced drought during the analysed time span (see table 2). While the 1987 Great Plains drought has no apparent influence on ET rates, the events in 2002 and 2012 are characterized by negative ET anomalies apparent in both the displayed clusters and in most of the diagnostic estimates (the anomaly signal is slightly more pronounced for ET_maize). In the case of western Russia, negative ET anomalies accompanying the 2010 heat wave are found in most of the diagnostic datasets, in the reanalyses, and most strikingly in the (model-dominated) clusters, which is in line with the ET signal discussed in Hauser et al (2016). For the European heatwave in 2003, the diagnostic estimates show a noticeable disagreement from the models and reanalyses. However, the month-to-month variations in ET anomalies from most of the datasets and simulations (figure 7) still resemble the event. The patterns agree both with theory Seneviratne (2012) and catchment-level observations Teuling et al (2013), as anomalously dry weather (accompanied by high net radiation and anomalously warm temperatures) commonly leads to an initial increase in ET, followed by a decrease later in the season (when soil moisture reaches critical levels). These temporal patterns can partly also be observed for the other events (although the initial positive ET anomalies are less pronounced or located earlier in the season).