Selection of Bias Correction Methods to Assess the Impact of Climate Change on Flood Frequency Curves

Climate projections provided by EURO-CORDEX predict changes in annual maximum series of daily rainfall in the future in some areas of Spain because of climate change. Precipitation and temperature projections supplied by climate models do not usually fit exactly the statistical properties of the observed time series in the control period. Bias correction methods are used to reduce such errors. This paper seeks to find the most adequate bias correction techniques for temperature and precipitation projections that minimizes the errors between observations and climate model simulations in the control period. Errors in flood quantiles are considered to identify the best bias correction techniques, as flood quantiles are used for hydraulic infrastructure design and safety assessment. In addition, this study aims to understand how the expected changes in precipitation extremes and temperature will affect the catchment response in flood events in the future. Hydrological modelling is required to characterize rainfall-runoff processes adequately in a changing climate, in order to estimate flood changes expected in the future. Four catchments located in the central-western part of Spain have been selected as case studies. The HBV hydrological model has been calibrated in the four catchments by using the observed precipitation, temperature and streamflow data available on a daily scale. Rainfall has been identified as the most significant input to the model, in terms of its influence on flood response. The quantile mapping polynomial correction has been found to be the best bias correction method for precipitation. A general reduction in flood quantiles is expected in the future, smoothing the increases identified in precipitation quantiles by the reduction of soil moisture content in catchments, due to the expected increase in temperature and decrease in mean annual precipitations.


Introduction
Climate change projections predict reductions in the availability of water resources, as well as changes in the frequency and severity of extreme hydrological events [1]. Consequently, climate change will affect floods, which are the natural hazard that generates the largest damages in Europe [2]. The impact of climate change on the hydrological cycle can be assessed by using both climatic and hydrological models. Several studies and models have been developed. However, most of such studies are focused on either monthly or annual scales [3]. Only a few studies have studied the impact of climate change on floods [4].
Global climate models (GCMs) are large-scale models that simulate the behavior of atmospheric circulation patterns with a gross resolution of approximately 100-250 km. Consequently, downscaling techniques are required to obtained results at finer resolutions of around 25-50 km [5]. Regional climate models (RCM) are the most usual downscaling technique.

Data and Methodology
The methodology consists of the following steps: (i) Calibration of the HBV model to adequately characterize rainfall-runoff processes in a changing climate; (ii) selection of the best bias correction method for precipitation and temperature projections, and (iii) assessment of the expected changes in flood quantiles in the future driven by climate change.

Study Area and Data
Four catchments in the Douro river basin in the northwestern part of Spain have been selected as case studies ( Figure 1 and Table 1). A dam is located at the outlet of each catchment. Consequently, inflow discharges in reservoirs were estimated from observations of mean daily reservoir water levels and dam releases, which were supplied by the Centre for Hydrographic studies of CEDEX.  Observations of daily rainfall at 42 rain-gauging stations and temperature at 22 thermometric stations were supplied by the AEMET. Gaps in time series were filled by using observations at nearby gauging stations, by using the inverse distance weighting (IDW) method with squared distances.
Climate change projections provided by 12 regional climate models of the EURO-CORDEX programme have been used (Table 2). Such projections are composed of daily rainfall and temperature time series in a grid with cells of 0.11°. The same control period , hydrological years) and future period under climate change (2011-2094, hydrological years) have been considered for all the climate models. Two representative concentration pathways (RCP) were considered: RCP 4.5 and 8.5.  Observations of daily rainfall at 42 rain-gauging stations and temperature at 22 thermometric stations were supplied by the AEMET. Gaps in time series were filled by using observations at nearby gauging stations, by using the inverse distance weighting (IDW) method with squared distances.
Climate change projections provided by 12 regional climate models of the EURO-CORDEX programme have been used (Table 2). Such projections are composed of daily rainfall and temperature time series in a grid with cells of 0.11 • . The same control period , hydrological years) and future period under climate change (2011-2094, hydrological years) have been considered for all the climate models. Two representative concentration pathways (RCP) were considered: RCP 4.5 and 8.5. Table 2. EURO-CORDEX climate models used in the study.

Acronym
CGM RCM

Generalized Extreme Value (GEV) Distribution
Flood quantiles are calculated by fitting a Generalized Extreme Value (GEV) distribution to the annual maximum series (Equation (1)). The GEV distribution with the L-Moment method has been selected following the national recommendations in the region where the catchments are located [23,24]. In addition, the GEV distribution is applied widely in earth system sciences and hydrology to study extremes of several natural phenomena, including rainfall, floods, wind speeds and wave heights [25].
The cumulative distribution function F(x) of the GEV distribution is expressed with the following equation: where u is the location parameter, α is the scale parameter, and k is the shape parameter. The GEV distribution function allows us to obtain the peak flow and precipitation quantiles for high return periods needed for the study.

