Spatiotemporal variability and empirical Copula-based dependence structure of modeled and observed coupled water and energy fluxes

In alpine and prealpine regions, eco-hydrometeorological variables and processes such as soil moisture, evapotranspiration (ET), vegetation type and dynamics, and surface heat fluxes exhibit rapid changes within short distances. This is mainly due to the heterogeneity in topography, soil hydraulic properties, landuse, and climate, as well as interactions between the earth surface and the atmospheric boundary layer (Kunstmann et al. 2004, 2006; Hingerl et al. 2016; Soltani et al. 2017). The energy and water budgets in such environments are, therefore, highly controlled by the soil-type properties (Pielke et al. 1998), landcover characteristics (Dirmeyer et al. 2010), and vegetation structure (Pielke et al. 2011). The accurate spatial prediction of hydrometeorological variables can only be achieved with a distributed high-resolution hydrologic modeling approach. Such models explicitly take into account all of the domain characteristics by simultaneously solving the water and energy balance over complex mountain terrain (Bronstert et al. 2002).

In recent years, sophisticated process-based hydrological models have been developed, i.e. ALPINE3D (Lehning et al. 2006), Distributed Hydrology-Soil-Vegetation Model (Cuo et al. 2008) and JGRASS-NewAGE (Formetta et al. 2011). In these models the full system of interactions between different environments is usually taken into account; however, the equations are simplified and parameterized for the process interactions. In models like GEOtop (Rigon et al. 2006; Endrizzi et al. 2014), modeling the interactions between various hydrological, ecological and atmospheric boundary-layer processes in an interdisciplinary research area is possible, as the model covers a wide spectrum of factors in the terrestrial hydrological cycle (Endrizzi et al. 2014).

To validate the performance of hydrological models, the simulation outputs of these models are compared against observation-based runoff data and, more recently, micrometeorological measurements derived from Eddy Covariance (EC) techniques (e.g. Rigon et al. 2006; Hingerl et al. 2016). EC-based information is now considered as a valid source for a model's calibration and validation (e.g. Rosero et al. 2010; Decker et al. 2012). Traditionally, linear statistical measures (e.g. correlation coefficient r) are used for model performance evaluation. However, using simple linear r-values between simulated and observed hydrometeorological datasets, which typically exhibit nonlinear characteristics, may not be an appropriate way to determine these complex relationships (Bárdossy & Pegram 2009). Copula functions can appropriately obtain underlying dependence structures of hydrometeorological variables (Dupuis 2007), including their complexities in time and space (Laux et al. 2011).

Copula-based models have been used in a variety of experimental studies for different purposes worldwide. In the field of hydrology, Bárdossy (2006) described the spatial variability of the groundwater quality parameters using bivariate empirical Copula (as an alternative to variograms and covariance functions). Using Copula functions, Sugimoto et al. (2016) made an attempt to detect the temporal changes of catchments independent from climate change by investigating the long-term discharge records. Li et al. (2016) estimated the bivariate flood quantiles by combinations of peak discharge and flood volume using Copulas in China. In a similar study, with regard to urban catchment applications, flood frequency curves were derived using bivariate rainfall distribution based on Copula functions by Balistrocchi & Bacchi (2017). The following studies are examples for applying Copulas in the hydrometeorological field: modelling the daily precipitation features in West Africa (Laux et al. 2009), spatial recorrelation of regional climate model (RCM) precipitation to generate unbiased climate change scenarios over the Rhine basin (Bárdossy & Pegram 2012), spatiotemporal patterns of precipitation extremes in China (Zhang et al. 2013), and bias correction of dynamically downscaled precipitation fields in Germany (Mao et al. 2015). However, to our knowledge, no study has used Copula-based models so far for evaluating the performance of hydrological simulations, as described and presented in this study.

Our investigation is performed for two heterogeneous catchments (ranging from small to mesoscale) within the TERrestrial ENvironmental Observatories (TERENO) prealpine region located in southern Germany using GEOtop2.0 to jointly simulate the water and energy budgets over two summer episodes in 2013 and 2015.

