Projections of runoff in the Vistula and the Odra river basins with the help of the SWAT model

Climate change, observed and projected, is not limited to the ubiquitous warming. In many areas, change of atmospheric precipitation is even more important, impact-wise. In much of Central and Eastern Europe, mean renewable, and easily available, surface water resources and river runoff in particular, are rather low. Hence, there is considerable concern that they can be adversely affected by climate change. However, while it is a robust, model-based finding that precipitation increases in the north of Europe and decreases in the south of Europe (Jacob et al. 2014), both annually and in summer, changes in Central and Eastern Europe are more complex.

The climate modelling community has been developing advanced general circulation models (GCMs) and regional climate models (RCMs) for many years now. In Europe, a major milestone has been achieved with launching new, high resolution, EURO-CORDEX projections (Jacob et al. 2014). The hydrological modelling community has been taking up the newest generation of available climate scenarios and updating their projections of runoff, droughts, floods, etc. Application of hydrological models for climate change impact assessments adds value to the existing climate projections, since simulated hydrological response is spatially and seasonally diverse. The models enable tracking of different hydrological processes, such as snow melt, evapotranspiration, soil water dynamics and runoff, that are often separated into sub-surface (slow) and surface (quick) components. However, due to a huge workload needed to process the data (e.g., perform the bias correction), develop and run models, there are still very few hydrological assessments of climate change impacts using these forcing data. To our knowledge, only two large-scale studies, covering Poland, using EURO-CORDEX data have been published so far (Papadimitriou et al. 2016; Roudier et al. 2016) for low flows and two studies (Alfieri et al. 2015; Papadimitriou et al. 2016) for floods. The EURO-CORDEX data were also used in the CHIHE project (Romanowicz et al. 2016) for studying future impacts on droughts (Meresa et al. 2016) and floods (Osuch et al. 2016) in small Polish catchments.

The objective of this paper is to assess climate change effects on spatiotemporal changes in annual and seasonal runoff and its components in the basins of two large European rivers, the Vistula and the Odra, for two future time horizons. This study makes use of the Soil and Water Assessment Tool (SWAT) model, set up at high resolution, and driven by an ensemble of nine bias-corrected EURO-CORDEX simulations under two representative concentration pathways (RCPs), 4.5 and 8.5. It also involves an assessment of robustness of runoff projections. To our best knowledge, it is the first study of this type covering and tailored for this region, filling the ‘scale gap’ between the global or continental climate change impact studies and small catchment studies.

The Vistula and the Odra basins (VOBs) are located in the temperate climatic zone in Central and Eastern Europe and drain into the southern Baltic Sea (Figure 1). Poland constitutes 88% of the entire VOB that covers 312,873 km². The area of the Vistula basin is 193,831 km² and 87% is in Poland and for the Odra the area is 119,041 km² with 88% in Poland. The remaining parts of the Vistula and the Odra basins are located in five neighbouring countries: Czech Republic and Germany for the Odra, and Ukraine, Belarus and Slovakia for the Vistula. Both large rivers have their sources in the mountainous areas in the south of the basins. The great majority of the drainage areas of both rivers extend through the European Plain. The VOB has cold winters and warm summers, with east–west temperature gradient characterized by more continental influence in the east and more maritime influence in the west.

The SWAT model is a process-based, semi-distributed, continuous-time hydrological model that simulates the movement of water, sediment and nutrients on a catchment scale with a daily time step (Arnold et al. 1998). This model has been used, worldwide, in a plethora of water management and climate impact applications (e.g., Krysanova & Srinivasan 2015; Awan et al. 2016; Tamm et al. 2016). A total number of 2,633 sub-basins and 21,311 hydrological response units were distinguished in the model setup for the VOB. In this application, SWAT was driven by the CHASE-PL Forcing Data: Gridded Daily Precipitation & Temperature Dataset 5 km (CPLFD-GDPT5; Berezowski et al. 2016). This product was constructed using geostatistical techniques based on daily data from over 700 stations (for precipitation) and is, to our knowledge, the best alternative for distributed modelling in this region (cf. Piniewski et al. 2017b). Since only daily precipitation and minimum and maximum temperature were available for the studied area, we selected the Hargreaves method to model potential evapotranspiration (PET) in SWAT.

