Projection of Droughts as Multivariate Phenomenon in the Rhine River

: Drought is a complex phenomenon whose characterization is best achieved from a multivariate perspective. It is well known that it can generate adverse consequences in society. In this regard, drought duration, severity, and their interrelationship play a critical role. In a climate change scenario, drought characterization and the assessment of the changes in its pattern are essential for a proper quantiﬁcation of water availability and managing strategies. The purpose of this study is to characterize hydrological droughts in the Rhine River in a multivariate perspective for the historical period and estimate the expected multivariate drought patterns for the next decades. Further, a comparison of bivariate drought patterns between historical and future projections is performed for di ﬀ erent return periods. This will, ﬁrst, indicate if changes can be expected and, second, what the magnitudes of these possible changes could be. Finally, the underlying uncertainty due to climate projections is estimated. Four Representative Concentration Pathways (RCP) are used along with ﬁve General Circulation Models (GCM). The HBV hydrological model is used to simulate discharge in both periods. Characterization of droughts is accomplished by the Standardized Runo ﬀ Index and the interdependence between drought severity and duration is modelled by a two-dimensional copula. Projections from di ﬀ erent climate models show important di ﬀ erences in the estimation of the number of drought events for di ﬀ erent return periods. This study reveals that duration and severity present a clear interrelationship, suggesting strongly the appropriateness of a bivariate model. Further, projections show that the bivariate interdependencies between drought duration and severity show clearly di ﬀ erences depending on GCMs and RCPs. Apart from the inﬂuence of GCMs and RCMs, it is found that return periods also play an important role in these relationships and uncertainties. Finally, important changes in the bivariate drought patterns between the historical period and future projections are estimated constituting important information for water management purposes.