Previous studies in the region have focused on how climate change impacts runoff generation, surface and sub-surface water balances, biosphere-atmosphere exchange (greenhouse gases), and energy balance closure (EBC) parameterization (Kunstmann et al. 2004; Ott et al. 2013; Unteregelsbacher et al. 2013; Eder et al. 2014; Wang et al. 2014; Wolf et al. 2016; Soltani et al. 2017; Zeeman et al. 2017). Inverse distributed hydrological modelling has been studied by Kunstmann et al. (2006) for the alpine/prealpine Ammer River catchment through coupling the Parameter ESTimation tool (PEST) in the WaSiM hydrological model. Hingerl et al. (2016) modeled the spatiotemporal variability of the water and energy flux components using GEOtop1.45 for this prealpine catchment in the TERENO region. In both cases, however, the hydrological models were calibrated against the runoff measurements only, and no soil moisture profile variation and/or radiative-turbulent flux variabilities have been accounted for in the simulation. Moreover, previous studies only attempted to evaluate the models' performances using the traditional linear statistical metrics (e.g. R², RMSE).

Given the above described gaps, in this study we address the following objectives: (i) quantifying the spatiotemporal variability of the hydrometeorological variables of the turbulent fluxes as well as the surface temperature and evapotransipration (ET) with respect to the elevation-gradient effect using high-resolution EC-based measurements and spatial hydrological simulations; (ii) simulating the coupled water and energy balances at a very high spatial resolution using the physically based hydrological model GEOtop; and (iii) estimating the underlying dependence structures of the observed and modeled water and energy fluxes using the nonlinear-based approach bivariate empirical Copula in the TERENO prealpine region.

The Rott catchment (∼55 km²) and the Upper-Ammer catchment (∼300 km²) are located in the TERENO prealpine region in southern Germany (Figure 1). The highest and lowest points in the Rott catchment range from 902 to 543 m, where the Raisting discharge gauge is located at the outlet of the basin. The corresponding values for the Upper-Ammer catchment range from 2,129 to 619 m, where the Peissenberg discharge gauge is located at the outlet. The micrometeorological measurements from the Fendt EC site (47.831°N, 11.061°E) and the Rottenbuch EC site (47.730°N, 10.061°E) are used for model validation.

The Fendt site is placed in an area characterized as pasture landuse (Figure 2(a)) with a Histosol soil type (Figure 2(b)). The Rottenbuch site is similarly located in pasture (Figure 2(a)), but has a soil profile with sandy loam soil (Figure 2(b)) at the surface, and a gravel layer at ∼50 cm depth (Soltani et al. 2017). Additionally, the southern portion of each catchment, especially the Upper-Ammer, has the highest elevations and slopes greater than 50° in the southern TERENO region (Figure 2(c)).