A large and representative set of 80 small, non-nested catchments with relatively unimpaired streamflow was selected for spatial calibration. Additionally, 30 gauges were selected for spatial evaluation of the model. Both calibration and evaluation were done using daily, continuous, flow data. The model performance was assessed as satisfactory, with the median Kling-Gupta efficiency equal to 0.7 and 0.62 in the calibration (1991–2000) and validation (2001–2010) periods, respectively. The use of a high number of gauges in calibration and validation, representing different landscape settings, climate conditions and flow regimes, justifies applying the model for studying spatial variability of runoff projections at high resolution.

Since calibration was based on a set of catchments with relatively undisturbed flow conditions, and since water management was, for purpose, neglected in the model setup, it is expected that the model driven by the future climate data will reflect pure effects of climate change, neither amplified, nor damped, by any other effects. For more details of the SWAT model, as used in this study, see Piniewski et al. (2017b), the study laying the groundwork for the present paper.

The multi-model ensemble (MME) of climate projections used in this study, the CHASE-PL Climate Projections: 5-km Gridded Daily Precipitation & Temperature Dataset (CPLCP-GDPT5; Mezghani et al. 2016) consists of nine RCM simulations (combinations of GCMs and RCMs) from the EURO-CORDEX (cf. Jacob et al. 2014) ensemble (Supplementary material, Table S1, available with the online version of this paper), assuming two RCPs, 4.5 and 8.5. RCPs are plausible pathways towards reaching specific target radiative forcing trajectories (Moss et al. 2010). Two RCPs used in this paper, corresponding to the 4.5 W m⁻² and 8.5 W m⁻² levels of radiative forcing in year 2100, seem to be the most widely used in impact studies worldwide.

The original EURO-CORDEX simulations were bias corrected onto a 5 by 5 km grid based on the CPLFD-GDPT5 (Berezowski et al. 2016) as reference data using the quantile mapping method. All bias-corrected values of climate parameters of concern were available for the following three time slices: 1971–2000, 2021–2050 and 2071–2100. In this study, the first three years of each period were truncated, since a warm-up period of three years is used for SWAT simulations. The corresponding time horizons will be hereafter referred to as ‘historical period’, ‘near future’ and ‘far future’, respectively.

As found by Piniewski et al. (2017a), projections of changes in annual mean of daily minimum (T_min) and maximum temperature (T_max) based on the areal-means of the MME mean are robust and are expected to increase by 1.2–3.7 °C for T_min and by 1.0–3.3 °C for T_max, depending on the future projection horizon and the RCP. Changes accelerated with time and with the RCP value, i.e., higher values of changes were projected for RCP 8.5 than for 4.5. Higher increases were consistently associated with T_min and the gradient of change goes from north-east to south-west. The magnitude of temperature (if not specified as minimum or maximum, denoted as T) increase varied seasonally, with the highest increase projected in winter, for all analysed combinations of parameters, horizons and RCPs. The differences in projected changes between other seasons were generally lower. A higher increase in winter temperatures was projected in comparison to other seasons, notably in the far future under RCP 8.5 (by 1.3 °C).

Projections of annual and seasonal precipitation (P) totals based on the MME mean (Piniewski et al. 2017a) showed much less consistent signals than those of temperature, even though the climate models agree well on the direction of change, consistently indicating increase of the areal mean P. As with temperature, changes were greater for later projection horizon and higher RCP, varying from 5.5% in the near future under RCP 4.5–16.2% in the far future under RCP 8.5. Spatial variability was substantial, although quite variable between individual climate model simulations. Albeit seasonal P totals are projected to considerably increase in all four combinations of RCPs and projection horizons for winter and spring, high model spread reduces considerably the robustness, especially for the far future. In most cases, P increase in winter and spring (7.3–27.9%) was significantly higher than in summer and autumn (1.4–9.1%). The MME members agree well that, overall, the summer and autumn P will not undergo statistically significant changes.

A complete overview of annual and seasonal projections of temperature and precipitation, including three ensemble statistics: mean, minimum and maximum, is presented in Table 1.

The main output of the SWAT model, interesting in the context of this paper, is sub-basin-level runoff (R), i.e., the fraction of P that ultimately reaches the stream network. Total runoff is composed in SWAT of surface runoff, lateral (sub-surface) flow and baseflow and is then routed via stream network as streamflow (Q) volume (Arnold et al. 1998).

