Introduction

Ecosystem climate change experiments are one of the key instruments to study the response of ecosystems to a change in climate. There are primarily four different factors that are altered in such experiments: air temperature, precipitation, CO2 concentration, and nitrogen deposition (Curtis and Wang 1998; Rustad et al. 2001; Lin et al. 2010; Wu et al. 2011; Knapp et al. 2018). More recently, multi-factor experiments are starting to emerge. In those experiments, different combinations of the four main drivers are altered (Kardol et al. 2012; Yue et al. 2017). What is common in the majority of climate change experiments is that while the drivers of interest are being altered, all other variables are being held equal between the different treatment groups. Consequently, differences in the response can be related to the change in the main driving factor (or multiple driving factors).

Altering only one or a limited number of climate change drivers allows for a straightforward analysis of the observed responses and has provided a wealth of mechanistic insights in ecosystem responses to environmental changes (e.g., Hovenden et al. 2014; Karlowsky et al. (Karlowsky et al. 2018; Terrer et al. 2018)). However, the resulting multivariate combination of climate variables may be physically unrealistic and may miss key aspects related to natural climate variability and the co-variance of multiple variables, linked to each other by synoptic conditions. This is particularly important for representing compound events, where the combination of non extreme drivers can lead to extreme impacts (Zscheischler and Seneviratne 2017; Zscheischler et al. 2018; Rineau et al. 2019). For example, droughts and heatwaves often co-occur (Zscheischler and Seneviratne 2017), and soil moisture conditions and precipitation occurrence are linked (Guillod et al. 2015; Moon et al. 2019). Incorporating the co-variability of key climate drivers is also important for the studied responses. For instance, heatwaves characterized by similar extreme air temperature can lead to different plant responses depending on the atmospheric conditions: under different shortwave radiation, relative humidity, and surface wind conditions, the leaf temperature and the potential for heat stress vary a lot (De Boeck et al. 2016).

Until recently, it was not possible to simulate plausible future climates in ecosystem climate change experiments (Korell et al. 2019), as these experiments require accurate manipulation of environmental variables to represent current and future climate conditions. Controlled environment facilities meet these requirements by providing systems to simultaneously manipulate as well as measure multiple parameters (e.g., Lawton 1993; Lawton 1996; Griffin et al. 1996; Steward et al. 2013; Clobert et al. 2018), especially in combination with an observation station in the field providing real-time observations of most of those parameters (Rineau et al. 2019). This approach is powerful especially when combined with a measurement station in the field providing real-time observations of most of these required parameters (Rineau et al. 2019). In such facilities, climate change experiments can be informed by meteorological forcing representing both present and future climatic conditions in a holistic manner. For instance, this forcing can include both realistic changes of climate variability as well as important drivers of changes in the frequency, intensity, and duration of meteorological extremes. This potential is especially interesting in gradient experiments covering a range of global warming levels, as this combination allows for the detection of non-linearities, thresholds, and possible tipping points in ecosystem responses to increasing climate change forcing (Kreyling et al. 2018; Rineau et al. 2019).

Sampling realistic climate information in a climate change context is challenging but can be achieved by using climate model output. Global climate models (GCMs) are generally used to assess the climate state and variability at global to continental scales with a resolution of 100 to 250 km. By dynamically downscaling GCMs, regional climate models (RCMs) typically resolve the climate on a regional scale with higher spatial resolutions of 1 to 50 km. As such, RCMs allow a more realistic representation of meso-scale atmospheric processes and processes related to orography and surface heterogeneities. As climate models realistically simulate the atmospheric state under past, present, and future climatic conditions with a high temporal resolution, they are suited to provide a holistic and physically consistent climate forcing for ecosystem climate change experiments. Generally, ensemble climate projections show a large spread for future climate conditions (Keuler et al. 2016), especially for variables relevant for ecosystem experiments such as extreme air temperature, droughts, and intense precipitation (Sillmann et al. 2013; Orlowsky and Seneviratne 2013; Rajczak and Schär 2017; Greve et al. 2018). This spread is related to (i) different climate sensitivities of the GCMs, (ii) structural differences between the models, and (iii) natural variability within the climate system. The Coordinated Regional Climate Downscaling Experiment in the European domain (EURO-CORDEX) provides an ensemble of high-resolution dynamically downscaled RCMs (Kotlarski et al. 2014) and is therefore highly suitable to serve as a base for the selection of representative climate forcing for climate change experiments. With a suite of GCM/RCM combinations available, a well-informed choice on the most adequate RCM/GCM simulation can be made based on (i) the model skill in representing the observed climatology and (ii) the air temperature sensitivity to future increases in greenhouse gas concentrations.

So far, statistically downscaled GCM output has only rarely been used as climate forcing in ecosystem experiments. Thompson et al. (2013) describe a process for generating air temperature forcing for experiments in which they use daily air temperature output from a GCM (MIROC) and a stochastic weather generator to generate hourly weather. They validated their method against statistical characteristics of air temperature observations. Likewise, the Montpellier CNRS ecotron facility is driven by multivariate statistically downscaled GCM projections (using the ARPEGEv4 model; Roy et al. 2016). They force their experiment with climatic conditions of an average climatological year of the period 2040–2060. During the summer months, they artificially simulate an extreme event by including drought and heatwave by reducing the irrigation amount to zero and increasing the air temperature artificially. However, by using a climatological year, possible extreme events are dampened by averaging. Both studies lack a thorough evaluation procedure for selecting the used climate model. Moreover, to the best of our knowledge, no study accounts for the co-variance between climate variables.