The meteorological dataset was obtained from three major sources for the calibration and validation periods from 1 May to 31 July 2013 and 2015. First, the EC-based turbulent fluxes and the micrometeorological measurements were obtained from the TERENO EC sites (Mauder et al. 2013). Second, the hourly forcing data for the Wielenbach and Hohenpeissenberg, Kohlgrub, and Oberammergau stations were obtained from the Deutscher Wetterdienst (DWD) (https://werdis.dwd.de). Third, the hourly data of precipitation for the Diessen rain gauge and runoff for the Raisting and Peissenberg discharge gauges were compiled from Bayerisches Landesamt für Umwelt (http://www.hnd.bayern.de).

The experimental period was selected for the simulation due to the fact that there usually remain uncertainties in the measured wintertime water and energy fluxes over the region (e.g. Hingerl et al. 2016). Additionally, the peak runoff volume and the highest EBCEBC were observed during summer periods (e.g. Soltani et al. 2017).

The TERENO prealpine region in Southern Germany is a network of observatories to investigate the impact of global environmental change on the terrestrial systems (Zacharias et al. 2011). It is a sophisticated set of EC-flux stations, boundary layer sounding systems, and hydrometeorological networks of precipitation gauges, soil moisture and soil temperature sensors (Wolf et al. 2016). This enables the investigation of various aspects of atmospheric–land-surface–subsurface interactions on water and energy fluxes, and also greenhouse gases.

GEOtop2.0 (Endrizzi et al. 2014) is a fully distributed physical-based model of the mass and energy balances, and is designed for simulations in mesoscale catchments over complex terrains (Bertoldi et al. 2004; Rigon et al. 2006; Endrizzi et al. 2014). The model deals with the interaction effects of topography on the energy balance and hydrological cycle. It describes the three-dimensional water flow in the soil and the energy exchange with the atmosphere, considering the radiative/turbulent fluxes. GEOtop2.0 (hereafter GEOtop) uses sophisticated integration methods, which, in turn, allow convergence even in cases where parameters have nearly discontinuous behavior (Endrizzi et al. 2014). The core components of GEOtop are fully described in Endrizzi et al. (2014).

The following static datasets were used for simulation by the GEOtop model. The DEM (digital elevation model) 90 × 90 m² (100 × 100 m²) was obtained from the Shuttle Radar Topographic Mission (STRM), http://srtm.csi.cgiar.org, (Kunstmann et al. 2006) for the Rott (Upper-Ammer) catchment. The river network, as well as terrain aspect, slope, and sky view were calculated using ArcMap GIS based on techniques described by Ghesla & Rigon (2006). The land cover, with a resolution of 250 × 250 m² (Rott) and 150 × 150 m² (Upper-Ammer), as well as the soil type with a resolution of 2 × 2 km² for both catchments were taken from Hingerl et al. (2016) and Kunstmann et al. (2006), respectively. Then, they were separately interpolated to the same resolution as the DEMs following Bertoldi et al. (2004). The soil is discretized in 13 layers, with thicknesses increasing from the surface to the deep layers. The top eight layers starting from the surface have thicknesses ranging from 0.1 to 0.5 m, with respect to the vertical gradients of water pressure and temperature (Endrizzi et al. 2014), while the lowest five layers have thicknesses ranging from 1.0 to 5.0 m.

The hourly meteorological data of precipitation, temperature, wind speed, wind direction, and radiation components were provided as input forcings for the model simulations. In the GEOtop code, air temperature and precipitation are spatially distributed according to Liston & Elder (2006). The integration of the model run is one-hour interval. Detailed information about the interpolation of the meteorological data in GEOtop is described in Endrizzi et al. (2014).

The GEOtop model is first applied to a summer episode starting from 1 May to 31 July 2013. A two-week spin-up period starting from 15 April to 30 April 2013 was conducted, as we found that this is sufficient for our experimental period. The simulation is then used for determining the parameters' sensitivities and model calibration. In our simulation experiments, the key model parameters acting on the water and energy budgets (Table 1) are iteratively estimated based on the most sensitive parameters identified by the previous hydrological models run in the study area (Kunstmann et al. 2006; Hingerl et al. 2016).

We used a trial-and-error procedure for the model calibration within an accepted range of values (Table 1). The effect of selected parameters on the model outputs were individually tested and optimized until the best fit was obtained between the simulated and measured data. The calibration results are given in Table 1. We found that changes in the parameter values from the northern (Fendt) to the middle parts (Rottenbuch) of our TERENO prealpine region are realistic and explicable based on changes in e.g. climatic-environmental conditions and land-surface properties for those regions. For example, the difference in the calibrated values of SoilAlbVisWet can be explained by differences in the soil properties for those two areas. In Fendt, the dominant soil-type is Histosol, whereas Rottenbuch is covered by a sandy-loam soil-texture.

The soil parameter VerticalHydrConductivity determines the highest capacity of infiltration and highly affects both the surface runoff generation and the magnitude of peak discharge. The parameter of HorizontalHydrConductivity impacts the flood tail as it determines the subsurface runoff quantity (Rigon et al. 2006). These parameters were also reasonably calibrated (see Table 1). Further, the initial condition on water table depth was set as: InitWaterTableHeightOverTopoSurface/InitWaterTableDepth = −5,000/1,000 mm (Rott catchment) and = −10,000/1,500 mm (Upper-Ammer catchment). The initial condition for the soil pressure, i.e. InitSoilPressure and InitSoilPressureBedrock, were kept constant by default for both catchments.

Model evaluation statistics are commonly utilized for comparing model outputs against measurements. There is a large variety of evaluation metrics, but no single metric encapsulating all aspects of interest exists. For this reason, the following statistical metrics were used to evaluate the performance of GEOtop in this study. (1) Mean Bias (MB), which provides a good indication of the mean over- or underestimate of simulations. (2) Root Mean Square Error (RMSE), which provides a good overall measure of how close modelled values are to observed values. (3) Coefficient of Determination (R²), which measures the strength of the linear relationship between two variables. (4) Nash-Sutcliffe Efficiency (NSE), which is a normalized statistic that determines the relative magnitude of the residual variance compared to the measured data variance (Nash & Sutcliffe 1970). The NSE ranges between −∞ and 1, with NSE = 1 being the optimal value. Values between 0 and 1 are generally viewed as acceptable levels of performance, whereas values ≤0 indicate that the mean observed value is a better predictor than the simulated value (Moriasi et al. 2007). (5) Coefficient of Efficiency (COE), based on Legates & McCabe (2012) and Legates & McCabe (1999), is used for measuring model performance. A perfect model has a COE = 1. Although COE has no lower bound, a value of COE = 0 has a fundamental meaning. It implies that the model is no more able to predict the observed values than the observed mean. Therefore, since the model can explain no more of the variation in the observed values than can the observed mean, such a model can have no predictive advantage. For negative values of COE, the model is less effective than the observed mean in predicting the variation in the observations (Carslaw 2015). (6) Index of Agreement (IOA) is used in model evaluation (Willmott et al. 2011). It ranges between −1 and +1 with values approaching +1 representing better model performance. An IOA of 0.5, for example, indicates that the sum of the error-magnitudes is one half of the sum of the observed-deviation magnitudes. When IOA = 0, it signifies that the sum of the magnitudes of the errors and the sum of the observed-deviation magnitudes are equivalent. When IOA = −0.5, it indicates that the sum of the error-magnitudes is twice the sum of the prefect model-deviation and observed-deviation magnitudes. Values of IOA near −1 indicate that the model-estimated deviations are poor estimates of the observed deviations (Carslaw 2015).

In addition, according to the fact that the relationships between hydrometeorological variables may be highly nonlinear, we employ the concept of the empirical Copula (described in the following section). Unlike the linear statistical measures, the nonparametric empirical Copula-based approach is able to describe the underlying joint behavior of water and energy fluxes (Genest & Favre 2007; Serinaldi 2009; Laux et al. 2011).

The calibrated GEOtop model simulated an independent period of 1 May to 31 July 2015 for validation.

Hydrometeorological variables usually show nonlinear behaviors and hence their relationships are very complex. Traditionally, the pairwise dependence between these variables has been described using classical families of bivariate distributions such as normal, log-normal, gamma, and extreme-value distributions. The main limitation of this approach is that the individual behaviors of the two variables or transformations thereof must then be characterized by the same parametric family of univariate distributions. Copula models, however, avoid this restriction (Genest & Favre 2007).

Figure 3 shows the temporal variations of surface turbulent fluxes on daily resolution as a calendar plot, where the wind angle is scaled to the wind speed to highlight both the direction and the strength of the wind of a particular day for the Fendt and Rottenbuch EC sites. At Fendt, the highest turbulent flux with a magnitude of more than 250 W m⁻² on 7th July was associated with northeasterly winds during the experimental period. The corresponding value for the Rottenbuch (>140 W m⁻²) was connected to northwesterly flows that occurred on 15th July. Our results indicate a clear increase in the turbulent flux concentration from May to July at both sites due to the increase in solar radiation; however, the temporal variation patterns are different at the individual sites.

The spatial variability of heat and water fluxes in the TERENO prealpine region is significantly affected by diversity in topography, radiation and wind components, soil moisture properties as well as land cover and vegetation types. The surface temperature and ET values (Figure 4) range from 9 to 21 °C and 17 to 160 mm, respectively, in the Upper-Ammer catchment. The corresponding values for the Rott catchment range from 15 to 20 °C and 42 to 153 mm, respectively. The maximum values are simulated at parts of the basin with: (i) the lowest elevation, slope (controls the net radiation e.g. Dubayah et al. 1990), and wind speed (Chow et al. 2006); (ii) Histosol or sandy-loam texture soil; and (iii) in-pasture or peatland vegetation cover (determines storage capacity (e.g. Bertoldi et al. 2010)). The different forest types tend to have the lowest ET values across the region. The mean surface temperature reflected the diverse surface topography of the TERENO prealpine region by exhibiting an elevation-dependent temperature variation of 10 °C in the Upper-Ammer area and 4 °C in Rott. The spatial patterns obtained herein should be interpreted with care, as a soil-type map of 2 × 2 km² used for model input is rather coarse with respect to other models' input layers and computational resolution.

The comparisons of measured and simulated discharge in the TERENO prealpine catchments are shown in Figure 5. Overall, GEOtop was capable of replicating appropriately the river discharge dynamics with a high efficiency (see Table 2). The NSE indicates a high performance of the model to represent the temporal variability of discharge (NSE = 0.86 for Rott and NSE = 0.81 for Upper-Ammer) in which the volumes of discharges are reasonably reproduced at both catchments. That is, the total simulated discharge volume for the 3-month period was 3,886 m³ in Rott and 53,884 m³ in Upper-Ammer. This indicates an underestimation of about 25% and 15%, respectively, compared to the measured discharge volumes of 5,168 m³ and 62,509 m³, respectively. This was also represented by MB values of −0.34 m³ for Rott and −2.01 m³ for Upper-Ammer. In both catchments, an increased runoff volume in early June peak flow highlights the importance of snow dynamics for runoff generation in the region.

The model captures the peak flow well in Rott, but underestimates it in the Upper-Ammer catchment. This might be explained by the lack of meteorological stations, which can result in considerable errors in the spatial interpolation by the model. These differences may further be explained by the rapid climate zone changes in a small spatial area or by the snow dynamics effect on the behavior of surface runoff during the springtime (Kunstmann et al. 2006). The values of RMSE are very low (i.e. Rott: 1.29; Upper-Ammer: 7.36) while R² are high (i.e. Rott: 0.80; Upper Ammer: 0.88), indicating that simulated and measured discharges have low residuals and strong linear relationships. According to Figure 5, the GEOtop model, similar to the calibration period, indicates a good performance to replicate the river runoff by capturing the peak flows at both catchments during the validation period (see Table 2 for detailed information about the statistical measures).

Figure 6 displays the time series of the site-scale, hourly measured and simulated soil moisture for the Fendt and Rottenbuch EC sites at three different depths from May to July for the calibration and validation periods. Over the calibration period, GEOtop indicates a strong linear relationship between the simulated and measured −5 cm soil moistures at Fendt (R² = 0.94) and Rottenbuch (R² = 0.84).

The daily fluctuation of observed soil moisture at Fendt ranges from about 25% to 75%, while it ranges approximately from 20% to 50% for the Rottenbuch site. The different soil texture and soil type could explain this discrepancy, where sandy-loam soil texture and Histosol type prevail. This is why the maximum soil moisture close to the surface layer in Rottenbuch was about 50%. Apart from a different soil type, which results in a high variation in the near surface soil water content in Fendt, the soil moisture also depends strongly on the variation of the groundwater table depth for that area (Wolf et al. 2016). When there is a deficit in rainfall, the rate of soil moisture drops quickly and reaches up to 20% at both sites. Due to scale differences between grid cells and point measurements it is not easy to compare precisely the simulated soil moisture with those of the measurements (Manabe et al. 2004).

Nevertheless, we found that GEOtop is appropriately capable of reproducing infiltration and daily cycle of soil moisture evaporation associated with the rainfall events with the lowest mean bias and error and the highest efficiency and agreements. The model shows a similar performance to replicate the soil moisture during the validation period for both EC sites (see Table 2 for further details).

The distributed soil moisture maps at −5 cm layer simulated by the model, as shown in Figure 7, are found to be highly influenced by the heterogeneous topography (terrain slope) and different landuse as well as soil texture patterns. With regard to topography, obviously the most (least) saturated soil water content greater than 70% (less than 30%) is mainly found in the low lands (steep gradients), which are mostly located in northern (southern) parts of the TERENO prealpine region. The spatial distribution of soil moisture strongly follows the spatial distribution of the ET pattern (see Figure 4), which was also confirmed by the site-scale results at both Fendt (70%) and Rottenbuch (45%).

Furthermore, in Rott, 87% (89%) of precipitation leaves the basin as ET and only 13% (11%) is consumed for runoff generation during the calibration (validation) period. In Upper-Ammer, however, the water balance condition is quite different. Approximately 68% (56%) of the precipitation leaves the catchment as discharge and ET consumes around 32% (44%) during the calibration (validation) period, which indicates the importance of runoff generation by increasing flooding potential in the lowlands during the summertime.

The comparison between simulated and measured site-scale diurnal surface energy fluxes of Rn (net radiation), LE (latent heat), H (sensible heat) and G (soil heat) in the TERENO prealpine EC sites is illustrated in Figure 8. Overall, the EC-based diurnal cycles of energy fluxes were well reproduced by GEOtop, which is confirmed by high values of correlation coefficients, efficiency and agreements between the measured and simulated Rn, H, LE and G (see Table 3). Some bias is observed in Rn (MB = −17.4 W m⁻² at Fendt and MB = −20.7 W m⁻² at Rottenbuch), which indicates that Rn is underestimated. To identify the possible reasons, the modeled radiation components (i.e. incoming short/longwave and outgoing short/longwave) were compared against the measurements. We found that the modeled Incoming Shortwave Radiation (ISR) shows considerable biases (MB = −115.6 W m⁻² at Fendt and MB = −120.9 W m⁻² at Rottenbuch) and was underestimated by the model.

The model slightly overestimates LE (MB = 14.3 at Fendt and MB = 24.8 at Rottenbuch), in particular during the first part of the day with a peak at around 11:00. The midday overestimation of the simulated LE flux at the EC sites might be explained by the lack of EBC in the EC-based measurements, where the imbalance (residual energy) at Fendt is 31% (Figure 9(c)).

At Rottenbuch, however, not only it is even higher (39%) as illustrated in Figure 9(d), but also LE flux partitioning to Rn is at least twice as low as that of in Fendt. A lower LE partitioning at Rottenbuch can be explained by a different soil texture for that area, where the sandy-loam at the surface and a gravel layer at about 50 cm depth is predominant, resulting in the precipitation infiltrating quickly into the deep ground, thus the surface dries out rapidly soon after the rainfall events (Soltani et al. 2017). The G flux is overestimated during early morning hours with a lower mean bias similar to Rn values. The main consumer of Rn is LE over the experimental period at the study sites.

The intercomparison of the simulated and measured energy balance components in the TERENO EC sites are shown in Figure 9(a) and 9(b). GEOtop overstimates LE flux at Fendt (with a slope of 1.04) and even more at Rottenbuch (1.07). This could be due to the fact that the EC-baesd technique usually underestimates turbulent fluxes, in particular LE measurement, as reported worldwide (e.g. Hendricks Franssen et al. 2010; Stoy et al. 2013; Imukova et al. 2016), also in the TERENO prealpine region (Eder et al. 2014; Soltani et al. 2017).

As shown in Figure 10(c) and 10(d), the EBC is estimated to be 0.69 and 0.61 with R² values of 0.92 and 0.88 at Fendt and Rottenbuch sites, respectively. In terms of energy balance ratio (EBR), a higher EBR value of 0.70 is calculated at Fendt indicating that the minimum heat and water vapor fluxes are lost for that area compared with the Rottenbuch site, where a lower slope of 0.61 and EBR value of 0.62 are found. Soltani et al. (2017) suggested that this imbalance at Rottenbuch is likely due to the advection of heat and vapor in that area (Figure 9(d)). Similar results were also found by Hingerl et al. (2016). Other studies have also shown an energy imbalance in the region (e.g. Kunstmann et al. 2013; Eder et al. 2014).

The spatial distributions of H and LE fluxes modeled by GEOtop are strongly determined by the heterogeneity of the TERENO prealpine region (Figure 10). In terms of topography and terrain slope, for example, the maximum H and minimum LE are found at the highest elevations at both catchments, which are in good agreement with the ET and surface temperature distributions (see Figure 4).

The energy partitioning of Rn, LE, H and G fluxes indicates that in Rott, approximately 78.5% (82.6%) of Rn leaves the catchment as LE, 17% (12.6%) leaves as H, and 4.5% (4.8%) enters the soil as G during the calibration (validation) period. The energy balance condition in the Upper-Ammer catchment, however, shows a different pattern, where the portion of H is significantly increased. Approximately 65% (69.5%) of Rn leaves the catchment as LE, 31% (27%) leaves as H and 4% (3.5%) enters the soil as G during the calibration (validation) period.

The empirical Copula densities of the measured and simulated variables of discharge and soil moisture are illustrated in Figure 11(a). For Rott, the empirical Copula density between the measured discharge and soil moisture indicates a strong symmetric dependence structure, where the highest density is found for the lower left, and a second maximum in the upper right. This means that the measured discharge and soil moisture show the highest correlations for both the low and extreme values, respectively. The density function for the simulated dependence structure shows a similar pattern to that for the measured values, but is different in terms of the lower and upper densities. This implies that the modeled discharge and soil moisture values represent lower dependencies (correlations) than those observed. In addition, a positive bias for the simulated discharge is observed at the low, but also at very high values. The Copula densities of both measured and simulated discharges and soil moistures in Upper-Ammer indicate a significant symmetric dependence structure with the highest densities in the lower-left for the measured, and in the upper-right for the simulated. In other words, the highest correlation between the measured discharge and soil moisture is found at the very low values, whereas the corresponding correlation for the simulations is observed at the high values. The empirical density between discharge and rainfall intensity was analyzed, as well. Here, no clear pattern of the dependence structure was found (not shown).

As shown in Figure 11(b), the distribution of measured against simulated discharge is asymmetrical, with the highest density in the upper-right corner. It indicates that measured and simulated discharges are strongly concordant in the higher ranks of the distribution. The concordance is weaker in the lower ranks in both catchments. This conveys that the calibrated model is more capable of replicating the high (peak) streamflow values than those of very low values at both Rott and Upper-Ammer catchments. However, GEOtop captured the low flow values for Upper-Ammer better than it did for the Rott (see Figure 5). This can be seen in Figure 11(b), where a second density maximum is found at the lower-left corner, indicating a reasonable performance of the model to simulate the low discharge at the Upper-Ammer catchment. However, the dependence structure between the measured versus simulated soil moisture is different. This means that the modeled low values show the highest agreements with those observed at both Fendt and Rottenbuch EC sites. Also, the model shows a good performance to replicate the higher values (especially in Fendt), as a second density maximum can be found at the upper-right corners at both sites. See Figure 6 for clarification, where the soil-moisture time series are plotted.

Figure 12(a) shows the empirical Copula densities of the measured and simulated turbulent fluxes. A nonlinear structure of the joint relationship for both measured and simulated latent and sensible heat flux is observed, with strongly varying densities in the different percentiles. The highest density is found for the high turbulent flux values in the upper-right corner, whereas the lower tails do not show a significant dependence structure in both Fendt and Rottenbuch. Hence, the modeled latent and sensible heat fluxes represent the highest correlations at the high values, whereas those observed show a lower correlation with the corresponding values.

The empirical Copula densities for the measured against simulated latent and sensible heat fluxes, according to Figure 12(b), are almost asymmetrical for both sites. This implies that the highest density functions between the measured and simulated latent (sensible) heat flux is seen in the upper-right (lower-left) corner, where the extreme (very low) values are found. In other words, GEOtop shows a better performance to replicate the high (low) values of LE (H) at both Fendt and Rottenbuch EC sites. The empirical Copula density functions for further hydrometeorological variables were estimated, as well. We could not identify clear dependence structures between the discharge – soil moisture and sensible heat flux – latent heat flux (not shown).

The spatiotemporal variability of water and energy fluxes was analyzed using both the GEOtop model and a Copula-based approach for two catchments and two recent episodes in the TERENO prealpine region.

Using GEOtop, the water and energy fluxes could be replicated with high performance and low bias. The model reasonably captured the peak discharge observed in early June for the Rott catchment but underestimated discharge in the Upper-Ammer catchment. Simulated streamflow was characterized by high values of R² and low residuals (RMSE) when compared with observations. The daily cycle of multiple-layer soil moisture variations was appropriately described by the model, which again indicates low biases and high efficiencies. The EC-based diurnal cycles of energy fluxes were well reproduced by GEOtop; however, the model slightly overestimated LE, especially during the early morning due to the lack of EBC in the EC-based measurements (imbalance of 31% at Fendt and 39% at Rottenbuch). The spatial distributions of water and energy fluxes revealed that in the Upper-Ammer catchment around 70% of precipitation leaves the catchment as discharge, compared to ∼10% in Rott. Around 30% of Rn leaves the catchment as H, while only ∼15% in Rott. The model results obtained for the validation period were satisfying, indicating that the estimated parameters are reasonably calibrated.

The linear statistical measures applied are assumed to not be capable of representing the interaction between the hydrometeorological variables, and we, therefore, employed a bivariate empirical Copula-based dependence structure analysis. We first found that the bivariate dependence structure patterns of both measured and simulated hydrometeorological variables considered in this study are very similar, representing a reasonable calibration of the GEOtop model. These non-linear features in the dependence structure of measured and simulated individual hydrometeorological variables are observed with the highest densities (or best fit between the modeled and observed values) either in the lower or upper ranks, i.e. in the low or high values, exhibit a worse model calibration for the middle ranks of the data. Finally, the Copula-based model performance analysis applied can be considered for model evaluation in the hydrological model community in addition to traditional model performance analyses.

This research was financially supported by the Helmholtz Research School on Mechanisms and Interactions of Climate Change in Mountain Regions (MICMoR). Access to TERENO and TERENO-ICOS infrastructures is gratefully acknowledged. The support by the landowners of the TERENO sites and technical staff of KIT/IMK-IFU is also appreciated. The contribution of Matthias Mauder was conducted within the Helmholtz Young Investigator Group ‘Capturing all relevant scales of biosphere-atmosphere exchange – the enigmatic energy balance closure problem’, which is funded by the Helmholtz-Association through the President's Initiative and Networking Fund. We thank Kayla Stan (Dept. of Earth and Atmospheric Sciences, University of Alberta) for proofreading the manuscript.