For each combination of climate model with RCP and period, and for each sub-basin, daily R and Q values simulated with the help of SWAT were then aggregated to monthly, seasonal and annual totals. Projections of future R were calculated as the relative change between the values in the respective future horizons and the control interval. For the purpose of map presentation, the change factors were calculated at mean annual and seasonal levels.

The MME mean change was used to aggregate the results driven by nine climate model simulations into a single map. Since projections from different models may considerably diverge, presenting exclusively a map of the MME mean or median change is not sufficient as it downplays the uncertainty. When analysing climate (impact) projections from multiple models for a given location, two properties can be assessed: agreement between models and statistical significance of change according to each model. A combination of these two properties was referred to as robustness by Knutti & Sedláček (2013). Their approach holistically considered the magnitude of change, the sign, natural variability and the inter-model spread. Here, we have adopted it, in the same way as in the assessment of robustness of T_min, T_max and P projections in the VOB (Piniewski et al. 2017a). Four categories were classified to indicate the robustness of R projections: (1) high robustness (associated with a high fraction of models with statistically significant changes and high agreement on a given magnitude of change); (2) lack of significance (low fraction of models with statistically significant changes); (3) inconsistent response (high fraction of models with statistically significant change despite low agreement between models on a given change); and (4) neither of the three above categories, i.e., a limited evidence for change of the indicated magnitude. Note that the thresholds chosen to delimit different categories (e.g., the fraction of model outputs without statistically significant change), were selected for this particular paper, as per suggestions by Knutti & Sedláček (2013).

For the purpose of a more holistic view, encapsulating spatial heterogeneity of the studied river basins, the MME median of mean monthly Q projections for the two most downstream gauges of the Vistula (Tczew) and the Odra (Gozdowice) rivers, situated close to the river mouths, were also analysed. Here, the uncertainty was depicted as the monthly Q band, whose lower and upper limits were the 10th and 90th percentiles among the MME, respectively.

Finally, in order to better understand the mechanisms of runoff change, the total runoff response is decomposed into the response of its three main components: surface runoff, lateral flow and baseflow, on both annual and seasonal levels.

The spatial distribution of runoff in the VOB for the historical period, resulting from a MME analysis carried out within this study, is presented in Supplementary material Figure S1 (annual runoff) and Supplementary material Figure S2 (seasonal runoff) (available with the online version of this paper). Spatial variability of mean annual R (Figure S1) is very high, reflecting spatial climatic gradients, such as from the south to the centre (the strongest one), from the north to the centre (intermediate), and from the east to the west (the weakest one), as well as some local-scale fluctuations related to differences in land cover, soil and hydrogeology. Topography is an important driver of spatial distribution of mean annual R. Its high values correlate with the mountains and uplands in the south. While the areal mean R is equal to 166 mm (which constitutes 24% of precipitation), the range depicted by the 1st and 99th percentiles (calculated across 2,633 sub-basins) is 44–563 mm.

Seasonal R dynamics is, in general, driven by two mechanisms: high surface runoff and baseflow as a result of snow melt in spring (areal mean R equal to 60.6 mm), and depleted groundwater contribution (areal mean R equal to 25.6 mm) in autumn, as a result of negative water balance in summer. The importance of these two mechanisms is controlled by local conditions, which results in spatial variability of seasonal R magnitudes. However, spatial trends described above for annual R are generally persistent also in the seasonal R maps (Figure S2).

The maps of projected MME mean changes in the mean annual R show that increases in R dominate across all RCPs and time horizons (Figure 2). The highest mean areal increase of 41.6% is projected under RCP 8.5 in the far future, when both T and P are projected to increase most. However, spatial variability of projected increase is quite pronounced. The north-western corner of the inner part of the VOB has consistently the highest rate of runoff increase, whereas the areas along the southern mountainous edge are characterized by the lowest change, in many cases classified as ‘an agreement on insignificant change’. It should be noted that the areas with the highest change correspond to regions with the lowest R magnitude in the historical period (cf. Figure S1), while the areas with the lowest change correspond to regions with the highest R magnitude. This shows that spatial variability expressed in absolute values (e.g., mm) would be lower than when expressed as per cent change in the maps. We refer to Table 1 for projected mean annual and seasonal runoff values expressed in mm.

