2.1 Precipitation

We used reconstructed seasonal precipitation (0.5^∘ × 0.5^∘) over Europe (30.25–70.75^∘ N, 29.75^∘ W–39.75^∘ E) from 1500 to 2000. Reconstructed precipitation (P) was derived by Pauling et al. (2006) through principal component regression based on documented evidence (i.e. memoirs, annals and newspapers), speleothem proxy records (Proctor et al., 2000) and tree-ring chronologies from the International Tree-Ring Data Bank (ITRDB; Jeong et al., 2021) .

2.2 Temperature

Reconstructed temperature (T) was obtained from Luterbacher et al. (2004), which relies on historical records and seasonal natural proxies (i.e. ice cores from Greenland and tree rings from Scandinavia and Siberia). Reconstructed temperature data are available at the same spatial and temporal resolution as precipitation (see Table 1). We refer to both of these datasets as reconstructed forcings or reconstructed precipitation/temperature fields.

2.3 Self-calibrating Palmer Drought Severity Index (scPDSI)

In addition, we used data from the Old World Drought Atlas (OWDA; Cook et al., 2015) which contains information regarding moisture conditions across Europe, specifically the self-calibrated Palmer Drought Severity Index (scPDSI) using summer-related tree-ring proxies over the period 0 to 2012 CE.

2.4 The Global Historical Climatology Network (GHCN)

The GHCN dataset (GHCN; Peterson and Vose, 1997) is one of the largest observational databases, collated by the National Oceanic and Atmospheric Administration (NOAA; Quayle et al., 1999). The GHCN-m dataset contains observed temperature, rainfall and pressure data from 1701 to 2010. Data for the majority of stations are, however, available after 1900. GHCN-m precipitation and temperature from GHCN V2, as well as from the GHCN V4 version (Menne et al., 2012), were used to assess the reconstruction accuracy of the P and T fields as an input into the considered models. We selected 113 precipitation and 144 temperature stations within the European domain (see Fig. 1) with records dating back earlier than 1875. Most stations are geographically concentrated in central Europe, and few stations are located in the eastern and northern areas of Europe (see Table 2). These data, hereafter, are referred to as the GHCN data.

2.5 Observed runoff