In this paper, we present a new method for creating realistic climate forcing for manipulation experiments. From an ensemble of dynamically downscaled climate model simulations, we select one simulation that well represents present-day climate conditions for four key variables in the region of interest and is representative of the multi-model mean of these variables in future projections. In this way, the new methodology accounts both for co-variance of climate parameters and for climate variability while naturally incorporating extreme events under present and future climate conditions. Furthermore, the method can be combined with a gradient approach. We apply the new methodology to generate climate forcing for the UHasselt Ecotron Experiment, an infrastructure consisting of 12 climate-controlled units, each equipped with a lysimeter containing a dry heathland soil monolith extracted from the National Park Hoge Kempen in Belgium (Rineau et al. 2019). In this experiment, six units are directly forced with regional climate model output along a global mean air temperature (GMT) gradient anomaly.

Data and methods

New methodology for generating climate forcing for ecosystem climate change experiments

In our methodology, variability and co-variance between variables are preserved by selecting the best-performing RCM simulation and subsequently extract the required variables from the grid cell covering the location of the experiment. By extracting a single grid cell of a single RCM simulation, climate extremes are not smoothed and the climate variability inherent to the model is fully preserved. The units in the ecosystem climate change experiments follow a gradient of increasing GMT anomalies. In this way, a given unit is forced with the climatic conditions consistent with, e.g., a 2∘C warmer world, and the units represent conditions associated with increasingly warmer climates.

The methodology presented here is deployed in three steps. First, the best-performing RCM projection needs to be selected based on two criteria: (i) the simulation should have high skill in reproducing mean and extreme present-day climatic conditions and (ii) the projected future air temperature anomalies should be close to the multi-model mean; that is, the selected simulation should be representative of the future mean projection (Fig. 1, step 1). To this end, the model performance is evaluated for four variables that are highly relevant for ecosystem climate change experiments: precipitation, air temperature, relative humidity, and surface wind speed. Precipitation is considered one of the most important variables, as water availability is likely to constrain plant growth the most.

Fig. 1
figure 1

Methodology for generating climate forcing along the GMT anomaly gradient

Full size image

Second, the time windows for the different units along the GMT anomaly gradient are defined based on the annual GMT projection of the driving GCM of the chosen RCM simulation (Fig. 1, step 2). To span a large range of climate change scenarios, we use projections following the Representative Concentration Pathway (RCP) 8.5, a worst-case scenario following an unabated greenhouse gas emissions pathway (Riahi et al. 2011). The UHasselt Ecotron experiment, including all units, is running for 5 years. We choose time windows corresponding to the experimental period and centered around the year in which the climatological GMT anomaly (averaged with a 30-year period) crossed the predefined thresholds for the first time. In the third step, the values of all necessary variables are extracted from the chosen RCM projection based on the defined time windows for the grid cell covering the experiment location (Fig. 1, step 3). These time series are then directly used to force the ecotron units, in the highest available temporal resolution.

The UHasselt Ecotron experiment

The UHasselt Ecotron experiment is an ecotron infrastructure consisting of replicated experimental units in which ecosystems are confined in enclosures. By allowing the simultaneous control of environmental conditions and the online measurement of ecosystem processes, the ecotron units are suited for experiments with highly controlled climate change manipulation of large intact parts of the ecosystem. The infrastructure allows intensive monitoring and control of key abiotic parameters on 12 large-scale ecosystem replicas, called “macrocosms.” These macrocosms had been extracted without disruption nor reconstitution of the soil structure from the same dry 6- to 8-year-old heathland plot in the National Park Hoge Kempen (50∘ 59’ 02.1” N, 5∘ 37’ 40.0” E) in November 2016.

The infrastructure is a W-E oriented, 100 m by 10 m wide, and 6 m tall building (Fig. 2a). Only 12 of the 14 units are used, excluding the outermost to avoid boundary effects. Each unit consists of three compartments in which the abiotic environmental variables are controlled: the dome, the macrocosm, and the chamber. The dome is transparent for photosynthetic active radiation (PAR) and long- and medium-wave ultraviolet radiation (UVa and UVb, respectively). Here, wind and precipitation are measured and generated, and CO2, N2O, CH4, PAR, and net radiation (NR; i.e., the difference in incoming and outgoing short-and longwave radiation) are measured. The second compartment, the macrocosm, contains the extracted soil column (the ecosystem) enclosed in a lysimeter. In this compartment, the soil water content, soil water tension, soil electrical conductivity, and soil temperature are measured and controlled. The chamber, the third compartment, the air pressure, air temperature, relative humidity, and CO2 concentration are controlled. The ecotron infrastructure is linked with an Integrated Carbon Observation System (ICOS) ecosystem station, which provides real-time information on local weather and soil conditions. These data are used to simulate the current weather conditions within the ecotron units with a frequency of at least once every 30 min (Rineau et al. 2019).

Fig. 2
figure 2

The UHasselt Ecotron experiment. a (picture: Liesbeth Driessen). b Scheme of a unit with the three compartments (1) denoting the dome; (2) the lysimeter, shown in detail on the right; and (3) the chamber. c An overview map with location of the infrastructure and reference weather observation stations

Full size image