Different combinations of RCPs and time horizons differ substantially with respect to the robustness assessment of R change projections (Figure 2). For large parts of the southern half of the VOB, the vast majority of models (at least eight out of nine) agree on the lack of statistical significance in the near future under RCP 4.5 (hatching). In the far future under RCP 4.5, around half of the area is classified as an inconsistent response (white overlay). This corresponds well to the P change map for this combination of future horizon and RCP (Supplementary material, Figure S3, available online). It can be noticed that SWAT projections driven by individual climate model simulations in this case provide quite substantially variable changes, ranging from 14% to 36% (Supplementary material, Figure S4, available online). Spatial patterns in individual maps also differ, although it is quite consistent that in the south the rate of increase is lower, or even decreases occur, over small areas. In the far future, under RCP 8.5, there are small areas in the centre and north, for which projections were assessed as robust (stippling). Areas with overlaid white in this map also correspond to the patterns in P change map for this combination.

It should be noted, though, that robustness spans through the continuum of possibilities, determined by the values of R_KS (robustness in the sense of Knutti & Sedláček 2013) and the fraction of models showing statistically significant change, ν_s. As an example, Figures S5 and S6 of the Supplementary material (available online) show the spatial variability of R_KS and ν_s for annual R projections, illustrating the complexity behind the maps shown in Figure 2.

An aggregate view on mean areal changes in average annual T, P and R is presented in Figure 3. First, it should be noted that all 36 climate model simulations project increases in T and P and all SWAT simulations driven by these climate models project increases in R. Second, precipitation and temperature change are positively correlated (R² = 0.56). Third, there is a very strong positive correlation between runoff change and precipitation change (R² = 91), while runoff change and temperature change are also positively correlated (R² = 0.5).

For all four RCP-horizon combinations, there is a similar seasonal pattern for average VOB conditions: the highest increase projected in winter, the lowest in spring, and the intermediate, similar values in summer and autumn (Figure 4). Winter is the only season in which runoff increases in all sub-basins and all horizon-RCP combinations. In contrast, for spring there is a horizontal divide, relatively stable across RCPs and horizons, north of which the runoff increases, and south of which the runoff decreases. This southern part embraces the areas where, presently, runoff values are higher – the mountain ranges of the Carpathians and Sudetes, as well as their foothills. Summer is another season for which runoff decreases are projected in the southern edge of the VOB. However, it should be noted that in most cases (with the exception of spring in the far future under RCP 8.5), the projected decreases in spring and summer are classified as statistically insignificant. Inconsistent model responses (white areas in the maps) are more frequent in the far future than in the near future, which is particularly visible for summer and autumn.

While the robustness method of Knutti & Sedláček (2013) lends itself well for analysing the uncertainty of climate models in a spatially explicit manner, box plots of areal mean runoff changes calculated across all ensemble members are another way of illustrating the spread of projections (Figure 5). The seasonal patterns of change are similar in all RCP-horizon combinations, with the highest model spread (and magnitude of change) in winter season, the lowest spread (and magnitude) in spring season, and intermediate positions of summer and autumn seasons.

Mean monthly hydrographs in the historical and future climate under both studied RCPs for the two most downstream gauges located close to the mouths of the Vistula and Odra rivers (Tczew and Gozdowice, respectively) show that there is a clear change in hydrograph pattern with changing RCP values and with more distant future periods (Figure 6). The occurrence of maximum monthly Q is steadily shifted from April to March for both rivers, and in the case of the River Odra in the far future under RCP 8.5, February is almost as likely to be the month with flow maximum as March. While in the near future the minimum of the monthly regime is projected to occur in October or November; as for the historical period, in the far future it is projected to be October rather than November, particularly under RCP 8.5. The most robust signals can be observed for all three winter months (December, January and February) in the far future. These are the only cases when the uncertainty bands of the future Q do not overlap at all with those of the historical Q. The highest, approximately two-fold increase is projected for January. In contrast, in the near future the overlap of uncertainty bands is relatively strong, which is consistent with the fact that large parts of the VOB area were characterized by a non-significant change in R for all seasons apart from winter (Figure 4).