Introduction
Droughts belong to the ranks of increasingly important phenomena, affecting water and food supply. A detailed analysis of droughts can help mitigate their effects and manage the water resources in a more efficient way. The knowledge and assessment of the repercussions of changes in the pattern Water 2020, 12, 2288 2 of 14 of drought events is a key issue in present and future decisions. In this context, several researchers have assessed a change in drought characteristics which are attributable to global warming [1][2][3].
However, a clear definition, a proper understanding, forecasting, and managing of droughts is a complex task due to the variety of temporal and spatial scales at which droughts occur, together with their diverse direct and indirect causes and consequences [4]. Different classifications of droughts have been proposed and are defined in different contexts [5][6][7]. For example, droughts can be expressed in terms of meteorological, agricultural, hydrological, or regional aspects, and for their definition, a number of indices have been introduced. Some of these indices depend on a single parameter such as the Standardized Precipitation Index (SPI) calculated from precipitation time series while others are defined in a multivariate perspective. The use and selection of an appropriate index to analyze droughts depends on the purpose of the research [8]. However, there seems to be a scientific consensus that there is no best drought index and that a quest for the best index is useless [9]. A comparison of different drought indices can be seen in [10], in which variations of droughts over several regions in the world during the last millennium are examined. Details can be found in [11][12][13][14]. In a climate change scenario, several studies addressing changes in droughts using the SPI have been performed. For instance, Ref. [15] analyzed future changes in meteorological droughts in the Lower Mekong River Basin concluding that the Mekong Delta is expected to experience a significant increase in drought events.
Hydrological drought refers to a lack of water in the hydrological system, manifesting itself in abnormally low levels in lakes, reservoirs, and groundwater [9]. The Standardized Runoff Index (SRI) has been used in various studies to characterize hydrological droughts. It can be understood as a standardized difference or anomaly in regards to a normal situation. In a study, Ref. [16] showed a comparison of SRI behavior with that given by the SPI. They conclude that on time scales ranging from monthly to seasonal, the SRI is a useful complement to the SPI for depicting hydrological aspects of droughts. A study comprising seven large river basins around the world analyzed the propagation of forcing and model uncertainties on hydrological drought characteristics [3] using the runoff index (RI). The drought characteristics' total drought magnitude and duration were considered independently. This means, droughts were addressed in a univariate realm. Analyzing the statistical properties of drought magnitude, the results showed supporting evidence of the influence of the Representative Concentration Pathways (RCP) on this characteristic. The authors also concluded that, on average, more severe droughts can be expected in the study regions under the RCP8.5 by the end of the 21st century. Furthermore, from all the basins investigated, the Rhine at Lobith exhibits the strongest increase in drought magnitude and duration under this scenario.
It is generally recognized that the discharge regime of the Rhine River will be affected by climate change [17]. Some studies show that important changes in discharge time series are projected due to climate change, including the mean, low, and high flows [18] for the mid-and end of the 21st century. Furthermore, mean and high flows are estimated to vary importantly between seasons (winter, summer). Others reported a change in the flow generation mechanism in the Rhine River under climate change [19]. Added to this, an increase in the frequencies of both low and high flows are projected.
In general, the analysis of droughts in the context of climate change are addressed by characterizing them by a single attribute. In other words, droughts are seen as a univariate phenomenon. An estimation of the possible changes of drought events due to climate change from a multivariate perspective has not yet been addressed. The aim of the present study is to fill this gap and to model hydrological drought events as a multivariate phenomenon in the context of climate change applied in the Rhine River at Lobith. To achieve this, we characterize droughts as a multivariate phenomenon in both a historical and projection period and analyze variations in their multivariate pattern, i.e., analyze the expected changes in the joint duration-severity behavior for different recurrence intervals (frequencies). Runoff simulation is performed by the application of the rainfall-runoff model HBV to simulate daily streamflow. Bias corrected time series from five General Circulation Models (GCM) and four RCPs from the latest IPCC assessment report AR5 are used (a detailed description is found in [20]. The (multivariate) drought patterns are then characterized for both the historical as well as the projection period and then compared. The number of projected drought events is calculated for each driving model. This characterization of drought patterns in both historical and projection periods will give an inside into the effect of climate change on the development of drought seen as a multivariate phenomenon

Study Area and Input Data
The study area corresponds to the Rhine basin at gauging station Lobith with an area of 160,000 km 2 . It is a warm temperate, humid area with warm summers with an average annual temperature of T = 8.7 • C and annual precipitation of P = 1038 mm. The dominant land covers in the area consist predominantly of forest (25%), cropland (38%), and grassland (9%). The altitude varies between 15 and 4075 m a.s.l. with a mean value of 497 m a.s.l. The average discharge over the area in the 1971-2000 period was 457 mm with an average runoff coefficient of 0.44. More information and description of the study area can be found in the introductory paper of the Inter-Sectoral Impact Model Intercomparison Project Phase 2 ISI-MIP2 [21]. Projected discharge time series were obtained by running the calibrated rainfall-runoff model HBV. Calibration and validation were carried out using the WATCH forcing data set and followed the ISI-MIP2 protocol (www.isimip.org). In this protocol, detailed information can be found as specification of periods for calibration and validation, spin-up, and projection periods, among others. WATCH data are based on the 40-year ECMWF Re-Analysis (ERA-40) and reordered reanalysis data [22,23]. The rainfall-runoff model is forced using the bias-corrected outputs [20] of the five general circulation models: GFDL-ESM2M, HadGEM2-ES, IPSL-CM5A-LR, MIROC-ESM-CHEM, and Nor-ESM1-M [24]. The uncertainty ranges for precipitation and temperature associated with these models are comparable with all models of the Coupled Model Intercomparison Project Phase 5, CMIP5 (protocol-report www.isimip.org).
Four future emission scenarios described by RCPs were considered in this study, namely, RCP2.6, RCP4.5, RCP6.0, and RCP8.5. Referring to the two extreme climate conditions, RCP2.6 represents a scenario in which greenhouse gas emissions peak in the 2010-2020 period, whereas the latter (RCP8.5) assumes that emissions continue rising in the present century [25].
Information is divided into three periods: historical , mid-century (2036-2065), and end-century (2070-2099) in order to model and characterize drought characteristics and their interrelationship in different periods of time.

Hydrological Model
The conceptually based semi-distributed HBV model was used in this study. The HBV model concept [26] was developed by the Swedish Meteorological and Hydrological Institute (SMHI) in the early 1970s and modified at the Institute of Hydraulic Engineering, University of Stuttgart, Germany. The HBV model comprises routines for the calculation of snow accumulation and melts, soil moisture, as well as runoff generation, runoff concentration, and flow routing in the river network. Soil moisture was calculated by balancing precipitation and evapotranspiration using field capacity and permanent wilting point. Runoff generation was simulated by a nonlinear function of the actual soil moisture and precipitation. The runoff concentration was modelled by two parallel nonlinear reservoirs representing the direct discharge and the groundwater response. The Muskingum method was used for flood routing between the river network nodes. The physical meaning of model parameters and their range is given in Supplementary Table S1. Additional information about the HBV model in general and the University Stuttgart subtype can be found in the specialized literature [27,28]. The HBV hydrological model was setup for the Rhine catchment at Lobith. The model was calibrated using the robust parameter estimation (ROPE) algorithm [27] where the Nash-Sutcliff efficiency (NSE) [29] coefficient was used as objective function. The calibration and validation of the model was performed for the periods 1991-1998 and 1999-2006, respectively. The modelling time steps were daily and the evaluation of objective function was performed on a monthly time step for the calibration and validation period.

Estimation of Drought Characteristics
Drought is characterized by two variables, namely, the severity of an event and the associated duration of that event, both derived from the calculated SRI set. The index (SRI) is calculated through the quantile function of the standardized normal distribution applied on the fitted univariate cumulative distribution function (CDF) of the data.
Drought duration is defined as the aggregated time in which the SRI is negative and preceded and followed by positive SRI values [30]. Drought severity is estimated as the sum of the SRI values during a drought event. Mathematically: S refers to severity and d denotes the duration of the drought event. As SRI takes only negative values in a drought event, S is a positive real number. Different time scales can be defined in the calculation of the SRI [30]. In this study a 3-month time scale was adopted. With this, we particularly intended to focus on seasonal components of droughts. Larger time scales (for example, 1-year) likely do not reflect some droughts occurrences. In agreement with this, [16] concluded that the use of monthly to seasonal scales for the calculation of SRI is useful for depicting hydrologic aspects of droughts. Descriptive statistics of duration and severity for the historical period are summarized in Supplementary Table S2. Minimum, maximum, and average duration values as well as minimum, maximum, and average severity is presented.

Bivariate and Probabilistic Model
In this section, only the basic concepts of the bivariate copula models are introduced. More details can be found in the given literature below. In general, copula functions represent a method for modelling the interrelation between different random variables (a comprehensive summary can be found in [31]. In the present case, it is aimed to find an adequate probabilistic model to represent the dependence structure of drought severity and drought duration. One of the advantages of using copulas to represent the interdependence between variables is that copulas display this interdependence in its purest or essential form [32]. A key point in using copulas and a fundamental characteristic is its ability of constructing the dependence structure between random variables independently of the choice of the marginal distributions [31]. Another important and interesting property of copulas is that it can express whether the corresponding dependence is different for different quantiles of the modelled variables [33]. Mathematically, a copula is a function C : [0, 1] d → [0, 1] with uniform marginals. d represents the dimension on which the copula is applied to. In other words, it is defined on the d-dimensional unit hypercube.
A key result which allows an expression of the interrelationship between drought duration and severity in terms of a two-dimensional copula is that any multivariate distribution function can be represented with a copula [34]. In a two-dimensional space, having fitted F 1,α 1 , F 2,α 2 , C θ representing the distributions for duration, severity, and copula, respectively, with parameters α 1 , α 2 , θ, a joint model can be constructed [31].
It is of common use in the literature to define the variables U, V to represent F 1 (x) and F 2 (y) to construct a 2-dimensional model. F 1 (x) and F 2 (y) represents the cumulative distribution function of the random variables X and Y, respectively. In a parametric realm and in order to construct a bivariate copula model, several families can be considered and have been used in recent studies in hydrological applications (see, for example, [35]). These models may depend on one or more parameters.
In order to find the optimum parameter values of the models, the (normalized rank-based) pseudo-likelihood function [31] l θ is used. Given a continuous model C θ with associated density c θ , l θ is given by Equation (2): in which the pairs (R i , S i ) n i=1 represent the ranks of the variables duration and severity and θ represents the parameter (set) of the bivariate model. In order to find the most appropriate kind of dependence structure from the different copula families, the empirical copula (Equation (3)) is used to approximate the theoretical copula in terms of the root mean square error (RMSE): in which 1(·) represents the indicator function. The model for the return period as a function of the joint relationship for drought duration and severity and of the corresponding marginals can be expressed as [37]: which represents the recurrence interval for the event severity S greater than a certain value s and duration D greater than a certain value d. E(L) represents the expected value of arrival time L.
after finding an optimum parameter θ and selecting the most appropriate model (Equations (2) and (3) respectively).

Hydrological Modelling and Projected Changes in Water Availability
The NSE for the calibration and the validation period are 0.90 and 0.84, respectively. The rainfall-runoff model HBV is able to reproduce the peak flow and the low flow conditions very well, similar to results presented elsewhere with NSE varying between 0.49 and 0.83 [38]. The model captures the dynamics of flow and reproduces the seasonal behavior in good agreement with observations in both the calibration and validation period (Supplementary Figure S1). Following, the calibrated model was forced with the bias corrected GCM precipitation and temperature time series from the five GCM models and the four RCPs to project future runoff. Figure 1 shows the percentage change in mean discharge from present to the mid-century (2036-2065) and end-century (2070-2099) for all scenarios (over all GCMs/RCPs). A strong variation among GCM for any given RCP is found, well in agreement with other reports for the Rhine River [3,39,40]. For example, RCP2.6 projects a percentage of change in mean discharge from −5% to +17% during mid-century, reflecting the important uncertainty associated to GCMs. The variation in the percentage in the mean discharge does not increase/change uniformly from low concentration pathways (RCP2.6) to high concentration pathways (RCP8.5). This is observed in all GCM outcomes for both the mid-and end-century. This previous result has also been observed in a number of catchments in Europe [41], in other continents [42,43], and globally [44].

Drought Characteristics and Copula Selection
The SRI is calculated using HBV rainfall-runoff model projections from all scenarios. From this index, drought duration and severity are derived. The statistics of duration and severity for the historical period are given in Supplementary Table S2. The correlation coefficient between drought duration and severity for this period is ρ = 0.94 , suggesting strongly the appropriateness of considering a bivariate model between these two variables.
After fitting each of the five copula models (Supplementary Table S3) through the pseudolikelihood function (Equation (2)), the RMSE associated to the optimum parameters was calculated from the empirical and theoretical copula. Each of the five GCMs was considered. Results for the historical period are given in Table 1. It is observed that there is not a unique copula which best fits in all cases. According to this, the differences between the Frank family and FGM family are rather similar. Nevertheless, the Frank copula performs slightly better in more cases. This is also observed when fitting the models in the projection periods in which the Frank copula performs (slightly) better in most cases. Given this, the Frank copula was adopted to model the bivariate duration-severity dependency in both the historical and projected period. In order to visually check the appropriateness of the selected copula in terms of replicating the observed bivariate pattern, scatterplots of the empirical copula and simulated pairs originating from the best fit model [31] are compared. This is accomplished through the conditional distribution of the fitted copula applied to randomly generated numbers in the unit interval (transformed space). N = 500 simulations are performed. Results for the historical period and end-century projection period for the scenario RCP8.5 and model GFDL-ESM2M are exemplarily displayed in Figure 2. A good agreement is observed between the empirical and theoretical copula. This means, the selected Frank model appropriately represents the bivariate interdependency of the two drought characteristics. The remaining scatterplots are omitted showing the same good agreement. However, a summary is given in Figure 3 in which sets of the same size of the scatterplot are simulated for comparison. This figure shows the theoretical against the empirical cumulative function for all driving models.

Drought Characteristics and Copula Selection
The SRI is calculated using HBV rainfall-runoff model projections from all scenarios. From this index, drought duration and severity are derived. The statistics of duration and severity for the historical period are given in Supplementary Table S2. The correlation coefficient between drought duration and severity for this period is ρ = 0.94, suggesting strongly the appropriateness of considering a bivariate model between these two variables.
After fitting each of the five copula models (Supplementary Table S3) through the pseudo-likelihood function l θ (Equation (2)), the RMSE associated to the optimum parameters was calculated from the empirical and theoretical copula. Each of the five GCMs was considered. Results for the historical period are given in Table 1. It is observed that there is not a unique copula which best fits in all cases. According to this, the differences between the Frank family and FGM family are rather similar. Nevertheless, the Frank copula performs slightly better in more cases. This is also observed when fitting the models in the projection periods in which the Frank copula performs (slightly) better in most cases. Given this, the Frank copula was adopted to model the bivariate duration-severity dependency in both the historical and projected period. In order to visually check the appropriateness of the selected copula in terms of replicating the observed bivariate pattern, scatterplots of the empirical copula and simulated pairs originating from the best fit model [31] are compared. This is accomplished through the conditional distribution of the fitted copula applied to randomly generated numbers in the unit interval (transformed space). N = 500 simulations are performed. Results for the historical period and end-century projection period for the scenario RCP8.5 and model GFDL-ESM2M are exemplarily displayed in Figure 2. A good agreement is observed between the empirical and theoretical copula. This means, the selected Frank model appropriately represents the bivariate interdependency of the two drought characteristics. The remaining scatterplots are omitted showing the same good agreement. However, a summary is given in Figure 3 in which sets of the same size of the scatterplot are simulated for comparison. This figure shows the theoretical against the empirical cumulative function for all driving models.

Drought Events
The number of drought events is calculated for every driving model and RCP. Figure 4 shows this information for the two projection periods 2036-2065 and 2070-2099. Differences in the projections are further quantified and compared for each case.

Drought Events
The number of drought events is calculated for every driving model and RCP. Figure 4 shows this information for the two projection periods 2036-2065 and 2070-2099. Differences in the projections are further quantified and compared for each case. This difference increases to = 17 in the 2070-2099 period for RCP4.5. On the other hand, the uncertainty associated to RCPs is observed to be smaller yet not negligible. By fixing every driving model (i.e., analyzed separately), the maximum difference is estimated to be = 10 for both projected periods and associated to the models IPSL-CM5A-LR and GFDL-ESM2M for the first and second period, respectively.

Duration-Severity Drought Patterns and Return Periods
Return periods associated to drought events were calculated through the fitted Frank copula (Supplementary Table S3). To each return period, all combinations GCM-RCP were considered. Results divided by GCM are summarized in Figure 5 showing the duration-severity relationship for the historical period and expected relationship for the mid-and end-century as a function of return periods 2, 5, and 10 years. For a better visual interpretation, results are also grouped according to RCPs in the Supplementary Figure S2. The number of drought events is observed to be clearly GCM and RCP dependent. For the first projection period, the number of events is expected to vary between n = 22 and n = 37 calculated over all GCMs and RCPs. Similar projections are observed in the 2070-2099 period for which the number of events varies within the range 23 ≤ n events ≤ 40.
Exploring the differences in the projections associated to GCMs (Figure 4), we observed a larger uncertainty in the number of events in the 2070-2099 period, especially for RCP4.5 and RCP6.0 compared to the 2036-2065 period. The maximum difference (d 1 ) between the predicted number of drought events in the 2036-2065 period over all GCMs occurs by RCP4.5 and RCP8.5 with d 1 = 11. This difference increases to d 2 = 17 in the 2070-2099 period for RCP4.5. On the other hand, the uncertainty associated to RCPs is observed to be smaller yet not negligible. By fixing every driving model (i.e., analyzed separately), the maximum difference is estimated to be d = 10 for both projected periods and associated to the models IPSL-CM5A-LR and GFDL-ESM2M for the first and second period, respectively.

Duration-Severity Drought Patterns and Return Periods
Return periods associated to drought events were calculated through the fitted Frank copula (Supplementary Table S3). To each return period, all combinations GCM-RCP were considered. Results divided by GCM are summarized in Figure 5 showing the duration-severity relationship for the historical period and expected relationship for the mid-and end-century as a function of return periods 2, 5, and 10 years. For a better visual interpretation, results are also grouped according to RCPs in the Supplementary Figure S2.
The characterization of drought events indicate that the association between drought duration and severity is case dependent ( Figure 5). We find that both low and high return periods present small and important variations in the bivariate patterns when compared. These variations are clearly GCM-RCP dependent showing different patterns for both mid-century and end-century. In regards to the historical period, Figure 5 shows that the (bivariate) duration-severity interdependence is expected to change in the future.
The differences in the projected bivariate patterns can be visualized comparing two models, namely, GFDL-ESM2M and MIROC-ESM-CHEM. Considering the same 2036-2065 projected period and same return period T r = 10, model 1 (GFDL-ESM2M) shows severity values between 12 and 14 (over all RCPs) for the duration in the range 0 ≤ d ≤ 10 months ( Figure 5). The driving model 2 (MIROC-ESM-CHEM) shows, on the other hand, a clearly larger spread among RCPs with values varying from 11.5 to around 15 for the duration in the same range. The characterization of drought events indicate that the association between drought duration and severity is case dependent ( Figure 5). We find that both low and high return periods present small and important variations in the bivariate patterns when compared. These variations are clearly GCM-RCP dependent showing different patterns for both mid-century and end-century. In regards to the historical period, Figure 5 shows that the (bivariate) duration-severity interdependence is expected to change in the future.
The differences in the projected bivariate patterns can be visualized comparing two models, namely, GFDL-ESM2M and MIROC-ESM-CHEM. Considering the same 2036-2065 projected period and same return period = 10, model 1 (GFDL-ESM2M) shows severity values between 12 and 14 To explore and compare the different projections between GCMs and related uncertainties, Supplementary Figure S2 groups the projected patterns according to RCPs. Important differences are observed between GCMs. In general, higher differences in the projections are estimated for T r = 10 compared to smaller return periods.
Overall, a clear and common tendency across all GCMs cannot be reported. Rather, results indicate there are important differences in the projected bivariate patterns. This is particularly clear for return period T r = 10. Further, uncertainty associated to GCMs exceeds uncertainty attained by RCPs. Similar to these results, other authors have also found larger uncertainties in projected droughts associated with GCMs compared to the uncertainty introduced due to RCPs [45,46].
Additionally to the last point, we underline that runoff projections from impact models for future drought studies are generally highly uncertain [9]. Uncertainty in GCMs is often the most relevant source, overruling the uncertainty introduced by the impact models. In a study on seven large river basins which includes the Rhine River at Lobith, GCM uncertainties mostly dominate over hydrological model uncertainty for the projection of runoff drought characteristics [3]. They also found for the Rhine that the GCM related uncertainty is much higher (more than four times) than that associated with impact (hydrological) models and that the difference grows continuously until the end of the century. In other regions, similar results have been found. Disagreements for GCM projections compared to hydrological model selection for the mid-and end-21st century analyzing runoff in a catchment in the Ecuadorian Andes [47] was also reported. In that study, an ensemble of seven different model structures and eight GCMs was utilized. Similar conclusions were drawn for drier environments in Southeast Australia [48] and in a snow-dominated area in Canada [49]. In a study addressing climate change impacts on European river floods under different warming conditions, [50] concluded that the contribution of GCM uncertainty to overall uncertainty is, in general, higher than the hydrological model-related uncertainty. Moreover, [51] showed that GCM uncertainty dominates over hydrological model uncertainty in Mediterranean and Atlantic regions in Europe for all analyzed warming levels (1.5, 2.0, and 3.0 • C). Hence, the conclusion is that the general patterns found for the variability of projections between GCMs and RCPs are likely to be similar for other impact models, overlying and masking some of projected changes in the mid-and long-term.

Discussion
Changes in the pattern of drought characteristics estimated here as a bivariate relationship can unfold as a result of several causes. For instance, precipitation patterns can dynamically and spatially change. A seasonal shifting and a regionally changing pattern of precipitation has already been shown for different latitudes [52]. Projections of the pattern of water supply in Germany have shown that wetter winters are expected for the next decades based on the SPI analysis. This can lead to the delay of soil drying which might extend into the summer period [53].
It is important to mention that changes in drought characteristics as a result of climate change is a phenomenon that has been reported to be already unfolding due to spatial-temporal changes in climate patterns [54]. Important in this regard are hydrological droughts which can cause devastating impacts on ecological systems and many economic sectors [9] and particularly in agriculture [54]. Although univariate drought analyses and projections are very valuable, we see the additional need to apply multivariate approaches for a more comprehensive understanding and assessment of droughts.
A key issue is the improvement of predictions and projections of drought characteristics and the consideration of both natural variability and anthropogenic climate change [4].
Drought is a complex phenomenon and there are many ways of defining it. This study addressed hydrological droughts. However, a more comprehensive understanding of droughts and their changes over time could be achieved by incorporating other drought definitions and indexes. Furthermore, there are other sources of uncertainties as well influencing droughts as, for instance, hydrological model choice, downscaling approaches, or predictive uncertainty, which are not addressed in this study and could eventually be investigated in other studies.
As previously stated, results presented here indicating a change in the bivariate drought patterns for different return periods brings with it some important implications with regard to occurrences of precipitation such as spatial and temporal (dynamical) changes. This can be translated as, for example, a change in seasonality together with different duration-severity interrelationships. A further study can address these expected temporal changes in precipitation occurrences for future projections as well as changes in projected water generation through different rainfall-runoff models. It is important to mention that the results presented here are valid only in the study region. However, the applied methodology can be performed in other areas and its validity will depend on fitting an adequate (optimal) probabilistic model for the bivariate duration-severity relationship and the calibration and validation of a rainfall-runoff model as carried out here.

Conclusions
Most studies addressing the relation of climate change impacts and droughts perform their analysis through univariate statistical approaches where individual drought characteristics are separately considered. However, droughts have multiple dimensions (e.g., magnitude, frequency, duration) and hence multivariate approaches provide a deeper insight and a more comprehensive characterization, especially for future applications. The present study addressed drought events as a multivariate phenomenon. Bivariate drought patterns were, first, characterized for current climate conditions and, second, for two projected periods, namely, during the mid-(2036-2065) and end-(2070-2099) part of the present century. Given the strong interrelationship of drought duration and severity, we developed a bivariate model by fitting an appropriate copula function. From different models applied, the Frank copula showed better performance in terms of the RMSE and hence was used to model the bivariate duration-severity interrelationship. Different GCM-RCP combinations were studied from which the differences in the joint duration-severity behavior were quantified. The clear identified differences between projections imply significant uncertainties in the projected joint patterns for the two periods.
Main results for the Rhine River can be summarized as follows: first, the joint behavior of drought duration and severity is expected to change in the next decades, but no single pattern amongst GCMs and RCPs can be found for both, mid-and end-century projection periods. Second, the interdependence drought duration-severity shows distinctly different associations depending on the return period, GCM, and RCP. By comparing them, both relatively small and important differences are observed. Third, uncertainty is present to different degrees depending on each particular case analyzed. This is valid for both projected periods studied.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/8/2288/s1, Figure S1: H Observed and simulated monthly discharge for calibration and validation period for the River Rhine at Lobith; Figure S2: Bivariate interdependence duration-severity for return periods 2, 5 and 10 years for the mid-and end-century and five GCMs separated by RCPs; Table S1: HBV model parameters and their meaning, Table S2: Statistics of drought events associated with the observation period 1971-2000, Table S3: Used copula models and associated parameter space.