HBV Model and Calibration
The hydrological response in the four catchments has been simulated with the continuous HBV rainfall-runoff model. Specifically, the HVB-light-GUI 4.0.0.7 version has been used.
The HBV model is a semi-distributed model, dividing a given catchment into elevation and vegetation zones, as well as into sub-catchments. The model supplies simulations of catchment discharge on a daily time step. The model consists of three routines: Snow, soil and groundwater. Each routine provides a flow time series depending on a set of parameters. The snow accumulation and snowmelt are computed by a degree-day method in the snow routine. In the soil routine, groundwater recharge and actual evapotranspiration are simulated as functions of actual water storage. In the groundwater routine, runoff is computed as a function of water storage. Finally, in the routing routine, a triangular weighting function is used to simulate the routing of the runoff until the catchment outlet [26]. Parameters values related to the soil routine have been estimated initially using land-use maps supplied by the CORINE land cover.
The model parameters have been calibrated using Monte Carlo simulations and GAP optimization. Both tools are integrated into the HBV software.
Model parameters have been calibrated using the coefficient of efficiency (Reff) goodness-of-fit function that is integrated into the HBV software. 'Reff' compares the prediction supplied by the HBV Water 2019, 11, 2266 5 of 16 model with the simplest possible prediction, which is a constant value equal to the mean value of observations over the entire period (Equation (2)): where Q sim (t) is the simulated discharge at time step t, Q obs (t) is the observed discharge at time step t and Q obs is the mean value of the observed discharges. A sensitivity analysis with 1,000,000 Monte Carlo simulations was conducted to identify the most important parameters in the four catchments. Maximum storage capacity in the soil (FC), maximum flow from upper to lower box (PERC), upper storage character coefficient 0 (K0) and upper storage recession coefficient 1 (K1) were found to be the characteri that had the most influence on the results. FC, PERC, K0 and K1 are associated with soil infiltration processes that are characterized with three boxes [27]. As expected, the snow routine is not important in the case studies.
Flood quantiles have been estimated by fitting a GEV distribution to the annual maximum series of streamflow observed and simulated by the HBV model. Flood quantiles obtained by simulation have been compared with flood quantiles obtained from observed data, for a set of return periods ( Figure 2). An iteration process has been used, prioritizing the similarity in the simulation of high return period quantiles, as the results of the study will be applied to dam design.
where Qsim(t) is the simulated discharge at time step t, Qobs(t) is the observed discharge at time step t and is the mean value of the observed discharges. A sensitivity analysis with 1,000,000 Monte Carlo simulations was conducted to identify the most important parameters in the four catchments. Maximum storage capacity in the soil (FC), maximum flow from upper to lower box (PERC), upper storage character coefficient 0 (K0) and upper storage recession coefficient 1 (K1) were found to be the characteri that had the most influence on the results. FC, PERC, K0 and K1 are associated with soil infiltration processes that are characterized with three boxes [27]. As expected, the snow routine is not important in the case studies.
Flood quantiles have been estimated by fitting a GEV distribution to the annual maximum series of streamflow observed and simulated by the HBV model. Flood quantiles obtained by simulation have been compared with flood quantiles obtained from observed data, for a set of return periods ( Figure 2). An iteration process has been used, prioritizing the similarity in the simulation of high return period quantiles, as the results of the study will be applied to dam design.

Bias Correction
Statistics of precipitation and temperature projections supplied by climate models in the control period do not usually fit exactly those of observations in the same period. Such errors may affect simulated flow results in the future period. Bias correction methods try to improve the fitting of climate model simulations to observations in the control period, in order to enhance the reliability of climate model results in the future period. First, temperature and precipitation series have been corrected separately with a set of methodologies. Second, temperature and precipitation bias correction techniques are combined to identify the best methodology.

Bias Correction
Statistics of precipitation and temperature projections supplied by climate models in the control period do not usually fit exactly those of observations in the same period. Such errors may affect simulated flow results in the future period. Bias correction methods try to improve the fitting of climate model simulations to observations in the control period, in order to enhance the reliability Water 2019, 11, 2266 6 of 16 of climate model results in the future period. First, temperature and precipitation series have been corrected separately with a set of methodologies. Second, temperature and precipitation bias correction techniques are combined to identify the best methodology.