Total runoff (R) simulated in SWAT is the sum of surface runoff, lateral flow and baseflow. On an annual basis, they constitute 42%, 13% and 45% of areal mean total runoff, respectively, according to the ensemble mean simulations in the historical period. Projected change in these three components is not uniform, as demonstrated in Figure 7, for the ensemble mean. Relative changes in both baseflow and lateral flow follow a similar seasonal pattern, with the largest increases in winter (more than double in the far future under RCP 8.5) and the lowest in summer. On an annual basis, this increase varies approximately between 20% and 60%. A different behaviour can be observed for surface runoff. Although winter is also the season with the highest increase, surface runoff in spring is projected to decrease by up to 26% for the high-end RCP trajectory in the far future. On an annual basis increases in surface runoff are relatively low and vary between 4% and 15% (ensemble means for four RCP-horizon combinations).

These patterns can be linked to changes in other hydrological processes (snow melt, infiltration) over winter and spring. Both these seasons are characterized by distinctly higher relative increase in precipitation than summer and autumn (cf. Table 1). Additionally, winter is the season with the highest warming, which, in combination with more precipitation, implies considerably more rainfall and thus higher infiltration and shallow aquifer recharge. Thus, rivers are fed with more abundant baseflow and lateral flow, not only during winter and spring, but also, due to a lag effect, during summer, a typical low flow period. Further consequences of climate change in winter and spring seasons are a shorter snow cover period and a shift in snow melt. As shown in Supplementary material Figure S7 (available online), snow melt between November and January is unlikely to undergo a significant change in all scenarios apart from RCP 8.5 in the far future. This suggests that the amount of snow fall in these months will remain relatively stable. A different situation occurs in March and April (February is a ‘transition’ month). Here, snow melt decrease is huge (for example, from 31 mm in March in the historical period to only 5 mm in RCP 8.5 – FF), which triggers some fundamental changes to the hydrology of the VOB. Much less snow melt translates into lower surface runoff in spring. Since snow melt is the main source of flood waters in Polish reference period climatology, advancement of snow melt has a direct consequence on the timing of floods (Figure 6).

The goal of this study was to assess projected changes in runoff in the Vistula and the Odra river basins using a MME of recent, bias-corrected climate model simulations with the help of the SWAT model. The added value of our investigation with respect to the current state-of-the-art, that is dominated by the European-scale studies on one end, and small catchment-scale studies on the other end, is a regional coverage along with a high spatial resolution of all projection maps. Furthermore, this study has quantified the robustness of runoff projections using the approach of Knutti & Sedláček (2013). Our maps point out areas with changes that are not assessed as statistically significant (hatching), as well as areas with a high robustness (stippling) or the lack of robustness (overlaid white colour). In general, at low warming levels the lack of statistical significance of changes dominates, while with increased warming these areas gradually disappear, which is consistent with the findings of Sedláček & Knutti (2014), who analysed the robustness of global projections of the water cycle variables derived from CMIP5 (Coupled Model Intercomparison Project – Phase 5) GCMs. It is noteworthy that spatial patterns of occurrence of different robustness categories in the annual and seasonal R change maps (Figures 2 and 4) follow to a considerable extent the corresponding maps with P change projections (Piniewski et al. 2017a). Nevertheless, some differences are visible, but their explanation would require analysis of other potential factors controlling spatial variability in runoff response, such as land cover, soils or hydrogeology. Such analysis could be done in future studies.

Maps of projected MME mean changes in the mean annual R for the Vistula and Odra river basins, across all nine climate model experiments, for the near and far future, and under two selected RCPs, 4.5 and 8.5, show that increases in runoff are dominating. The seasonal patterns of change are similar in all RCP-horizon combinations, with the highest model spread (and magnitude of change) in winter season and the lowest spread (and magnitude) in spring season. The projections show a clear change in hydrograph pattern in the future: both high spring Q and low autumn Q advance.