The aim of the UHasselt Ecotron experiment is to study the ecological and societal impacts of climate change, by manipulating climatic variables alone or in combination and, across a wide range of predicted values, while monitoring as many soil biota and processes as possible and to translate them into socio-economic values using heathland as a case study (Rineau et al. 2019). Examples of measured ecosystem processes are evapotranspiration, net ecosystem exchange, CH4 or N2O emissions. The main research questions of this multi-disciplinary experiment are how climate change will affect the transitioning of the heathland ecosystem to alternative stable states like pine forest or acid grassland and what the consequences are for ecosystem services. The experiment will run uninterrupted for a period of at least 5 years. Six units will be used to simulate a gradient of increasing variability in precipitation regime. They are driven by the ICOS station and a perturbed precipitation time series following a gradient of increasingly long periods with no precipitation (2, 6, 11, 23, 45, and 90 days). In the remaining six units, atmospheric conditions along the GMT anomaly gradient will be simulated as described in “New methodology for generating climate forcing for ecosystem climate change experiments.” The 3-hourly RCM output is linearly interpolated to a 30-min time resolution to force the ecotron units. For soil temperature and soil water tension however, the 30-min ICOS data is used. This is because leaving the lysimeter uncontrolled would lead to (i) an overestimation of soil temperature variability as the lysimeter is exposed to air temperature in the chamber (despite being thermically insulated), and (ii) accumulation of water at the bottom of the lysimeter, hence considerably overestimating soil water level, as soil water movements are mimicked by suction from the bottom. Following the gradient design, each ecotron unit represents the local climate conditions of a globally 0∘C (historical), + 1 ∘C (present day), + 1.5 ∘C (Paris Agreement), + 2 ∘C, + 3 ∘C, and + 4 ∘C warmer world. The climatology of the unit forced by + 1 ∘C can thereby be directly compared to the unit driven by the ICOS station and thus representing the present-day observed conditions. In this regression design, there is no experiment replication. To minimize the noise in initial ecosystem responses, the units are allocated to the two gradient experiments based on a cluster analysis of the variance of the 14 variables measured during a test period of 11 months (Rineau et al. 2019).

Meteorological data

EURO-CORDEX

The best-performing RCM simulation compared to observations is selected from the Coordinated Regional Climate Downscaling Experiment in the European domain (EURO-CORDEX), an ensemble of high-resolution dynamically downscaled simulations available at a horizontal resolution of 12 km (0.11∘ on a rotated grid; Jacob et al. 2014; Kotlarski 2014). The simulations, hereafter referred to as GCM downscalings, cover the historical period (1951–2005) and the three RCP scenarios (RCP 2.6, 4.5, and 8.5, for the period 2006–2100) by using GCMs as initial and lateral boundary conditions. Additionally, for each RCM, a reanalysis downscaling is provided in which the RCM is driven by the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim as initial and lateral boundary conditions for the period 1990–2008 (hereafter referred to as reanalysis downscalings). These reanalysis-driven simulations allow to evaluate the skill of the RCMs themselves by comparing them to observations (Kotlarski et al. 2014).

In this study, we use the variables for daily mean, minimum, and maximum air temperature, precipitation, mean surface wind, and relative humidity of all available simulations (Table 1). We consider the values of the 12 km by 12 km pixel covering the location of the reference station providing the observations. As relative humidity is not directly available for all simulations, we converted specific humidity to relative humidity using the mean air temperature and surface air pressure for every simulation. Comparing the applied conversion with the simulations for which relative humidity is available proves this conversion is applicable. Neither specific nor relative humidity is publicly available for the simulations with RegCM4-2 and ALARO-0 and the mean surface wind speed variable is not available for ALADIN53 and ALARO-0; therefore, we do not analyze these variables for the respective simulations.

Table 1 Bias in annual precipitation (P bias) and rank based thereof (from 1-best to 18-worst) for the EURO-CORDEX GCM downscalings for the period 1951–2005 over Maastricht Airport
Full size table

Once the EURO-CORDEX ensemble member is selected, the relevant variables (precipitation, mean air temperature, surface air pressure, surface up-welling latent heat flux and sensible heat flux, wind speed, and relative humidity) are extracted from the 3-hourly RCP 8.5 simulation for the pixel covering the ecotron location for the time windows in which the GMT anomalies are crossed for each dome. These 3-hourly values (except for surface up-welling latent heat flux and sensible heat flux) are then linearly interpolated to 30-min resolution and used to drive the climate controllers in the ecotron units. For precipitation, one additional step was added where drizzle (precipitation of less than 1 mm) was postponed and accumulated until it reached 1 mm to start a rain event in the ecotron. The surface air pressure is calculated from the mean sea level pressure using the elevation of the ecotron facility (43 m a.s.l.) and assuming hydrostatic equilibrium. The concentrations of the controllable greenhouse gases (CO2, CH4, and N2O) are determined based on the annual values calculated by van Vuuren et al. (2011) according to RCP 8.5. These correspond to the prescribed concentrations of the RCM simulations.

Weather station observations

Reference station data is obtained from the European Climate Assessment and Dataset (Klein Tank et al. 2002). The three operational weather stations closest to the UHasselt Ecotron experiment are Maastricht Airport (11 km), Aachen (37 km) and Heinsberg-Schleiden (29 km; Fig. 2c). These weather stations provide daily observations from the end of the 19th century (Maastricht Airport and Aachen) or mid 20th century (Heinsberg-Schleiden) until the present-day, thereby covering both the EURO-CORDEX GCM and reanalysis downscaling periods. All stations record air temperature (∘C), precipitation (mm day− 1), relative humidity (%) and surface wind speed (m s− 1) at daily resolution, except for the Heinsberg-Schleiden station where there are no surface wind observations available.

