1. Introduction
In recent years, several studies have focused on the impact of climate change on the hydrological cycle. Understanding the impacts of future climate on river flood risks and water resources is important to plan effective adaptation strategies to manage the expected changes in extreme event risks [
1]. Flood risks are generally connected to the intensity and frequency of rainfall events in urban areas either located downstream of small catchments or subject to pluvial floods. Several studies have shown an increase pattern in precipitation intensities and the number of extreme events in a warmer climate [
2,
3,
4]. Therefore, flood risks can increase in urban areas in the future [
5,
6], imposing high costs to aquatic and terrestrial ecosystems, human societies, and the economy [
7]. The Intergovernmental Panel on Climate Change (IPCC) [
8] identified statistically significant increasing trends in the number of heavy precipitation events in some regions, claiming that the frequency of extreme precipitation events or the proportion of total rainfall of such events will likely increase in the 21st century over many areas of the globe due to anthropogenic influences that have contributed to the intensification of extreme precipitations at the global scale.
Nevertheless, such increasing frequency and magnitude of extreme rainfall events do not imply that the long-term statistical trends of flood peaks will be also increasing [
9]. The expected impact of climate change on the water cycle can be analyzed by using two approaches. First, trend analyses of observations recorded in the past with a set of statistical models can predict what will happen in the future. Second, rainfall-runoff models that use climate projections as input data can extract the flood change signal predicted by climate models [
10]. Regarding the spatial scales, while the first approach is usually applied to large-scale studies, the second usually regards either large or small and finer spatial resolutions, such as catchment scales.
The input data in the second approach consist in a combination of global climate models (GCMs) and regional climate models (RCMs). GCMs simulate the climate behavior under a set of given representative concentration pathways (RCPs). RCMs downscale such climate projections on a finer spatial scale. Indeed, it is necessary to reproduce rainfall projections with adequate spatial and temporal scales specially to force hydrological models for small and medium catchment areas [
11]. An ensemble of RCMs and GCMs is preferred to a single climate model to remove potential biases in the outputs of RCMs or GCMs [
12]. In climate change studies, the most usual emission scenarios are the RCP 4.5, with the peak of greenhouse gas emissions around 2040 and a decline in the rest of the century, and the RCP 8.5, with increasing greenhouse gas emissions throughout the 21st century [
13].
Many studies have analyzed the expected changes in flood risks in Europe by using both trend analyses and hydrological models as described above [
14]. Recent studies to identify large-scale trends in observed time series have shown that a changing climate in the last decades both increases and decreases European river floods [
15]. The results depend mainly on the part of Europe that is considered. While southern and eastern Europe show negative trends in river floods, the northern part of Europe has increasing trends. Nevertheless, some studies that use climate projections in hydrological models point out that, on average, flood peaks with return periods above 100 years can double their frequency within three decades in Europe [
16]. Other studies confirm the increasing trend for the 100-year flood, particularly in western and eastern Europe, though they found decreasing trends in northeastern Europe [
17,
18]. Several large-scale studies in Europe that use hydrological models differ in their results, as they consider different climate models, downscaling techniques, or bias-correction methods [
19]. For instance, results regarding the future design rainfalls of the Iberian Peninsula obtained from studies with larger scales are conflicting. Indeed, some findings obtained at the European scale show an increase of river flood risks in Spain [
20], though global scale studies show the opposite in the same area [
21]. In addition, the results can be different for a given region even using the same scale. Hence, at the European scale, the results obtained for Spain by Rojas et al. [
18] are in contrast with Alfieri et al. [
16] and Roudier et al. [
20]. This lack of agreement cannot be read easily by decision makers. An additional source of uncertainty that could affect all these studies is that RCMs cannot reproduce extreme rainfall events in the future adequately. Indeed, Herrera et al. [
22] show that RCMs have a good agreement with the observed mean precipitation in Spain, though upper percentiles that represent the amount of total rainfall from extreme events are underestimated. Furthermore, large-scale studies cannot provide flood risk assessment under climate change conditions at the catchment scale, despite that the quantification of flood quantile changes in urban areas would be interesting for municipalities and urban planners. Moreover, large-scales studies are mainly conduct for the 100-year flood, though other return periods are not usually considered.
Consequently, several studies have considered the impact of climate change in smaller river catchments, avoiding the extrapolation of the results from large-scale studies. They evaluate the climate change impact on the water cycle with different objectives, as climate change can have implications in sustainable management of ecosystem services [
23], daily, monthly, and yearly streamflow patterns [
24], seasonal precipitation [
25], or engineering hydraulic design [
10]. Indeed, the use of RCMs in the Mediterranean area at the river basin scale also demonstrated that different parts of southern Europe could have problems related to water scarcity in the future, as precipitation will decrease by the end of the 21st century in the central part of Greece [
25], and a reduction of water yield mainly forced by decreasing precipitations has been detected in the Taro river basin (Italy) [
23]. In Spain, similar results have been found for the time period 2021–2050. In this case, they are mainly due to the increase in both maximum and minimum temperatures at the national scale [
26].
Findings in flood quantile changes in small-scale studies either can differ from the results obtained in the same zone by large-scale studies or can have different change signals for all the considered climate models. For instance, Ducharne et al. [
27] found that low-flow quantiles could decrease, and high-flow quantiles could not significantly change in the Seine and Somme catchments (France), where European scale studies found an increase of design discharge. Statistically significant changes could be neglected when the mean of the results of the different trends in the ensemble of regional climate models is used [
28]. The uncertainty associated with the results is not due only to the downscaling techniques but also on hydrological models. There are several studies that analyze the ensemble of climate models with an ensemble simulation approach [
27,
29,
30,
31].
This study aims to quantify expected flood quantile variations under climate change conditions in Pamplona (Spain) using climate change projections and a distributed rainfall-runoff model. Expected variations in precipitation quantiles extracted from climate change projections in the Iberian Peninsula in a recent study [
32] are used as input data of the Real-time Interactive Basin Simulator (RIBS) model [
33,
34]. The Arga River catchment located in northern Spain is considered as case study, as it is the main driver of the largest floods observed in recent years in the River Ebro [
35]. In addition, Pamplona has suffered frequent fluvial floods in last decades. This study focuses on a single medium-sized basin and on a single distributed hydrological model, trying to minimize the uncertainties associated with the results. Rather than using an ensemble simulation approach, a part of the work is dedicated to conducting a statistical calibration of the RIBS distributed hydrological model to reduce the model uncertainty. Furthermore, this region has been also chosen because it shows low biases associated with the climate models in reproducing the annual maxima series in the control period. This study also aims to analyze the results of the 12 climate models to identify given climate models that always provide the greatest or smallest flood quantile changes in the future.
2. Methodology
The methodology, proposed to identify flood quantile modifications driven by climate change, is divided into five parts: description of the RIBS hydrological model (
Section 2.1), calibration methodology of the RIBS model (
Section 2.2), hydrological modeling to quantify flood peaks in the current scenario (
Section 2.3), spatial distribution of the future design storms (
Section 2.4), and flood quantile change quantification (
Section 2.5).
2.1. RIBS Model
The Real-time Interactive Basin Simulator (RIBS) is a distributed hydrological rainfall–runoff model that is typically used for real-time application on medium-size river basins. RIBS requires the information in a raster format, in this case with a cell size of 50 m. A digital elevation model (DEM) is used to determine the flow direction and flow accumulation in each cell of the domain. The soil information is used to estimate the parameters in the Brooks–Corey equation (Equation (1)) to calculate the part of the rainfall depth that is transformed into runoff in each cell.
where
Ks(y) is the saturated hydraulic conductivity (mm/h);
K0n is the saturated hydraulic conductivity at the soil surface in normal direction (mm/h);
f is he saturated hydraulic conductivity decay in depth (mm
−1),
y is the soil depth (mm);
is the soil moisture content,
is the residual soil moisture content; i.e., the minimum value under which the humidity cannot be extracted by capillary forces;
is the saturated moisture content; and ε is the index of soil porosity [
36]. The Brooks–Corey equation assumes that
Ks(y) has a maximum value at the soil surface (
and then exponentially decrease in the normal direction in the soil depth
with the parameter
.
When rainfall intensity exceeds the infiltration capacity of the soil, surface runoff is generated. Runoff is propagated in the domain through Equations (2) and (3). The catchment is divided into two parts based on a given flow accumulation threshold. Cells with flow accumulation values under such a threshold represent hillslopes in the model and have a runoff velocity equal to
. Cells with flow accumulation values over the threshold represent streams and have a water velocity equal to
. It is assumed that
depends on the discharge in the catchment outlet with respect to a reference flow rate
(Equation (2)).
where
is the water velocity in streams (m/s);
is a model parameter (m/s);
is the discharge at the catchment outlet in time step
t (m
3/s);
is a reference discharge (m
3/s); and
r is a model parameter. If
equals 0,
equals
, remaining constant throughout the simulation. If
is greater than 0,
is greater than the parameter
when the discharge at the catchment outlet is higher than
. At a given time step, the velocity in hillslopes,
, is related to the velocity in streams,
, with the dimensionless parameter
(Equation (3)).
Both velocities can be considered uniform in all the stream and hillslope cells at a given time step, as
depends only on the discharge at the catchment outlet. More information about the RIBS model can be found in Mediero et al. [
37] and Garrote and Bras [
33,
34]. The parameter
in Equation (4) has been set equal to 0 to reduce the computational cost of the model, as the number of simulations with climate change scenarios is high. Therefore, velocity in streams (
) is always equal to
in all the time steps. The model parameters that are considered during the calibration procedure are reduced to three:
(Equation (1)),
(Equation (2)), and
(Equation (3)).
2.2. Calibration of the RIBS Model
The calibration procedure aims to find the model parameter values that minimize the errors between the model results and the observations of the real system that it is intended to represent [
38]. Such errors in model simulations can be represented with Equation (4) [
39].
where
is the measured discharge in a given location
at the time step
;
is the modeled discharge that is obtained using a set of parameters
; and
is the error in the same location and time step with a set of parameters
.
The model calibration can be based on optimization methods, seeking to minimize model errors. Model errors can be quantified by using different objective functions. A set of objective functions are needed to consider a set of aspects in the hydrograph, such as the hydrograph shape, the timing and magnitude of the peak and the magnitude of errors in low flows, among others.
In this study, five objective functions are used to calibrate the RIBS model. Three of them quantify the errors considering the complete hydrograph duration: mean absolute error (
MAE), root mean square error (
RMSE), and the Nash–Sutcliffe efficiency coefficient (
NSE) (Equations (5)–(7)). Two additional parameters are considered to assess the fitting in the upper part of the hydrograph: time to peak (
TP) and the magnitude of the peak (
MP) (Equations (8) and (9)).
where
is the measured discharge at a given time step
;
is the simulated discharge with a set of model parameters
at the time step
; and
is the mean value of the observed discharges.
The
NSE objective function [
40] quantifies the improvement of using the simulated hydrograph compared with the average discharge of the observed hydrograph. The best value of all the objective functions is 0, except for
NSE that supplies a value equal to one in a perfect match between the simulated and measured hydrographs. Indeed,
MAE and
RMSE values equal 0 and
NSE values equal to 1 mean that there are not biases between the measured and the simulated discharge. A
TP value equal to 0 means that the peak discharge of the simulated hydrograph occurs in the same time step as in the observed hydrograph. A
MP value equal to 0 means that there is no difference between the magnitude of simulated and observed peaks.
2.3. Current Scenario
The design floods in the current scenario have been obtained by using the RIBS hydrological model with a set of design hyetographs as input data, obtained through an extreme frequency analysis. Daily observations at eight rain gauge stations have been considered. Both the GEV and the Gumbel distribution functions have been used to fit a frequency curve to annual maximum series. The return periods of 2, 5, 10, 50, 100, 500, and 1000 years are considered. An Areal Reduction Factor (ARF) is considered to reduce the mean areal precipitation in the catchment obtained from the precipitation quantiles estimated at rain gauge locations, as annual maximum rainfalls in the eight rain gauge stations are unlikely to occur on the same day in all the years.
The mean rainfall intensity for a given subdaily duration is obtained from the intensity–duration–frequency curve in the region. A storm duration of 24 h is considered. The 24 h design rainfall for a given return period is obtained by scaling a dimensionless hyetograph by the T-year daily rainfall multiplied by ARF. Therefore, the design hyetographs in the Arga River catchment have a fixed shape with the peak intensity in the central part of the event regardless the return period.
Finally, at-site rainfall precipitations at each time step given by the design hyetographs in the eight gauging sites are spatially distributed in the catchment, by using a grid with the cell size used by the RIBS model. The Thiessen polygon technique is used to obtain the precipitation fields in each time step. Therefore, 24 precipitation fields in the catchment are obtained for each design event, one for each time step of the hyetograph.
2.4. Climate Change Scenario
Expected changes in daily precipitation quantiles in the future under climate change conditions were extracted from rainfall climate projections supplied by 12 combinations of GCMs and RCMs (
Table 1) of the EURO-CORDEX program [
32]. Delta changes in daily precipitation quantiles are supplied for seven return periods (RP = 2, 5, 10, 50, 100, 500, and 1000 years), two representative concentration pathways (RCP 4.5 and RCP 8.5), and three time windows in the future (2011–2040, 2041–2070, and 2071–2100). A delta change represents the expected variation in rainfall quantiles for each combination of return period, climate model, time window, and RCP scenario. The delta changes are provided in a grid with a cell resolution of 0.11°. Hence, the Arga River catchment is covered mainly by three points that have the greatest influence on the design rainfall variation, though additional seven points are also considered, as they cover smaller parts of the river catchment (
Figure 1). Therefore, 504 values of delta changes are considered in each cell that are the possible combinations with the 7 return periods, 2 emission scenarios, 3 time windows, and 12 climate models considered in the study.
The delta changes are spatially distributed with the 50 m grid cell size that is used in the RIBS hydrological model simulations. Though the best geostatistical method depends on each case characteristics, ordinary kriging has been proved to provide better results than the inverse distance weighting (IDW) method in several studies regarding precipitation fields [
41]. Nevertheless, in this case study the IDW technique has been used with a high exponent to maintain the same original squared shape of the delta changes. Furthermore, the use of the kriging method is not possible with these data as a fitting model cannot be determined for the semivariogram. Indeed, all the points have a fixed reciprocal distance. Therefore, for the same distance in the x-axis, there are several semivariance values in the y-axis that prevent any model from being able to converge (gaussian, exponential, and spherical).
The spatial distributions of the future design precipitations are obtained combining the spatial distributions of both the current design rainfalls and the delta changes.
Figure 1 shows the spatial distribution of the peak intensity in the 2 year hyetograph and the delta changes for the CNR-CCL climate model for the same return period in the 2040 time window. Considering the 504 combinations in each cell described above and the 24 time steps of the design hyetographs, a total of 12.096 rainfall spatial distributions have been considered as input data in the RIBS model.
2.5. Quantification of Expected Changes in Flood Quantiles
The expected changes in flood quantiles in the future are obtained from the comparison between flood quantiles in the future periods and the current period for each climate model. The flood quantile change supplied by a given climate model and return period (Δ
QT) are obtained using the ratio between the hydrograph peak discharge in the climate change scenario
and the hydrograph peak discharge in the current scenario
(Equation (10)). Δ
QT quantifies the change between a future time window and the current period, extracting the signal change supplied by a climate model forced with a given emission scenario, regardless the flood frequency curve obtained from observations. Therefore, the potential biases between the model simulations and observations do not have influence in the results.
where (
Qp,T)fut is the hydrograph peak magnitude for the return period
T in the future period;
(Qp,T)curr is the hydrograph peak magnitude for the return period
T in the current period; and Δ
QT is the expected change in the flood quantile for the return period
T.
If ΔQT is greater than one, flood quantiles are expected to increase in the future. If ΔQT is smaller than one, flood quantiles are expected to decrease in the future. A value of one means that no changes in flood quantiles are expected in the future.
3. Data and Case Study
Pamplona is crossed by the Arga River that has a catchment area of 510 km
2. The Arga River is in the Navarre region and rises in the Urquiaga pass, located in the Paleozoic Quinto Real massif, one of the rainiest areas in Northern Spain. It is a tributary of Aragón River that has a total drainage area of 2759 km². In this study, the catchment outlet is located in the Pamplona city, where the A323 streamflow gauging station is placed, as shown in
Figure 2.
Two types of rain-gauge stations are used in the study. For the RIBS model calibration procedure, six 15 min rain-gauge stations are used (purple circles in
Figure 2). For estimating the design rainfall hyetographs, daily rain-gauge stations with longer time series are used (green triangles in
Figure 2). Therefore, a time step of 15 min has been considered in the RIBS model simulations for both the calibration and the calculation with current and future design rainfalls.
Table 2 contains a summary about the streamflow- and rain-gauge station network.
The soil data were supplied by the Spanish National Geographic Institute IGN (Madrid, Spain), “Instituto Geográfico Nacional” in Spanish, in a shapefile format. The DEM was also supplied by IGN in a raster format. All the data in the RIBS model are provided in a raster format with a cell size of 50 m. Therefore, the soil raster in
Figure 3 was obtained sampling the shapefile with a 50 m grid.
Table 3 reports the Brooks–Corey parameter values (see
Section 2.1) for each soil class.
The greatest flood events of the Arga River in Pamplona in the last decade have been selected for the calibration of the RIBS hydrological model. A peak-over-threshold (POT) analysis has been applied to the streamflow time series recorded at the A323 hydrometer located at the catchment outlet. A threshold of 300 m
3/s has been considered, as floods that exceed such a threshold usually drive significant direct losses in the Pamplona city (
Figure 3). Eight flood events have been selected for the calibration and validation methodology of the RIBS hydrological model. They are named from EPI01 to EPI08 in chronological order.
However, EPI01 has been removed from the calibration process, as two rain-gauge stations were out of service for this event. Hence, five events were used in the calibration procedure (EPI02, EPI03, EPI04, EPI05, and EPI06) and two in the validation (EPI07 and EPI08). A summary of the selected flood events is included in
Table 4.
5. Discussion and Conclusions
An analysis of the expected changes in flood quantiles in the future for Arga River in the city of Pamplona has been carried out using precipitation projections supplied by 12 climate models of the EURO-CORDEX program. Delta changes in daily rainfall quantiles obtained in a previous study have been used as input data in the RIBS distributed hydrological model to transform the expected changes in precipitation quantiles into changes in flood quantiles.
For this reason, the RIBS hydrological model has been calibrated with the seven greatest flood events in the city of Pamplona in the last decade. The results of calibration and validation show that the simulated hydrographs have a good fit with the observations, obtaining acceptable residual model errors both in terms of hydrograph shape (RMSE, MAE, NSE) and flood peak (TP, MP).
The results point to a decrease of the design peak discharges for both RCPs in the three time windows for return periods below 10 years. On the contrary, flood quantiles for return periods greater than 50 years are expected to increase for the three time windows in the RCP 8.5 scenario. In addition, the 100-year flood quantile will not change in the future if environmental policies cut greenhouse gasses emission (RCP 4.5). On the contrary, design peak discharges above 50 years are expected to increase 10–40% in the emission scenario RCP 8.5, especially for the time windows 2041–2070 and 2071–2100. GCMs seem to have more influence on the results than RCMs. The greatest increases are predicted by the results obtained with climate models that downscale from the CNRM-CM5 and IPSL-CM5A-MR GCMs.
The results of this work agree with Alfieri et al. [
16] that show an increase of 100-year floods across Europe in the RCP 8.5 on average. The variation of flood quantiles in the Arga River could also be a signal for northern Spain, where a downtrend of the annual maximum flows has been found by Mediero et al. [
43] and Bloesch et al. [
15]. All these findings could be summarized together pointing that without a greenhouse gas emission abatement a generalized increase of the hydraulic risk for greater return periods could be expected also in the areas where the mean annual precipitation is decreasing.
Garijo and Mediero (2018) [
35] found a slight decreasing trend in flood quantiles in Pamplona for RCP 4.5 and a slight increasing trend in flood quantiles for RCP 8.5 with the EURO-CORDEX climate models. In this study, clear increasing trends were obtained mainly for the RCP 8.5. The differences between both studies can be attributed to differences in the methodological approach. Garijo and Mediero (2018) used the continuous HBV model with climate projections using a daily scale, while in this study a distributed hydrological model was used with subdaily precipitations estimated with an approach based on design hyetographs and delta changes obtained from climate projections. Furthermore, the spatial distribution of the future design storms has been considered. Indeed, flood quantile delta changes have not been directly applied to the nearest rain gauges, though they have been interpolated spatially and then combined with design rainfall fields for the current scenario. This methodological approach generated 12,096 spatial distribution of rainfall for the future design storms that were used to force the RIBS hydrological model.
Comparing the results of this paper with similar studies carried out in other parts of Europe, they agree with Hennegriff et al. [
30] who obtained a 15% increase in 100-year flood in Bavaria (Germany), though differing the results about the most frequent floods, as they found an increase up to 75% in the two-year floods. The results also agree with Hattermann et al. [
31] who found a stronger significant increasing trend in flood frequency quantiles, reaching up to 60% increases in the 50-year floods. In addition, Sharma et al. [
44] found similar increases in the most extreme floods, though analyzing climatic areas completely different from the Mediterranean area.
The most interesting result obtained in this study consist in the clear relationship between the flood hazards expected in the future the city of Pamplona and the temporal evolution of emission scenarios. Indeed, the greatest expected changes in the flood quantiles coincide with the peak of the greenhouse gasses emission in both emission scenarios: 2040 for the RCP 4.5 and in 2100 for the RCP 8.5. Therefore, results show how policies that aim to reduce greenhouse gases emissions in the future could lead to a reduction in the future flood risks. This could reduce costs related to flood damages of more extreme events and costs related to the oversized defense infrastructures.
However, the results and considerations of this study cannot be easily extrapolated to other geographical contexts. Nonetheless, they show the importance of the future hydraulic risk assessment for small and medium river basins. Indeed, small river catchments could have differing results from those obtained in large-scale studies that can only analyses the largest river basins in which they are contained.
The results also underline how important it would be to consider the increase of the future hydraulic risks due to climate change as soon as possible, especially for higher return periods and emission scenarios in which environmental policies do not aim to reduce emissions. Considering the expected increases in hydraulic risks for the 1000-year return period, decision makers could design structures such as dams avoiding spillway underestimation in the future. Similarly, the use of flood risk maps with underestimated 100-year floods could also underestimate the areas prone to flooding in urban cities.