Temperature Correction
Temperature projections have been corrected by using simple seasonal bias correction [6]. Mean monthly temperatures are obtained in each climate model in the control period for its comparison with observations in the same period. The difference between climate model simulations and observations is added to the projections in the future period for the two emission scenarios (Equation (3)): where T i,j is the mean temperature in month j for the climate model ©, ∆T is the difference between the mean temperature of the climate model i and the observations in month j, and T i,j,corr is the corrected temperature in month j for the climate model i.

Precipitation Correction
There are several methodologies for correcting bias in precipitation projections [9]. Some methods are simple, consisting of obtaining a rainfall threshold to correct rainfall from 0.1 mm to dry days. Other methods are more complex and require a separate precipitation model. Finally, some methods correct the precipitation based on empirical techniques [28][29][30].
This study is focused on dam safety. Consequently, extreme events are the variable of interest. The following bias correction techniques for precipitation projections are proposed. First, quantile mapping (QM) linear correction (Equation (4)). Second, QM second-order polynomial correction (Equation (5)) [30,31]. Both methods require to estimate a set of parameters that are adjusted by comparing the observed data with each climate model simulations in the control period. With such parameter values, precipitation projections are corrected in the future periods for the two emission scenarios.
where P i,j is the raw precipitation supplied by the climate model i in day j, P i,j,corr is the precipitation corrected for the climate model i in day j, and a, b, c, d and e are quantile mapping parameters.

Results and Discussion
The best bias correction method has been identified for: (i) Temperature projections, in terms of monthly averages; (ii) extreme precipitation in climate projections, in terms of frequency curves for annual maximum series; and (iii) extreme simulated discharges.