The seasonal cycles of the observations for the different stations follow a similar annual course (Fig. 3). For air temperature, the curves overlay and for precipitation they are similar. Relative humidity (RH) has a small offset between the three stations, possibly owing to the differences in absolute height and local topography. The difference in surface wind speed between Maastricht Airport and Aachen is considerable but plausible considering the large spatial variability in wind speed. Given that the model evaluation showed very little sensitivity to the choice of the reference station, we hereafter present the results with the reference station closest to the ecotron facility (Maastricht Airport).

Fig. 3
figure 3

Seasonal cycles of observed mean air temperature (a), mean daily precipitation (b), mean relative humidity (c), and mean surface wind speed (d) in the weather stations of Maastricht Airport, Aachen, and Heinsberg-Schleiden (monthly averages based on daily data from 1963 to 2018). For Heinsberg-Schleiden no surface wind observations are available. The curves for air temperature are overlaying

Full size image

Metrics and diagnostics

The evaluation of the EURO-CORDEX ensemble members is performed using different metrics accounting for performance of representing the climatic means, distributions, and extremes.

A ranking is made of the reanalysis downscalings, assigning the lowest ranks to the best-performing models and higher ranks to the least-performing models (1-best to 9-worst). First, the bias is calculated as the difference between the averages of the daily modelled and observed variables. The second metric, the Perkins skill score (PSS), is a quantitative measure of how well each simulation resembles the observed probability density functions by measuring the common area between two probability density functions (Perkins et al. 2007). The mean absolute error (MAE) is calculated by taking the means of the absolute differences between the modelled and observed seasonal cycles, calculated based on the whole series. This is done for the whole series and to capture the potential errors in the extremes, also for the 1st, 10th, 90th, and 99th percentiles which are calculated based on the daily time series of both observed and modelled time series. Next, the root mean square error (RMSE) is calculated by taking the root of the squared errors. The Spearman rank correlation (hereafter referred to as Spearman) coefficient shows the correlation of the observed and modelled series, calculated based on daily values. Finally, the Brier skill score (BSS) is calculated, which gives an indication of the improvement of the Brier score (an index to validate probability forecasts) compared to a background climatology in which each event has an equal occurrence probability (Brier 1950; Murphy 1973). For the GCM downscalings, we use the same ranking method and scores, except for the RMSE, Spearman rank correlation, and BSS because the internal variability, inherent to individual simulations with a coupled climate model, cannot be predicted on multi-decal timescales, and can therefore not be compared to observations on a day-by-day basis (Fischer et al. 2014; Meehl et al. 2014).