Our results suggest, overall, a wetter future for this region, both in terms of P and R. It is well known that according to the EURO-CORDEX projections under RCPs 4.5 and 8.5, the annual total P is expected to increase in central, northern and eastern parts of Europe, an area delimited approximately by the latitude 45°N and the longitude 10°E, and these results have been assessed as significant and robust (Jacob et al. 2014). A similar assessment focused on the VOB using an ensemble of bias-corrected EURO-CORDEX projections (Mezghani et al. 2016) confirmed the high level of agreement on direction of change in P over the VOB, although their robustness and significance were questionable (Piniewski et al. 2017a). SWAT-derived runoff projections for the VOB in the present study show 100% agreement on the direction of change in annual R and a very high correlation between the areal means of annual R and P changes (R² = 0.91, N = 36).

It would be interesting to know how different hydrological models would respond to a similar forcing in this region, but it is hard to find suitable studies for comparison, due to a general scarcity of climate change impact assessments on hydrology driven by EURO-CORDEX. Alfieri et al. (2015) presented a map of mean annual European streamflow change simulated with LISFLOOD under RCP 8.5 for the 2080s, indicating that the majority of rivers in the VOB would undergo a 15–25% increase in streamflow. Papadimitriou et al. (2016) presented a map of change in mean annual runoff in Europe obtained using the JULES land surface model driven by the bias-corrected EURO-CORDEX projections for the +4 °C horizon. Their map showed that in most of the VOB mean annual runoff was projected to increase by up to 25%, but the robustness was high only in the northern half of this region. The last European-scale assessment of hydrological impacts of climate change in Europe with EURO-CORDEX was derived using three models (VIC, LISFLOOD and HYPE) for the +2 °C horizon, but focused on high and low flows rather than mean runoff or streamflow (Roudier et al. 2016). It can be expected that since in their study both low and high flows were projected to increase over the majority of the VOB, mean runoff would undergo the same direction of change. Thus, all available pan-European studies corroborate our findings on the projected runoff increase in the VOB. Additionally, modelling studies for the Upper Western Bug catchment in Ukraine (part of the VOB) using SWIM (Didovets et al. 2017), for three Estonian catchments using SWAT (Tamm et al. 2016) and for nine Polish catchments using HBV (Osuch et al. 2016), all driven by EURO-CORDEX models, also suggest a wetter future. The seasonal patterns of change in runoff as well as the key role of winter and spring precipitation increase influencing baseflow increase throughout the whole year, reported in the study of Didovets et al. (2017), also corroborate our findings. Furthermore, flood projections of Osuch et al. (2016) showed a rapid acceleration of hazard with time, and slightly higher changes for RCP 8.5 than for RCP 4.5, in accordance with results of our study. It should be noted though, that there is a plethora of specific reasons why hydrological projections differ (cf. Piniewski et al. 2017c) so the evidence gathered here cannot be regarded as general. Higher value can be ascribed to evidence gathered from model-intercomparison studies that follow consistent protocols, such as the Inter-Sectoral Impact Model Intercomparison Project (ISI–MIP, www.isimip.org). Hence, it would be useful to investigate the hydrological model choice uncertainty, as yet unexplored in the VOB, in the future climate change impact assessments in this region.

Surprisingly, the rate of R increase was also positively correlated with the rate of T increase, which, in our opinion, can be explained by two reasons: (1) a positive correlation between temperature and precipitation increase in the EURO-CORDEX projections; and (2) a low response in evapotranspiration change in response to climate warming, as simulated by SWAT. While the first aspect is consistent with the fact that globally, about 1.5–2% increase in P per °K surface warming is estimated in response to greenhouse gas forcing (Allen & Ingram 2002), the second one requires more attention. In this study we used the Hargreaves method for PET estimation, which resulted in rather modest changes in PET (areal mean MME mean increase ranging from 2.9% for RCP 4.5 in the near future to 8.8% for RCP 8.5 in the far future). Until recently, the uncertainty of PET and runoff response to climate change due to different PET calculation methods was largely ignored, whereas lately it has attracted more serious attention (Hosseinzadehtalaei et al. 2016; Koedyk & Kingston 2016; Williamson et al. 2016). While the outcomes of these studies vary from region to region, all of them stress the importance of the choice of PET method in estimation of hydrological impacts of climate change. Hosseinzadehtalaei et al. (2016) used the Penman–Monteith FAO 56 (PMF-56) method as a benchmark for comparison with other PET methods in their assessment of future PET projections over Belgium, reporting that Hargreaves estimates of PET increase were on average lower than those derived with PMF-56 by 24–42%, depending on the season. On the other hand, Koedyk & Kingston (2016) showed that even though the choice of PET method influenced PET estimates in their New Zealand study, it hardly affected runoff estimates. Due to a site-specific character of the role of PET for hydrological modelling and the fact that some authors advocate the use of energy-based (PMF-56) rather than temperature-based PET methods, such as Hargreaves (Hosseinzadehtalaei et al. 2016; Williamson et al. 2016), it is recommended that future studies investigate this source of uncertainty in the VOB, which has not been done so far.