Bias Correction in Temperature Series
In the four catchments, monthly mean temperatures supplied by climate models are significantly lower than observations (

Bias Correction in Precipitation Series
Climate models usually supply differing precipitation magnitudes in the control period compared to observations (Figures 4-7). In the Barrios de Luna catchment, climate models supply larger extreme precipitations than observations. However, in the other three catchments, climate models supply lower precipitations.
As mentioned previously, in the Barrios, the Luna catchment, the observed precipitation data are smaller than climate model projections ( Figure 3). Both bias correction techniques reduce errors. However, the QM polynomial correction technique is the best method for this case.

Bias Correction in Precipitation Series
Climate models usually supply differing precipitation magnitudes in the control period compared to observations (Figures 4-7). In the Barrios de Luna catchment, climate models supply larger extreme precipitations than observations. However, in the other three catchments, climate models supply lower precipitations.

Bias Correction in Precipitation Series
Climate models usually supply differing precipitation magnitudes in the control period compared to observations (Figures 4-7). In the Barrios de Luna catchment, climate models supply larger extreme precipitations than observations. However, in the other three catchments, climate models supply lower precipitations.
As mentioned previously, in the Barrios, the Luna catchment, the observed precipitation data are smaller than climate model projections ( Figure 3). Both bias correction techniques reduce errors. However, the QM polynomial correction technique is the best method for this case.  In Camporredondo, Porma and Riaño catchments, frequency curves fitted to observations are slightly larger than the median of the climate model projections (Figures 4-6). The best results are obtained for the QM linear bias correction technique, as the QM polynomial correction leads to precipitation quantiles larger than observations, especially for high return periods.    In Camporredondo, Porma and Riaño catchments, frequency curves fitted to observations are slightly larger than the median of the climate model projections (Figures 4-6). The best results are obtained for the QM linear bias correction technique, as the QM polynomial correction leads to precipitation quantiles larger than observations, especially for high return periods.    slightly larger than the median of the climate model projections (Figures 4-6). The best results are obtained for the QM linear bias correction technique, as the QM polynomial correction leads to precipitation quantiles larger than observations, especially for high return periods.    As mentioned previously, in the Barrios, the Luna catchment, the observed precipitation data are smaller than climate model projections (Figure 3). Both bias correction techniques reduce errors. However, the QM polynomial correction technique is the best method for this case.
In Camporredondo, Porma and Riaño catchments, frequency curves fitted to observations are slightly larger than the median of the climate model projections (Figures 4-6). The best results are obtained for the QM linear bias correction technique, as the QM polynomial correction leads to precipitation quantiles larger than observations, especially for high return periods.
Summarizing, bias correction by using both the lineal and polynomial techniques reduce errors compared to raw precipitation data in the four catchments, in terms of precipitation frequency curves. The QM linear correction technique is the best bias correction method in three catchments. The QM polynomial correction technique is the best bias correction method in the Barrios de Luna catchment.

Bias Correction in Terms of Flood Frequency Curves in the Control Period
Simulations of the continuous HBV model have been conducted with a set of combinations of raw and bias corrected temperature and precipitation time series in the control period as input data, in order to compare the bias correction techniques. The best combination of bias correction techniques is identified in terms of the smallest errors in the flood frequency curve estimated with observations. The higher return periods have been considered, due to its importance in dam design.
In the Barrios de Luna catchment, the best fit to the observed data is obtained with raw precipitation and temperature projections with no bias correction ( Figure 8). Bias correction in temperature time series leads to similar results, but with larger errors for high return periods, which are used in dam safety. QM precipitation corrections improve the fitting for low return periods, but increase errors significantly for high return periods. Finally, similar results are obtained for combinations of bias correction techniques in both precipitation and temperature series. Summarizing, bias correction by using both the lineal and polynomial techniques reduce errors compared to raw precipitation data in the four catchments, in terms of precipitation frequency curves. The QM linear correction technique is the best bias correction method in three catchments. The QM polynomial correction technique is the best bias correction method in the Barrios de Luna catchment.

Bias Correction in Terms of Flood Frequency Curves in the Control Period
Simulations of the continuous HBV model have been conducted with a set of combinations of raw and bias corrected temperature and precipitation time series in the control period as input data, in order to compare the bias correction techniques. The best combination of bias correction techniques is identified in terms of the smallest errors in the flood frequency curve estimated with observations. The higher return periods have been considered, due to its importance in dam design.
In the Barrios de Luna catchment, the best fit to the observed data is obtained with raw precipitation and temperature projections with no bias correction (Figure 8). Bias correction in temperature time series leads to similar results, but with larger errors for high return periods, which are used in dam safety. QM precipitation corrections improve the fitting for low return periods, but increase errors significantly for high return periods. Finally, similar results are obtained for combinations of bias correction techniques in both precipitation and temperature series.  In the Camporredondo catchment, a significant improvement was found by using bias correction techniques in precipitation and temperature data (Figure 9). The polynomial correction in precipitation and temperature series led to the smallest errors. In the Camporredondo catchment, a significant improvement was found by using bias correction techniques in precipitation and temperature data (Figure 9). The polynomial correction in precipitation and temperature series led to the smallest errors. The Porma catchment shows larger errors, as the observed data are significantly larger than simulations provided by the HBV model ( Figure 10). However, the smallest errors are obtained for temperature bias corrected and precipitation with QM polynomial correction, mainly for high return periods, where the smallest errors are found. The Porma catchment shows larger errors, as the observed data are significantly larger than simulations provided by the HBV model ( Figure 10). However, the smallest errors are obtained for temperature bias corrected and precipitation with QM polynomial correction, mainly for high return periods, where the smallest errors are found.
In the Riaño catchment, the temperature correction leads to small changes in flood frequency curves. In this case, the smallest errors are obtained for precipitation data with QM polynomial correction ( Figure 11). The technique that leads to the best fit for high return periods is the QM polynomial correction in precipitation data with the raw temperature time series.
In general, the smallest errors are obtained with the polynomial bias correction technique. More specifically, the methodologies with the smallest absolute error for higher return periods are: Raw temperature and raw precipitation supplied by climate models in Barrios de Luna; QM lineal correction for precipitation and mean monthly temperature correction in Camporredondo; QM polynomial correction for precipitation and mean monthly temperature correction in Porma; and QM polynomial correction for precipitation and raw temperature in Riaño.
Finally, the errors of the selected methods in each catchment are summarized in Table 3. In general, smaller errors are obtained for higher return periods, that are usually used for dam design. In the Riaño catchment, the temperature correction leads to small changes in flood frequency curves. In this case, the smallest errors are obtained for precipitation data with QM polynomial correction ( Figure 11). The technique that leads to the best fit for high return periods is the QM polynomial correction in precipitation data with the raw temperature time series.

Expected Changes in Flood Frequency Curves in the Future Period
Precipitation and temperature projections in the future period (2011-2095) have been corrected with the best bias correction technique identified in each catchment in the previous step. Such corrected precipitation and temperature time series are used as input data in the HVB model. Annual maximum series of streamflow series simulated by the HBV model is extracted, and flood frequency curves are estimated by using a GEV distribution function. Delta changes in peak flow quantiles expected in the future are obtained by comparing the flood frequency curves in the future and control periods ( Figure 12).
In the Barrios de Luna catchment, significant differing results are found between RCP 4.5 and RCP 8.5. In In general, the smallest errors are obtained with the polynomial bias correction technique. More specifically, the methodologies with the smallest absolute error for higher return periods are: Raw temperature and raw precipitation supplied by climate models in Barrios de Luna; QM lineal correction for precipitation and mean monthly temperature correction in Camporredondo; QM polynomial correction for precipitation and mean monthly temperature correction in Porma; and QM polynomial correction for precipitation and raw temperature in Riaño.
Finally, the errors of the selected methods in each catchment are summarized in Table 3. In general, smaller errors are obtained for higher return periods, that are usually used for dam design. Figure 11. Comparison between flood frequency curves fitted to annual maximum series simulated by the HBV model with observed data and HBV simulations with climate projections in Riaño catchment, using a set of combinations of bias correction techniques. Red lines represent the median of the 12 climate models considered. (a) Raw data with no bias correction; (b) raw temperature and precipitation with QM linear correction; (c) raw temperature and precipitation with QM polynomial correction; (d) raw precipitation and temperature bias corrected; (e) temperature bias corrected and precipitation with QM linear correction; (f) temperature bias corrected and precipitation with QM polynomial correction. Table 3. Absolute and relative errors of HBV simulations of peak flow quantiles in the control period for the best bias correction method identified in each catchment. similar to the control period in 2011-2040. In the RCP 8.5, flood quantiles will be larger in the period 2071-2100, smaller in the period 2011-2040 and similar to the control period in 2041-2070. In the Porma catchment, similar results are found for all periods and emission scenarios. Flood quantiles will decrease for high return periods, though they will be slightly higher for return periods smaller than 25 years.

Barrios de Luna
Finally, in the Riaño catchment, the results are similar to the Barrios de Luna catchment. In the RCP 4.5, flood quantiles will be larger in the period 2041-2070, smaller in the period 2071-2100 and similar to the control period in 2011-2040. In the RCP 8.5, flood quantiles will be larger in the period 2071-2100, smaller in the period 2011-2040 and similar to the control period in 2041-2070.
Summarising, simulations with the HBV model by using precipitation and temperature projections corrected with the best bias correction techniques identified previously show that, in general, flood frequency curves will decrease in the future, though increases can be seen in some cases.

Conclusions
The results of this study confirm the need for bias correction in temperature and precipitation projections to improve their fitting to observations. Mean monthly temperatures supplied by climate models in the control period are significantly lower than observed data. For precipitation projections, annual maximum series are slightly smaller than observations in three of the four catchments considered. In terms of precipitation frequency curves, bias correction reduces errors compared to raw data in the four catchments. The quantile mapping linear correction technique is the best bias correction method in three of the four catchments. It confirms the results of previous studies that found that the quantile mapping technique is the best approach to correct extreme values of precipitation.
A set of combinations of bias correction methods for precipitation and temperature were analyzed, by comparing the outputs of the HBV model with the observations in the control period. Flood frequency curves fitted to simulations and observations were compared. The best bias correction method for precipitation projections, in terms of precipitation frequency curves, differs from the best method in terms of flood frequency curves. It was found that bias correction of precipitation time series is more important than temperature correction, as precipitation has a larger influence on streamflow outputs of the HBV model. This finding agrees with previous studies that highlighted the importance of correcting biases in precipitation projection with advanced techniques, while correcting only mean values in temperature projections. In terms of flood frequency curves, the best bias correction technique for precipitation is the quantile mapping polynomial correction. Both raw temperature projections and bias corrected temperature time series lead to similar results in terms of flood frequency curves.
Simulations with the HBV model in the future period, using temperature and precipitation projections corrected with the best combination of bias correction techniques identified in each catchment, show a general reduction in flood quantiles, smoothing the increases identified in precipitation quantiles. For example, in the control period, when precipitation quantiles are larger than observations, flood quantiles are similar to observations. In general, the period 2071-2095 presents the smallest reductions and, in some cases, the largest increases. This finding agrees with previous studies that stated that floods could decrease in many parts of southern Europe under climate change conditions. Author Contributions: E.S. wrote the manuscript, conducted the analyses and was the person in charge of the review process. C.G. supplied the climate projection data in the case studies that was processed previously. L.M. coordinated the analyses conducted in the paper and edited the manuscript.

Funding:
The authors acknowledge that this study was supported by the project CGL2014-52570-R 'Impact of climate change on the bivariate flood frequency curve' of the Spanish Ministry of Economy and Competitiveness.