In addition to the performance metrics computed on the actual time series, the RCM performance is also evaluated based on the bias in climatological diagnostics related to air temperature and precipitation. To this extent, the average diurnal air temperature range (DTR; K; the difference between the daily maximum and minimum air temperature) is calculated for the whole year, for the winter (December–February) and summer (June–August) season. Next, the number of wet days (defined as days during the year for which precipitation is larger than 0.1 mm or larger than 1 mm, to account for differences in drizzle; (Casanueva et al. 2016) and the number of frost days (days with a minimum air temperature below 0∘C) are calculated. Furthermore, the monthly maximum 1-day precipitation (Rx1day; mm day− 1) and the number of consecutive dry days (CDD; days); the annual maximum number of days for which precipitation is below 1 mm and consecutive wet days (CWD; days); and the annual maximum number of days for which precipitation is equal to or more than 1 mm are included in the analysis. All indices are calculated for the simulated and observed time series, and consequently the ranking is established based on the difference between the model and observed diagnostic. Next, the correlation between the different variables is evaluated by comparing them to the observed correlation. This is done both on the annual time scale and for the summer and winter seasonal averages, as correlations are expected to differ in sign and magnitude between the two seasons (e.g., negative correlation between air temperature and relative humidity in summer reflecting heatwave conditions, and a positive correlation between wind speed and precipitation in winter reflecting storm conditions).

After choosing the best-performing simulation based on the evaluation of both the reanalysis and GCM downscalings, the climate change signals for this simulation are investigated by calculating changes in various climate change indices, based on the Expert Team on Climate Change Detection and Indices (ETCCDI; see http://etccdi.pacificclimate.org/list_27_indices.shtml) for the 5-year periods defined by the GMT anomalies relative to the reference period (1951–1955). These indices are widely used for analyzing changes in extremes (e.g., Zhang et al. 2009; Orlowsky and Seneviratne 2013; Sillman et al. 2013). The air temperature indices are (i) Δ T (∘C), the mean daily air temperature change; (ii) Δ TXx (∘C), the difference in the annual maximum value of daily maximum air temperature; (iii) Δ TNn (∘C), the difference in the annual minimum value of daily minimum air temperature; (iv) Δ frost days, the difference in the number of frost days (with a minimum air temperature below 0∘C); (v) Δ summer days, the difference in the number of summer days (with the maximum air temperature above 25∘C); and finally (vi) Δ GSL (days), the difference in growing season length, defined as the annual count between the first span of at least 6 days with a daily mean air temperature higher than 5∘C and the first span after July 1st of 6 days with a daily mean air temperature lower than 5∘C. The precipitation indices are (i) Δ PRCPTOT (mm), the difference in annual accumulated precipitation (as simulated over the 5-year period); (ii) Δ Rx1day (mm), the difference in monthly maximum 1-day precipitation; (iii) Δ R10mm (days), the difference in the number of days per year with more than 10 mm precipitation, (iv) Δ CDD (days), the difference in the maximum length of a dry spell (measured as the maximum number of consecutive days with less than 1 mm precipitation); and finally, (v) Δ CWD (days), the maximum length of a wet spell (measured as the maximum number of consecutive days with more than 1 mm precipitation).

Applying the new methodology for the UHasselt Ecotron experiment

The best-performing RCM simulation is identified by elimination based on expert judgment based on the performance of the two selection criteria. Next, we define the time windows for the different units along the gradient based on the 30-year averaged GMT anomaly of the driving GCM under RCP8.5 relative to 1951–1955. Based on these time windows, we extract the 3-hourly data for all necessary variables from the simulation for the 12-km by 12-km grid cell covering the location of the experiment.

Results

Identification of the best-performing model simulation

First criterion: skill in present-day climate

Overall, model skill strongly varies across RCMs (Fig. 4). While the annual air temperature cycle is generally well represented by all RCMs, biases may reach up to 2 degrees in individual months for some RCMs. The biases in precipitation are generally positive (up to factor 2.4) and vary across RCMs. Only CCLM4-8-17 simulates precipitation in the same range as the observed climatology (nearly no bias (100.22 %) on annual mean precipitation amounts), while the other RCMs overestimate the total precipitation amounts from 114 % up to 182 %. For relative humidity and surface wind speed, all RCMs generally succeed in representing the seasonal cycle, but exhibit deviations in amplitude and absolute values (e.g., amplitude biases of RCA4 (− 37.8 %), ALADIN53 (23.3 %) and CCLM4-8-17 (+ 16.3 %) for relative humidity, and annual mean biases for WRF331F (+ 15.6 %) and HIRHAM5 (− 9.1 %) for surface wind speed). Overall, these seasonal cycles indicate that for all simulations, the relative bias in precipitation is large compared to biases in other variables.

Fig. 4
figure 4

Seasonal cycle of the reanalysis downscalings for mean air temperature (a), mean daily precipitation (b), mean relative humidity (c), and mean surface wind speed (d). (The RegCM4-2 and ALARO-0 simulations are not available for relative humidity and the ALADIN53 and ALARO-0 simulations are not available for surface wind speed)

Full size image

The rankings of the reanalysis downscalings for the four variables (Fig. 5) indicate that, overall, CCLM4-8-17, RACMO22E, REMO2009, and HIRHAM5 are performing best. CCLM4-8-17 and RACMO22E show the highest relative skill for precipitation, while REMO2009 and HIRHAM5 demonstrate high skill for air temperature. CCLM4-8-17 is the best-performing model based on the bias and total MAE metrics for air temperature and precipitation but is ranked in the mid range for the metrics related to the shape of its air temperature distribution (PSS and percentile MAE). This can be attributed to an overestimation of the amplitude of the seasonal air temperature cycle in this model (too cold in winters, too hot in summers; Fig. 4a; Kotlarski et al. 2014). For relative humidity and surface wind speed, RACMO22E generally demonstrates the highest skill. Considering the climatological diagnostics (Fig. 6a), CCLM4-8-17 shows the highest relative skill for precipitation-related diagnostics (wet days, monthly maximum 1-day precipitation, length of dry and wet spells), while RACMO22E and RCA4 show higher relative skill for the annual, winter, and summer diurnal air temperature range. While RCA4 is highly ranked for air temperature-related diagnostics, it is one of the models with the lowest relative skill for precipitation-related diagnostics. The correlation ranking shows a more scattered image, for the annual correlation as well as summer and winter correlations (see Suppl. Fig. A1). Overall, as the reanalysis-driven simulations with ALADIN53, RegCM4-2, WRF331F and ALARO-0 show the lowest skill compared to the other RCMs, we take them out of consideration to serve as ecosystem forcing.

Fig. 5
figure 5

Ranking of the reanalysis downscalings based on performance on mean air temperature (a), mean daily precipitation (b), mean relative humidity (c), and mean surface wind speed (d) compared to observations from Maastricht. The metrics shown are the bias, Perkins skill score (PSS), mean absolute error (MAE) for the entire time series and the 1st, 10th, 90th, and 99th percentiles, root mean square error (RMSE), Spearman rank correlation (Spearman), and Brier skill score (BSS). Rankings are from 1-best to 9-worst. Gray colors indicate that the variable is not available for the considered model

Full size image
Fig. 6
figure 6

Ranking of the reanalysis (a) and GCM (b) downscaling for the historical period based on climatological diagnostics. Diurnal air temperature range (DTR) in summer (July–August) and winter (December–February), number of wet days defined as days with precipitation > 0.1 mm and precipitation > 1 mm, number of frost days defined as days with mean air temperature < 0 ∘ C, monthly maximum 1-day precipitation (Rx1day), consecutive dry days (CDD), the maximum length of a dry spell, and consecutive wet days (CWD), the maximum length of a wet spell. Next to the diagnostic name, its value as observed in Maastricht Airport is shown. Rankings are from 1-best to 9 or 18-worst for the reanalysis and GCM downscaling, respectively

Full size image

Second, we evaluate the GCM downscalings for the period 1951–2005. The seasonal cycles of the air temperature, precipitation, relative humidity, and surface wind speed show a similar pattern as the reanalysis downscalings, with again a strong wet bias for precipitation in most models (see Suppl. Fig. A2). The rankings show a mixed pattern for the different variables: there are no simulations which rank high for all considered variables (Fig. 7). For precipitation, the simulations with CCLM4-8-17 and RACMO22E have better relative skill compared to the other simulations, which is in line with the high ranking of these models in the reanalysis downscalings. Furthermore, it is remarkable that the simulations which show a high skill for precipitation, typically show lower skill for relative humidity and vice versa, e.g., CCLM4-8-17 driven by HadGEM2-ES (high ranking in precipitation, lowest in relative humidity) and REMO2009 driven by MPI-ESM-LR (high ranking in relative humidity and lower in precipitation). The three MPI-ESM-LR-driven simulations appear to be better in reproducing the air temperature climatology compared to the other simulations. For the climatological diagnostics, generally CCLM4-8-17 is scoring best for the precipitation-related diagnostics, whereas simulations with RCA4 are ranked the best for DTR (annual, summer, and winter).

Fig. 7
figure 7

Ranking of the GCM downscalings based on performance on mean air temperature (a), mean daily precipitation (b), mean relative humidity (c), and mean surface wind speed (d) compared to observations from Maastricht. The metrics showed are the bias, Perkins skill score (PSS), mean absolute error (MAE) for the total and 1st, 10th, 90th, and 99th percentile. Rankings are from 1-best to 16, 17, or 18-worst for surface wind speed, relative humidity, precipitation, and air temperature, respectively. Gray colors indicate that the variable is not available for the considered model

Full size image

Based on the ranking of the GCM downscalings, the following simulations are considered potential candidates to serve as climate forcing: CCLM4-8-17 driven by CNRM-CM5, EC-EARTH and MPI-ESM-LR, HIRHAM5 driven by EC-EARTH and HadGEM2-ES, and RACMO22E driven by HadGEM2-ES (Figs. 5, 7, and 6). Since precipitation biases strongly differ among RCMs (Table 1), and since precipitation is a critical variable for the ecosystem experiments (Van der Molen et al. 2011; Vicca et al. 2014; Estiarte et al. 2016), we prioritize a minimum relative bias for precipitation over a lower bias for air temperature, relative humidity, and surface wind speed. The precipitation biases for the considered simulations are + 150 mm year− 1 for CCLM4-8-17 driven by CNRM-CM5, + 8 mm year− 1 for CCLM4-8-17 driven by EC-EARTH, + 24 mm year− 1 for CCLM4-8-17 driven by MPI-ESM-LR, + 323 mm year− 1 for HIRHAM5 driven by EC-EARTH, 101 mm year− 1 for HIRHAM5 driven by HadGEM2-ES, and 36 mm year− 1 for RACMO22E driven by HadGEM2-ES. Based on this, the CCLM4-8-17 EC-EARTH-driven simulations has the best chance to be chosen as forcing, followed by the CCLM4-8-17 MPI-ESM-LR, and the RACMO22E HadGEM2-ES-driven simulation.

Second criterion: Representativeness of multi-model mean in future projections

To verify the second requirement, we look at anomalies from the mean signal of the four variables for the future period of the simulations under RCP 8.5. The EC-EARTH-driven CCLM4-8-17 simulation is representative of the multi-model mean for all four variables (Fig. 8), and even the median simulation for the mean air temperature anomaly. For precipitation and relative humidity, however, the CCLM4-8-17 EC-EARTH simulation shows decreasing anomalies after 2050 and underestimates the multi-model mean anomaly. The other selected simulations have a larger positive bias in precipitation for their GCM downscalings. A possible reason is that these simulations overestimate precipitation and simulate a more intensive hydrologic cycle, which also implies stronger changes in the future.

Fig. 8
figure 8

Anomalies for the CCLM4-8-17 EC-EARTH simulation following RCP 8.5 at the ecotron site for mean air temperature (a), mean daily precipitation (b), mean relative humidity (c), and mean surface wind speed (d). The reference period is 1977 to 2006, the anomalies of the CLM4-8-17 EC-EARTH simulation are calculated compared to its own values in the reference period. In gray, the envelope of all EURO-CORDEX RCP 8.5 simulations is showed

Full size image

The remaining five simulations from step 1 (CCLM4-8-17 driven by MPI-ESM-LR, HIRHAM5, and RACMO22E driven by HadGEM2-ES) all systematically underestimate or overestimate other variables (Suppl. Figs. A4, A5, A6, A7 and A8). For instance, the mean air temperature anomaly of CCLM4-8-17 driven by MPI-ESM-LR simulation (1.46∘C) is lower than the 10th percentile of all simulations (1.51∘C) and the air temperature anomaly for CCLM4-8-17 driven by CNRM-CM5 is the 30th percentile (1.67∘C). HIRHAM5 driven by HadGEM2-ES overestimates relative humidity anomalies compared to the multi-model mean, with a mean value (1.26 %) around the 80th percentile. Finally, the HadGEM2-ES-driven RACMO22E simulation overestimates relative humidity and air temperature anomalies, up to the 90th percentile for air temperature. Overall, we conclude that the EC-EARTH-driven CCLM4-8-17 simulation is the most appropriate candidate for serving as climate forcing for the UHasselt Ecotron experiment.

Characterization of the selected meteorological forcing

Based on the selection criteria, we single out the EC-EARTH (ensemble member r12i1p1) driven CCLM4-8-17 simulation as climate forcing for the UHasselt Ecotron experiment. The climatic conditions in the six units along the gradient represent an increasing signal of climate change, representing the local climatic conditions of 6 future climate states corresponding to an increasing GMT. The overall trend of the local air temperature anomaly compared to the reference period (0∘C) increases monotonically with the corresponding GMT anomalies (Fig. 9a). No clear trends are visible for precipitation, relative humidity, and surface wind speed anomalies, but very clear for the minimum and maximum air temperature anomalies which are both increasing (Fig. 9). The mean daily air temperature is increasing at a similar rate compared to GMT anomaly, and minimum and maximum air temperature show a larger increase (Table 2). None of the air temperature indices show a linear increase, reflecting the difference between global and local climatic conditions and the influence of decadal internal variability. The ecotron unit representing a + 4 ∘C world is the most extreme case, with increases of TXx of + 6.30 ∘C and an increase of TNn with + 10.21 ∘C (Table 2). The number of frost days decreases with about − 76.2, while the number of summer days with air temperature above 25∘C increases with about 36.6 days. The annual growing season length is extended with 80 days on average, leaving only 59.4 days of the year not favorable for growth. The indices for precipitation show a less clear trend (Table 2). The total precipitation amount varies for the five units, without any trend and shows a substantial decadal variability in all seasons (see Fig. 9). Rx1day has positive anomalies for the + 1.5 ∘C, + 2 ∘C, and + 3 ∘C units (+ 0.35 mm day− 1, + 1.92 mm day− 1, and + 2.34 mm day− 1, respectively). These + 2 ∘C and + 3 ∘C units also know an increase in R10mm (+ 3.2 and + 3.6 days) compared to the other units. Finally, there is no clear trend in CWD, but there is an increase in CDD up to + 11.8 days for the + 4 ∘C unit. The + 1.5 ∘C unit spans a drier time window, with an average CDD of + 9.6 days. Figure 9 further shows a systematic decrease of relative humidity during summer with increasing warming and a strong decadal variability of surface wind speed especially in winter.

Fig. 9
figure 9

Annual cycles of the CCLM4-8-17 EC-EARTH ecotron unit forcing for the + 1 ∘C, + 1.5 ∘C, + 2 ∘C, + 3 ∘C, and + 4 ∘C units compared to the 0∘C reference period with a daily mean air temperature, b mean daily precipitation, c mean relative humidity, d mean surface wind speed, e daily maximum air temperature, and f daily minimum air temperature. Curves were smoothed using Savitzky-Golay filtering (order = 2 frame = 301; Savitzky and Golay 1964)

Full size image
Table 2 Extracted 5-year periods and air temperature and precipitation indices based on ETCCDI for the CCLM4-8-17 EC-EARTH simulation at the ecotron location
Full size table

Discussion

The presented methodology exhibits some challenges, which are addressed in the following section.

We extract all climate variables from one grid cell of the RCM simulation to conserve the most realistic, non smoothed signal of the climate models. However, the extracted time series of the grid cell can differ a lot between different models and time periods, reflecting the natural climate variability. GCMs and RCMs provide robust signals when aggregated over a larger spatial area (Seneviratne et al. 2016; Fischer and Knutti 2015). By taking the spatial mean, a more robust estimate of the mean climate is obtained, including robust signals of climate change. This explains the difference in local climate change signals (Fig. 8; Table 2) and non-linearities compared to the GMT anomaly obtained by global averaging (Seneviratne et al. 2016). It is however necessary to use actual time series from a single grid cell to capture, e.g., the extreme precipitation event occurring in the considered grid cell, but not in the neighboring grid cells. The grid cell values also reflect strong inter-annual to decadal variability which is of high relevance for a realistic forcing of the ecosystem.

Climate model simulations are often biased, which is mostly related to structural model deficiencies (Flato et al. 2013). Applying bias adjustment is a standard way to deal with biases (Gudmundsson et al. 2012; Vanderkelen et al. 2018), but such methods face several challenges and need to be chosen carefully to not increase biases in the co-variability of variables (Zscheischler et al. 2019). In the proposed method, we therefore directly use the “raw” model output, as such preserving climate variability and the physically consistent co-variance of the different meteorological variables. In this way, the Ecotron experiment will study ecosystem responses to multivariate drivers as compound controls. For instance, it will provide a unique opportunity to study the impact of realistic compound events (Zscheischler et al. 2018), e.g., events similar to the drought-heat event of 2018, which caused massive heather die-off both in the field and in the ecotrons, forced by conditions like they happened in the field.

The gradient for the different ecotron units does not follow a monotonic trend for some of the key indicators (Fig. 9 and Table 2), due to the high local and inter-annual natural climate variability of the climate system. This issue could be alleviated by running the experiment for a longer period. Comparing different time frames, all extracted based on 30-year averaged GMT anomaly thresholds, shows that choosing longer time windows of 10 or 20 years leads to more clear monotonic trends (Figs. 10 and 11), which is more pronounced for air temperature-derived indices than for precipitation-derived indices. For shorter time windows of 1 to 2 years, the inter-annual and local natural variability leads to larger variations in trend for the different GMT anomaly levels. Therefore, the experiment would have to run for a long period, but the experimental time frame is constrained by the experimental setup and possible renewal. As a compromise, here, we use a 5-year experimental period. Ideally, the entire gradient should be replicated several times with different climate trajectories to average out the natural climate variability. This approach is however constrained by the high cost of the experimental setup.

Fig. 10
figure 10

Annual anomalies per GMT anomaly for increasing time window lengths (ranging from a 1-year period to a 20-year period) of the CCLM4-8-17 EC-EARTH simulation following RCP 8.5 for air temperature indices: mean air temperature anomaly (Δ T; a), annual maximum air temperature (Δ TXx; b), annual minimum air temperature (Δ TNn; c), anomaly in annual number of summer days (d), frost days (e), and the anomaly in growing season length (f). Note the different y-axis scales

Full size image
Fig. 11
figure 11

Same as Fig. 10, but now for precipitation indices: the annual accumulated precipitation anomaly (Δ PRCPTOT; a), anomaly of monthly maximum 1-day precipitation (Δ Rx1day; b), anomaly of annual number of days with more than 10 mm precipitation (Δ R10mm; c), anomaly of annual maximum length of a dry spell (Δ CDD; d) and anomaly of maximum length of a wet spell (Δ CWD; e). Note the different y-axis scales

Full size image

In the different ecotron units, we assume that the controlled variables (CO2 and CH4 concentration, air temperature, precipitation, atmospheric humidity, wind, ...) are in equilibrium with the warming level, by extracting the 5-year period in which the GMT anomaly in the driving GCM is reached. While this is a reasonable assumption, several components in the climate system will not yet be in equilibrium with the GMT anomaly at the time of simulation (e.g., glaciers, ice sheets, sea level; Zekollari 2019; Church et al. 2013. Therefore, we cannot rule out that changes in these slower components may still affect the meteorological conditions until these reach equilibrium too. For instance, a delayed melting of sea ice could alter the polar circulation and thereby affecting the mid-latitude circulation (Coumou et al. 2018), whereas ice sheet melting may affect oceanic pole-ward heat transport (Caesar et al. 2018). However, to select the time windows, we follow the same approach as the transient response to cumulative emissions (TRCE) as presented in the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (IPCC 2013). This concept describes the warming per unit of carbon emissions, which largely follows a linear relationship independent of the emission scenario (Knutti and Rogelj 2015).

The setup of the UHasselt Ecotron experiment implies that the incoming shortwave radiation will follow current weather conditions and not the weather conditions as prescribed by the RCM forcing. It is thus possible to have, for instance, clear-sky conditions and associated high incoming shortwave radiation in the field, while in the ecotron unit, a heavy precipitation event is simulated consistent with the RCM forcing. In this example, the system receives more incoming shortwave radiation than in the simulated climate. Likewise, the surface fluxes will be higher, but the resulting air temperature and moisture are corrected within the ecotron unit by the controlling devices to fully follow the boundary layer conditions as they are prescribed by the RCM.

The UHasselt Ecotron experiment allows to investigate ecosystem responses to different levels of climate change. This allows to study subtle changes in ecosystem responses such as impacts of decreased frost frequency on plant mortality (Berendse et al. 1994) and the interactions between the occurrence of mild droughts and plant acclimation for longer droughts (Backhaus et al. 2014). Although climate variables are prescribed, ecosystem-climate feedbacks originating from interactions between the biosphere and atmosphere can by partially diagnosed. For instance, heatwave reinforcements by occurring droughts (Seneviratne et al. 2010; Zscheischler and Seneviratne 2017) as well as soil moisture effects on precipitation events (Guillod et al. 2015) may be assessed by calculating imbalances in the energy budget.

Conclusions

Ecosystem experiments investigating climate change responses require a holistic, realistic climate forcing, reflecting not only the changes in the mean climate, but also representing physically consistent co-variance between climate drivers, natural variability, and changes in extreme events. To this extent, we presented a new method for creating realistic climate forcing for manipulation experiments using a single RCM simulation, and subsequently applied it on the UHasselt Ecotron experiment. To account for co-variances between variables and to fully capture the climate variability including extreme events, we selected an RCM simulation from the EURO-CORDEX ensemble based on the following criteria: (i) high skill in the local present-day climate and (ii) representative of local changes in the multi-model mean.

Based on a thorough evaluation of four key variables (air temperature, precipitation, relative humidity, and wind speed), we found that there is no single RCM-GCM combination outperforming all others for all considered variables and metrics. We made a selection of the six best-performing simulations as potential candidates and verified whether they represent the multi-model mean for the considered variables. As precipitation is considered, the most important variable in ecosystem experiments, and as most GCM downscalings have a large bias for this variable, we use the precipitation bias as the decisive factor to single out the simulation which will serve as forcing: CCLM4-8-17 driven by EC-EARTH.

The units of the UHasselt Ecotron experiment are forced with climate conditions along a global mean air temperature (GMT) anomaly gradient, representing conditions of a 0∘C (historical), + 1 ∘C (present day), + 1.5 ∘C, + 2 ∘C, + 3 ∘C, and + 4 ∘C warmer world. Five-year time windows corresponding to these warming levels are defined based on when the 30-year averaged GMT anomaly of EC-EARTH, the driving GCM, crosses these air temperature thresholds. Subsequently, the ecotron forcing is extracted from the 3-hourly RCM simulation according to the time windows.

Our new methodology provides realistic climate forcing, accounting for co-variances between climatic variables and their change in variability, well representing possible compound events. This is particularly interesting for controlled environment facilities, as their setup allows to realistically simulate future climate by controlling and measuring multiple parameters. Other controlled environment facilities could also benefit from the proposed methodology, depending on the posed research questions. The protocol for selecting a suitable regional climate simulation and extracting time series for the needed variables based on the time window defined by a global mean air temperature threshold provides a framework for different types of manipulation experiments aiming to investigate ecosystem responses to a realistic future climate change, even without a gradient approach.