Since the VOB covers over 90% of Poland, the significance of this study, i.e., the first almost country-wide assessment of climate change impact on runoff, for the water management planning in Poland is high. The shortage of water resources is of serious concern in Poland and several other Central and Eastern European countries (Kundzewicz 2001), thus the perspective of increasing trend in annual runoff might seem comforting. Yet, climate projections are more complex, as they also show: (1) major seasonal changes (huge winter R increase, earlier occurrence of annual Q maximum); (2) increased flood and reduced drought hazard (Piniewski et al. 2017c); (3) increased nutrient loss hazard, especially for nitrogen (cf. Marcinkowski et al. 2017); (4) impacts on various water users, such as hydropower companies (increased potential, cf. Tamm et al. (2016)) and coal power plant operators (water used for cooling), inland navigation sector, etc.; (5) challenges of meeting the water needs of freshwater (instream and riparian) ecosystems, i.e., the silent water users (cf. Piniewski et al. 2014). In the context of the European Union's Water Framework Directive, it is important that this knowledge mainstreams into Poland's water policy, especially within the second cycle (2015–2021) of the River Basin Management Plans, whose recent updates for the Vistula and the Odra river basins (http://www.apgw.kzgw.gov.pl/) lacked consideration of any projections beyond 2030. Despite the inherent uncertainty, projections for the future are important to inform the process of adaptation to climate change, sketching the range of possible futures. Decision-makers responsible for water management and climate change adaptation have to be aware of the added uncertainty introduced by increased greenhouse forcing (Kundzewicz et al. 2017).

In our analysis for the VOB, runoff increases dominate, regardless of RCP pathway and future horizon. The southern mountainous belt of the VOB is characterized by the lowest magnitude of change, which has been characterized as the lack of statistically significant change in projection maps. The inner, lowland part of the VOB is, in contrast, subject to the highest magnitude of relative change. High correlation between areal mean runoff and precipitation changes, together with similar patterns of projection maps for these variables, suggest that precipitation change, and not evapotranspiration change, is the major driver for runoff change in the VOB. Runoff increase is consistent through all seasons, although variable: very distinct, and even robust in one case (RCP 8.5 in the far future), in winter, the lowest and the least uncertain in spring, and medium, yet with uncertain magnitude, in summer and autumn. These changes are reflected in streamflow hydrographs simulated for the main outlets of the investigated rivers: advancement of maximum flows from April to March or even February, as well as a sharp increase in January/February flow are two most important characteristics. The main mechanism leading to runoff and streamflow increase is increased infiltration and sub-surface flow components (lateral flow and baseflow), whereas increases in surface runoff are low. A huge decrease in snow melt over March and April leads to a decrease in surface runoff in spring.

Assuming that the most recent generation of GCMs and RCMs, with high resolution outputs available, constitutes the most credible projections to date, an overall agreement of our findings with runoff or streamflow projections reported in several other studies, based on various hydrological models forced by EURO-CORDEX data, is encouraging.

Support of the project CHASE-PL (Climate change impact assessment for selected sectors in Poland) of the Polish-Norwegian Research Programme operated by the National Centre for Research and Development (NCBiR) under the Norwegian Financial Mechanism 2009–2014 in the frame of Project Contract No. Pol Nor/200,799/90/2014 is gratefully acknowledged. The first author is grateful for support to the Alexander von Humboldt Foundation and to the Ministry of Science and Higher Education of the Republic of Poland. The Institute of Meteorology and Water Management – State Research Institute (IMGW-PIB) is kindly acknowledged for providing the hydrometeorological data used in this work. The referees are acknowledged for providing useful comments that helped to improve the quality of the manuscript.