The Global Runoff Data Center (GRDC; http://www.bafg.de/GRDC/EN/Home/homepage_node.html, last access: 24 November 2016) provides data for more than 2780 gauging stations in Europe, with the oldest records starting from 1806. Only the GRDC runoff time series with at least 25 years of data prior to 1900 were selected. In total, there were 21 such stations predominantly available in central Europe: 11 in Germany, 2 in France, 2 in Switzerland, 1 in the Czech Republic, 1 in Sweden, 1 in Finland, 1 in Lithuania and 1 in Romania (see Fig. 1). These stations cover 12 European river basins (Rhine, Loire, Elbe, Danube, Wesser, Main, Glama, Slazach, Nemunas, Gota Alv, Inn and Kokemaenjoke), with areas ranging from nearly 6100 km² (Kokemaenjoki, Muroleenkoski, Finland) to 576 000 km² (Danube, Orsova, Romania). The mean annual discharge (Q_mean) varies from 50 to 5 600 m³ s⁻¹ and spans different time periods for each catchment.

Table 2Selected study catchments.

Download Print Version | Download XLSX

The most extensive records were available in Sweden (Vargoens KRV) and Germany (Dresden), containing the longest discharge series of 212 and 208 years, respectively. The gauging station in Köln also provided 195 years of data for the Rhine River. Note that some of the gauging stations are located nearby and therefore have a greater degree of similarity in their runoff time series (e.g. two stations in Basel, Rhine). Detailed information relating to all selected stations is provided in Table 2.

2.6 Study area

In the first part of the study, the grid-based reconstruction of precipitation and temperature was verified against the available GHCN data across the European region bounded by (30.25–70.75^∘ N, 29.75^∘ W–39.75^∘ E). The second part focused on 21 specific central European catchments, corresponding to the available long-term GRDC discharge records. The study area and the observational data of the hydroclimatic variables are shown in Fig. 2.

3.1 Data preprocessing

Two databases were considered for the analysis and development of the annual runoff reconstruction. The first one was used for evaluating the accuracy of meteorological forcing reconstructions used for hydrological simulations and consists of observed GHCN data for all available European stations with long records (see Sect. 2.4) and values of corresponding grid cells from the reconstructed forcings dataset.

Figure 2A schematic representation of workflow carried out in the study.

Download

The second database was created as the basis for runoff reconstruction containing the observed runoff data for 21 selected catchments (Table 2) and the corresponding input variables of the models used to generate the multi-century runoff reconstructions. Several input variables were considered for inclusion in models such as reconstructed precipitation and temperature and Old World Drought Atlas scPDSI. The catchment average precipitation, temperature and scPDSI were estimated from the corresponding (gridded) datasets by averaging the relevant grid cells over the catchments. This database was further divided into two parts, calibration (1900–2000) and validation (before 1900), to assess the model's accuracy and to select an appropriate model. The data preprocessing, model selection, and evaluation of the models are depicted in Fig. 2.

3.2 Hydrologic model (GR1A)

We applied the annual timescale hydrologic model, GR1A (Mouelhi et al., 2006), to simulate annual runoff in each catchment. GR1A is a conceptual lumped hydrologic model (Manabe, 1969), considering dynamic storage and antecedent precipitation conditions. The model consists of a simple mathematical equation with a single (optimized) parameter:

where Q, E and P represent annual runoff, basin average potential evapotranspiration and basin average precipitation, respectively, and i denotes the year. The parameter X is optimized individually for each catchment by maximizing the Nash–Sutcliffe efficiency (NSE) between observed and simulated runoff. Default gradient-based optimization from the R package airGR (Coron et al., 2017) was used. The potential evapotranspiration was calculated using the temperature-based formula (Oudin et al., 2005). Compared to other conceptual models from the GR family (GR4J, GR5J), GR1A is simple to use, and it allows for analysing many variants, particularly defining the best antecedent rainfall, and is potentially useful to predict wet and dry hydrologic conditions (Mouelhi et al., 2006).

3.3 Data-driven models

Artificial neural networks (ANNs; Senthil Kumar et al., 2005; Kwak et al., 2020) can describe non-linear relationships and are widely used for rainfall–runoff prediction. The ANNs consist of artificial neurons organized in layers and connections that route the signal through the network. Each connection has an associated weight that is optimized within the calibration (in the context of ANNs, known as training). There are many types of ANNs which differ in terms of structure and type of connections, as well as direction and functional forms used for neuron activation or training. In the present study, we considered two approaches: long short-term memory (LSTM) neural networks and Bayesian regularized neural networks (BRNNs). These approaches have been commonly used in past rainfall–runoff modelling studies (Hu et al., 2018; Kratzert et al., 2018; Xiang et al., 2020; Ye et al., 2021). We considered combinations of reconstructed forcing, OWDA-based scPDSI and lagged forcing as an input into the network for both model types. Specifically, the network using only reconstructed precipitation and temperature fields is referred to as [P,T], the network with reconstructed forcing and OWDA scPDSI is termed as $[P,T,\text{PDSI}]$; and finally the network which includes 1-year lagged P and T forcing in addition to actual P and T is referred to as $[P,T,\text{lag}]$. We also considered and explored lag times longer than 1 year. However the correlation between precipitation and runoff drops significantly at lag times longer than 1 year and therefore was not included in presented analysis.

Figure A1 shows the architecture of LSTM, which is a modified version of the recurrent neural network (Hochreiter and Schmidhuber, 1997), using backpropagation in time (Werbos, 1990). LSTM is known for efficient simulation of time series with long-term memory (Van Houdt et al., 2020). It generally consists of two unit states (hidden and cell states) and three distinct gates (hidden, input and output). In this process, the cell state saves the long-term memory at the previous unit, while hidden states act as a working memory to process information inside the gates. These gates can determine which information needs to be processed, remembered and transferred in the next state. With LSTM, different activation functions, such as hyperbolic tangent and sigmoid, can be used to update unit states. The implementation of the LSTM is carried out by applying the R packages “keras” (Arnold, 2017) and “tensorflow” (Abadi et al., 2016).

The training process of the LSTM is time-consuming due to its inherent complexity. Therefore we also considered the BRNN models that provide fast learning and convergence and were already used to tackle the complex relationship between rainfall and runoff (Ye et al., 2021). BRNNs are based on recurrent neural networks, which are often used to model time-series data (Wang et al., 2007), and the methods are extended with Bayesian regularization (Okut, 2016) to account for uncertainty related to network parameters and input data (Zhang et al., 2011). We trained this model in R using the “brnn” function of the “caret” package (Kuhn, 2015). More details are available in Appendix A3.

To set the optimal hyperparameters of the models (such as the number of neurons and activation functions) and to reduce the likelihood of overfitting during the calibration/training, the model performance was cross-checked considering an independent (or so-called “testing”) set. The testing set was for each learning exercise extracted from the calibration data (1900–2000) as a random fraction (25 %). This process of the model development was repeated several times, minimizing the root mean square error (for BRNN) and mean square error (for LSTM) for each catchment individually. The model with the best performance was then chosen for further evaluation.

3.4 Goodness-of-fit assessment

We used a set of seven statistical metrics to assess the performance of simulated runoff, namely Nash–Sutcliffe efficiency (NSE), Pearson correlation (R), standard deviation ratio (rSD), Kling–Gupta efficiency (KGE), root mean square error (RMSE), mean absolute error (MAE), bias (BIAS) and relative bias (relBIAS). The mathematical formulations of these metrics are provided in Appendix A1.

3.5 Runoff drought identification

To check the utility of our reconstruction, we finally explore how well the annual runoff droughts are represented in the simulations. Our study considers annual hydrological droughts, defined as the streamflow/runoff deficit, following the threshold level approach (Yevjevich, 1967; Sung and Chung, 2014; Rivera et al., 2017). This approach is typically used for daily or monthly timescales, considering 0.1 or 0.2 quantile threshold levels. To accommodate the annual scale used here, we defined the start of the drought, when the annual runoff anomaly falls below the 0.33 quantile (regular drought) and the 0.05 quantile (extreme drought). The drought persists until the runoff rises above the threshold again. After that, the difference between runoff and the threshold was determined for each identified drought year, called the runoff deficit. Hydrological drought series can be further assessed to understand the critical aspects of runoff (temporal) dynamics and to classify past droughts in Europe (Wetter and Pfister, 2013; Cook et al., 2015).

4.1 Evaluation of reconstructed precipitation and temperature fields

The 500-year annual paleoclimate reconstructions of precipitation (P) and temperature (T) were validated against the GHCN observation data. The map showing the comparison is given in Figs. 3 and 4. The reconstructed data are evaluated against observational P and T across 99 and 94 European sites, respectively. Figure 3 shows that for most of the sites the correlation coefficient (R) of P reconstruction at most of the sites is above 0.5; the relative bias (relBIAS) is between −0.1 and 0.1; KGE and NSE are showing values below 0.5 and 0.6, respectively; the rSD is between 0.7 and 0.9; and RMSE varies between 0 and 150.

The performance of the temperature reconstruction was relatively better, as depicted in Fig. 4. In this case, RMSE between reconstructed and observational T is around 0.2 ^∘C, and rSD fluctuates between 0.95 and 1.05, while R is higher than 0.84, and BIAS is less than 0.5 ^∘C, except for stations located in the Alps. The NSE and KGE values were above 0.5 at the majority of the stations. Low skill observed at some locations can be explained by the unresolved variability of grid-cell average temperature, especially in regions with complex terrain.

It is worth noting that the large spread of goodness-of-fit (GOF) statistics is mainly due to the outlying values at the grid cells located along the boundary of the domain (i.e. the interface between land and sea/ocean) and high elevations (see also Figs. 3 and 4). In general, reconstructed precipitation exhibits greater differences from observations than temperature. This may be because the proxies considered in the reconstruction rely on different seasons and climate conditions. Additionally, the shortest available instrumental data before the 20th century could encounter certain technical errors, such as problems with instrumental tools, station relocation and dating issues (Dobrovolný et al., 2010). Moreover, other studies (e.g. Ljungqvist et al., 2020) stated that the precipitation series employed for the reconstructions were relatively shorter and more erroneous than the temperature series before the 20th century (Pauling et al., 2006; Harris et al., 2014). Finally, the chosen statistical technique (principal component regression) could also possibly contribute to variance inflation with larger timescales (Pauling et al., 2006).

Figure 5The correlation coefficient (top) and NSE (bottom) for calibration (left) and validation (right) of the considered models for 21 study catchments. The vertical axis represents the catchments (station name and river) and the horizontal axis the considered models. The rectangular black frames represent the catchments with NSE > 0.5 over the validation period.

Download

4.2 Assessment of the reconstructed runoff simulations

The GR1A conceptual hydrological model was driven by catchment average P and T and calibrated using observed annual runoff for each catchment separately. The simulated annual runoff series were then compared to the corresponding GRDC observations (for calibration and validation periods), and the results were summarized by means of GOF statistics. As can be seen in Fig. 5, the correlation and NSE statistics for calibration achieve reasonable results at most of the catchments, with a few exceptions (i.e. Kokemenjoki, Goeta, Nemunas and Inn). The catchments with relatively poor skills are located in northern Europe, which is in line with the previous findings by Seiller et al. (2012), who noted that the lumped hydrological models often exhibit larger uncertainties and fail to capture the extreme catchment values (both high and low) in those regions. The low skill for some of the catchments cannot be easily attributed only to bias in reconstructed precipitation and temperature (described in Sect. 4.1) but rather to low station and proxy coverage in some (especially northern) parts of Europe, leading to biased basin-average precipitation and temperature estimate. Another study of Fathi et al. (2019) suggested that the performance of the GR1A model is less efficient than the new Budyko-framework-based SARIMA model in simulating the annual runoff across the Blue Nile and the Danube catchment. This may be due to the simplified nature of the model that does not easily capture the complex relationship between rainfall and runoff variability.

In general, statistical values presented in heat maps (Fig. 5) indicate that the neural network algorithms are more skilled for runoff prediction than the GR1A model. The NSE and R statistics for the BRNN and LSTM models indicate a significant improvement in runoff prediction, as compared to the results obtained through the GR1A model. For instance, for Basel Rheinhalle the NSE increases from 0.27 to 0.73 (BRNN) and 0.75 (LSTM) for calibration and 0.2 to 0.54 (BRNN) and 0.52 (LSTM) for validation. Moreover, including scPDSI from OWDA with reconstructed forcing $[P,T,\text{PDSI}]$ increases the performance slightly more (NSE 0.76 for calibration and 0.57/0.59 for validation, for BRNN/LSTM, respectively), and considering the lagged forcing results in the best performance (NSE 0.75/0.8 for calibration and 0.6/0.54 for validation, for BRNN/LSTM).

Similarly for all sites, the data-driven methods exhibited a strong correlation with the observed runoff, with the GR1A simulations resulting most frequently in lower correlation values. Other metrics (RMSE, MAE, KGE, rSD and relBIAS) are shown in Figs. S1–S5 in the Supplement. Across many study locations, the combination of reconstructed forcings and their 1-year lag performed the best in terms of rapid convergence (the number of iterations needed) and high accuracy from all input combinations for both data-driven models (BRNN, LSTM). For the validation period, the mean NSE (across all catchments) for the GR1A model is 0.16, for the BRNN $[P,T,\text{lag}]$ it is 0.68 and it improves to 0.73 for the LSTM $[P,T,\text{lag}]$. In the case of the mean KGE, GR1A yields 0.62, BRNN $[P,T,\text{lag}]$ is 0.73 and LSTM $[P,T,\text{lag}]$ is 0.78.

Figure 6Comparison between the models for the station with the best (Blois–Loire River, top) and the worst (Smalininkai–Nemunas River, bottom) model fit.

Download

To further demonstrate the differences between the individual models, we show the simulated runoff series for all models for those catchments with the highest (Blois–Loire) and lowest (Smalininkai–Nemunas) performance in Fig. 6. The performance of the models is comparable during the calibration period for the Loire River. Clearly, all data-driven models are capable of mimicking the observed runoff, while the GR1A model exhibited certain minor deviations, primarily until 1930. In the validation period, the differences between the models are more visible, in particular, for above-average flows. This can be attributed to different generalization skill of individual models. At the beginning of the validation period (1870–1880), all models failed to simulate the high annual flows.

In the case of Nemunas catchment, the GR1A simulation deviates extremely from the observed data and cannot capture the mean flow level. However, the calibration is poor, even for the data-driven models, and does not simulate the year-to-year variability appropriately. Interestingly, for the validation period, the error in the GR1A model decreases. The performance of the data-driven models is similar in validation and calibration periods. Looking at the GOF statistics, the models considering OWDA-based scPDSI or lagged forcings (e.g. P_t−1) perform slightly better in terms of KGE than the other model configurations.

4.3 The annual runoff reconstruction datasets

As a first step, we excluded the catchments that exhibited poor performance in validation (see Fig. 5). As a threshold, we considered validation NSE greater than 0.5 for at least one model, following the approach used by Ayzel et al. (2020). In this step, we excluded 7 catchments (Vlotho–Wesser, Decin–Elbe, Burghausen–Salzach, Smalininkai–Nemunas, Vargoens KRV–Goeta, Elverum–Glama and Muroleekoski–Kokemenjoki) out of 21, ending up with a set of simulations for 14 catchments (highlighted by the rectangular box in Fig. 5).

Secondly, we identified the best candidate models for each of the 14 selected catchments, considering the GOFs based on the validation NSE and R greater than 0.5 and 0.7, respectively. The best model for each catchment was finally selected from those models considering the remaining validation measures (relBIAS, rSD, KGE, RMSE and MAE) as well. Specifically, we picked the models with consistent good validation measures. This choice is partly subjective, and more formal selection should be explored further. On the other hand, the candidate models all performed comparably in most cases.

Figure 7Distribution functions for BRNN $[P,T,\text{lag}]$, LSTM $[P,T,\text{PDSI}]$, i.e. the best two models, GR1A [P,T]- and GRDC-observed data for the periods 1500–1800, 1800–1900 and 1900–2000 over the Basel Rheinhalle–Rhine catchment. The values on the horizontal axis are transformed using the “probit” function. The coloured labels indicate the most extreme drought years according to each model.

Download

The resulting selected models are shown in Table 3. The combination of reconstructed forcing with 1-year time lags results in the best performance over nine catchments, of which seven employed the BRNN and the remainder the LSTM model. The LSTM with reconstructed forcing and OWDA-scPDSI was the best in just one case, and the remaining time-series reconstructions were most appropriately simulated with the BRNN [P,T] and BRNN $[P,T,\text{PDSI}]$. It should be noted that the differences between the models performing well are small, as noted in Fig. 6 and further demonstrated in Fig. 7. The latter figure compares the cumulative distribution functions of annual runoff for the periods 1500–1800, 1800–1900 and 1900–2000, as simulated by the BRNN $[P,T,\text{lag}]$ and LSTM $[P,T,\text{PDSI}]$ – the two best-performing models – and the GR1A (the most deviating simulation from the best model) with the distribution of the observed annual runoff for the Basel Rheinhalle–Rhine catchment. For the calibration period (1900–2000) in Fig. 7, the models perform well except the GR1A, which generally overestimated the observed maxima. The cumulative distribution of BRNN- and LSTM-simulated runoff values is very similar for the validation period (1800–1900), except for the top and bottom 5 % in 1500–1800. The GR1A simulation showed significant differences for the entire distribution, thus overestimating/underestimating the maxima/minima. Our finding shows that GR1A simulates a Rhine minimum of 279 mm yr⁻¹ in Basel, whereas the observed minimum in the past century is greater than 532.6 mm yr⁻¹, inferring that the cumulative distribution function (CDF) has significantly lower/higher runoff values between 1500 and 1800 for BRNN and GR1A, whereas LSTM appears to extrapolate less. The difference from the best model can be expressed in terms of KGE – even here, it was evident that the GR1A model deviated considerably (KGE 0.6–0.7), while the LSTM is very similar to the BRNN (KGE 0.92–0.96). The most severe drought year identified by the models in the period 1500–1800 appears to be 1669 and the year 1921 in the past century (1900–2000) (Fig. 7 left and right panels), while for 1800–1900 the models identified either 1865 (GR1A, LSTM) or 1858 (BRNN). Please note that the 1858 low-water mark is available at Laufenburg Pfister et al. (2006) near Basel and was regarded as one of the worst winter droughts in the last 200 years.

Figure 8Reconstruction of runoff series for Köln–Main and Hann. Münden Wesser rivers. The blue line corresponds to the reconstructed series, and the black and red lines represent the observed runoff for the calibration and validation period, respectively.

Download

Figure 9Observed and simulated runoff for 14 selected catchments in the calibration and validation periods. The solid line represents the 1:1 relation, and the dashed line corresponds to fitted regression between observed and simulated runoff.

Download

Figure 10The observed and simulated runoff deficit based on the 33rd percentile threshold for 14 selected catchments during the calibration and validation period. The solid line represents the 1:1 relation, and the dashed line corresponds to fitted regression between observed and simulated runoff.

Download

The resulting 14 annual runoff reconstructions are available at https://doi.org/10.6084/m9.figshare.15178107 and are shown in Figs. S6–S8 in the Supplement. As an example, we present only two runoff reconstructions here (Fig. 8). As an additional validation for the reconstructed series, we inspected the scatter plots of the observed and reconstructed runoff (Fig. 9). The simulated series are generally consistent with the observed runoff, especially for the Montjean–Loire, Köln–Rhine and Basel Schifflaende–Rhine catchments, which exhibit the best relationship between the observed and the simulated runoff.

Finally, to check the consistency of our reconstructed dataset, we compared the skill of our simulation with respect to the GRDC runoff observation and the GSWP3-forced GRUN monthly runoff (Ghiggi et al., 2019) datasets. The gridded GRUN datasets were averaged over the respective catchments to enable comparison (Figs. S9 and S10 in the Supplement). Our reconstruction outperforms GRUN data in terms of RMSE, MAE, relBIAS and NSE across the majority of the catchments, while the correlation (reproduction of interannual dynamics) to GRDC runoff is slightly higher for GRUN compared to our reconstruction. The variability, which our data-driven models underestimate (on average by 16.5 %), is overestimated by GRUN (on average by 17.2 %). Since the correlation compensates for the relBIAS, the KGE for our reconstruction and GRUN is comparable. This suggests that GRUN could be used for data-driven model training, provided at least some information on flow characteristics is available in the catchment.

Table 3Selection of best model for runoff in individual catchments.

Download Print Version | Download XLSX

4.4 Identification of low flows, significant hydrological drought events and trends

In the final step of the analysis, we compared the droughts identified in the reconstructions with the GRDC-observed series (Fig. 10). The agreement between the simulated and observed runoff deficit is lower compared to the annual runoff time series. For most of the stations, the simulated deficit is lower than the corresponding observed estimates. This suggests that the reconstructed precipitation and temperature fields do not represent the inter-annual variability correctly. Despite a widespread issue with the representation of inter-annual persistence, Fig. 10 shows that the runoff deficits are simulated reasonably well for the Rees–Rhine and Köln–Rhine catchments.

Table 4Simulated runoff droughts since 1500. Years in bold indicate extreme droughts below 5 % quantile.

Download Print Version | Download XLSX

Furthermore, we contrasted reconstructed drought patterns over the last 500 years with data available from documentary evidence and other sources. In the case of extreme droughts, we considered the q_0.05 threshold before 2000 CE. Low-flow analysis since 1500 and the large deficit values for catchments (below 5th percentile) are shown in Table 4. In the 16th century, the years 1536, 1540 and 1590 are associated with significant runoff deficits. The event of 1540 has already been reported (Brázdil et al., 2013; Cook et al., 2015; Brázdil et al., 2019) as the worst event of the 16th century and more severe in terms of changing hydrologic conditions. In 1540, almost 90 % of the Rhine and Elbe River catchments (Basel and Cologne) experienced low yearly discharge, which ranked as the greatest low flows in the last 5 centuries (Leggewie and Mauelshagen, 2018). The seasonal precipitation was also deficient and was evident primarily in central Europe and England (Dobrovolný et al., 2010). Wetter and Pfister (2013) stated that the spring and summer of 1540 were likely to have been warmer than the comparable period during the 2003 drought. The simulation shows that the drought during 1540 was evident in most study catchments, such as the Rhine, Main, Wesser, Loire and Danube, except Wasserburg–Inn.

In the 17th century, the years 1603, 1616, 1631, 1666, 1669, 1676, 1681, 1684 and 1686 were simulated as exceptionally low-flow years. Furthermore, two events (1669 and 1686) were associated with the largest water deficit across several study catchments. Basel Schifflaende–Rhine catchment is a good example of this, which appears to have experienced an extreme runoff deficit during 1669. In the Köln–Rhine catchment, 26 remarkable droughts have been captured over the past 500 years, and the year 1686 reached the largest runoff deficit (156 mm yr⁻¹). The 1616 is considered the driest year of the 17th century, the so-called “drought of the century” (Brázdil et al., 2013), which significantly impacted the major rivers in Europe (e.g. Rhine, Main and Wesser). Brázdil et al. (2018) identified three unusual drought periods (1540, 1616 and 1718–19) over the Czech lands, highlighting the 1616 drought, which caused widespread famine, dried up the Elbe river watershed and altered the climate of neighbouring nations (Switzerland and Germany). The hunger stone of the Elbe River also revealed the exceptionally dry year of 1616 (Brázdil et al., 2013). During the 18th century, a similar level of runoff deficit was simulated in the years 1706 and 1719.

During the 19th century, the years 1863, 1864, 1874, 1893 and 1899 were recognized as drought years in all catchments, while in the 20th century, the driest periods occurred in 1921, 1934, 1949 and 1976. The 1921 drought in the Blois–Loire, Rees–Rhine, Köln–Rhine, Orsova–Danube, Basel Rheinhalle–Rhine and Basel Schifflaende–Rhine catchments was ranked as the most exceptional drought in the 20th century. Three catchments (Basel Rheinhalle–Rhine, Basel Schifflaende–Rhine and Blois–Loire) exhibited a large runoff deficit during the year 1921. A noticeable increase in temperature was experienced across Europe, and certain areas were notably affected by a heatwave in July of that year. The majority of central Europe, southern England and Italy were affected by this drought, where the rainfall was found to have decreased around 50 % to 60 % relative to the average (Bonacina, 1923; Cook et al., 2015). The precipitation totals were recorded as the lowest since 1774, and the year was also ranked top (in terms of deficit rainfall) in the Great Alpine region (Haslinger and Blöschl, 2017), where the rainfall deficit began in winter 1920/21 and lasted until autumn 1921. Also reported in newspapers, the Rhine River (Switzerland), Molesey Weir on the Thames River (United Kingdom) and Loire River (France) all had low river flows in 1921 (Van der Schrier et al., 2021). Monthly runoff anomalies analysed from the GRUN dataset (Ghiggi et al., 2019) show that August 1976 was the fifth driest month between 1900 and 2014, in agreement with some of our catchment reconstructions signalling the 1976 as a yearly drought in the Köln–Rhine, Hann. Münden–Wesser and Bodenwerder–Wesser.

In summary, the reconstructed annual runoff corresponded well to the majority of extreme drought years (e.g. 1540, 1616, 1669, 1710, 1724 and 1921, as highlighted in Table 4) and previously demonstrated in the OWDA-based PDSI tree-ring reconstructions and previous works (Dobrovolný et al., 2010; Brázdil et al., 2013; Wetter and Pfister, 2013; Cook et al., 2015; Markonis et al., 2018). It is important to note that the presented runoff reconstructions might have missed notably documented dry events, e.g. 1894 (Brodie, 1894), which was associated with unprecedented low levels of rainfall and excessive temperature rises in the south of England, the British Isles and other European regions (Brodie, 1894; Cook et al., 2015; Hanel et al., 2018).

Finally, we assessed the linear trends in the decadal runoff series for several time periods. The reconstructed annual runoff for 1500–2000 for each catchment was first aggregated to 10-year averages and divided by the mean annual runoff. The resulting series are shown in Fig. A2. Although significant negative trends were found for all catchments except for one considering the whole 1500–2000 period, the signal is not clearly linear. Instead, for a number of catchments, there is a period of sustained above-average (Orsova–Danube and Dresden–Elbe) or below-average (Blois and Montjean Loire) annual runoff during approx. 1600–1800, while for the rest, the persistence is weaker, although a low runoff signal is still visible (Basel Rheinhalle, Basel Schifflaende and Köln Rhine). When only the last 50 years is considered, the trends are significantly negative (positive) for seven (two) catchments, with the rest being insignificant.