Return Level Estimation of Extreme Rainfall over the Iberian Peninsula: Comparison of Methods

Different ways to estimate future return levels (RLs) for extreme rainfall, based on extreme value theory (EVT), are described and applied to the Iberian Peninsula (IP). The study was done for an ensemble of high quality rainfall time series observed in the IP during the period 1961–2010. Two approaches, peaks-over-threshold (POT) and block maxima (BM) with the generalized extreme value (GEV) distribution, were compared in order to identify which is the more appropriate for the estimation of RLs. For the first approach, which identifies trends in the parameters of the asymptotic distributions of extremes, both all-days and rainy-days-only datasets were considered because a major fraction of values of daily rainfall over the IP is zero. For the second approach, rainy-days-only data were considered showing how the mean, variance and number of rainy days evolve. The 20-year RLs expected for 2020 were estimated using these methods for three seasons: autumn, spring and winter. The GEV is less reliable than the POT because fixed blocks lead to the selection of non-extreme values. Future RLs obtained with the POT are greater than those estimated with the GEV, mainly because some gauges show significant positive trends for the number of rainy days. Autumn, rather than winter, is currently the season with the heaviest rainfall for some regions.


Introduction
The precipitation regime in the Iberian Peninsula (IP) is highly variable due to its complex topography, for which reason no wide areas are left without coverage. The spatial variability of the IP rainfall is such that certain regions receive more than 3000 mm/year, while others, for example in the southeast, receive on average less than 200 mm/year, the lowest values in Europe. The atmospheric circulation patterns over the IP change depending on the season [1]. Although precipitation is generated by different physical processes during the different times of the year [2], the intensification of rainfall during the rainy seasons is due to frontal systems coming from the Atlantic Ocean [3,4], which cause persistent rainfall. For the eastern IP, rainfall is produced by easterly flows leading to heavy convective precipitation over the Mediterranean area, especially when there is colder air at high levels. Furthermore, over extensive regions of the IP, a few rainy days concentrate much of the annual precipitation [5]. This variability leads to a study of extreme rainfall over the IP being particularly interesting.
Previous extreme rainfall studies covering particular regions in Spain, the whole of Spain or the whole IP have considered a variety of indices, all taken from the Expert Team on Climate Change Detection and Indices (ETCCDI) set of indices recommended for use by the World Meteorological Organization [6]. For the region of Andalusia, among others [7], the most noteworthy results were found in winter, with a predominance of decreasing trends for intensity indices in central and western (CW) Andalusia, but positive trends in the south-eastern (SE) area. Another study of the whole IP [8] analyzed annual and seasonal trends of extreme precipitation indices using observations for Portugal and Spain and the ERA-driven simulation.
Some studies of extreme rainfall trends for the IP have used the statistical extreme value theory (EVT) and some the block maxima (BM) [9] or the peaks-over-threshold (POT) approaches [10,11]. These methods each have advantages and disadvantages, and it is important to compare the results they give. Recently, a new method has been developed for the calculation of non-stationary return levels (RLs) for extreme rainfall in a southwestern region of the IP [12]. The aim of this present study was not to propose a method to estimate the Rls, but to compare two ways of dealing with possible trends in rainfall so as to estimate near-future extremes. This objective is tackled for a set of complete daily rainfall time series from 76 gauges over the period 1961-2010 for the whole IP.
The organization of the paper is as follows: the data used are described in Section 2; an overview of the method is given in Section 3; and the main results are presented and discussed in Section 4. Finally, the main conclusions are presented in Section 5.

Data
The study area was the Iberian Peninsula. Most of the time series were taken from an extensive database of daily rainfall time series provided by the Spanish National Meteorology Agency (In Spanish: Agencia Estatal de METerología, AEMET). The quality requirements were to choose time series with no missing data and to cover the orographic diversity of the IP.
The final choice was a set of 76 daily rainfall time series (Table 1) corresponding to gauges as regularly spaced as possible over the IP. Their locations are shown in Figure 1. The study period was 1961-2010. Most of the time series were selected from the daily rainfall time series of AEMET according to the quality requirements described above. Four of the series were provided and quality controlled by the European Climate Assessment and Dataset service (ECA, available online at http://eca.knmi.nl), and one series was from Portugal (Gafanha da Nazaré) provided by the National System of Water Resources Information (SNIRH, available online at http://snirh.apambiente.pt/, managed by the Portuguese Institute for Water) database.  Table 1.
It is often important to determine whether a set of data is homogeneous before any statistical technique is applied to it. Inhomogeneities in long-term climatological time series could be due to various factors, such as changes in instruments, observation methods, station relocation, etc. Data homogeneity was evaluated using the R-based program RHTestV3, developed at the Climate  Table 1.
It is often important to determine whether a set of data is homogeneous before any statistical technique is applied to it. Inhomogeneities in long-term climatological time series could be due to various factors, such as changes in instruments, observation methods, station relocation, etc. It is based on a two-phase regression model with a linear trend for the entire base series [13,14]. This analysis, together with the metadata of the stations showed that none of the 76 time series had change-points significant at 5%, with all of them being homogeneous in the cited period of study.
For this study of extreme precipitation over the IP, in view of the highly seasonal nature of rainfall, each season was studied separately. The definition of the seasons used was that common in climatological studies: winter was December, January and February; spring was March, April and May; and autumn was September, October and November. Summer was not considered because of the low number of rainy events during that season in most parts of the IP. For the BM approach, each season was considered for the selection of its maximum value. For the POT approach, the threshold used in defining extreme rainfall for each season separately was the 98th percentile of all the daily rainfall values for the whole period and the 95th percentile of the rainy-days-only values.

Methods
The 20-year RLs for the precipitation in the IP were studied using two different EVT approaches-the BM and the POT-considering both all-days (rainy and non-rainy) and rainy-days-only datasets. The BM approach is frequently applied to the study of extreme events of meteorological variables in the framework of EVT. Its focus is on modeling the BM of the variables with the generalized extreme value (GEV) distribution. The GEV theory is a kind of "law of large numbers". It states that the maximum of n independent and identically distributed variables with probability distribution function of type F tends to follow a GEV distribution when n tends to infinity (i.e., in practice, for large enough n).
In a similar way, the POT approach is based on the asymptotic convergence of the exceedances of a high threshold u to a general Pareto distribution (GPD) when u tends to infinity. This holds true again for independent and identically distributed variables. Making an optimal choice for the threshold u is a difficult problem. A very high value of u leads to few exceedances and consequently high variance estimators. On the contrary, a very low value of u is likely to violate the asymptotic basis of the model, leading to biases [15].
An independent and identical distribution is assumed by the theory, but is generally not verified for meteorological variables. Weak dependence can be dealt with: annual maxima can generally be considered as independent, and values above the threshold can be made independent, for example by applying a run declustering procedure [16]. An identical distribution is not true because of seasonality, interannual variability and possible trends. To handle seasonality, estimations are made for each season rather than for the whole year. Interannual variability and trends are more difficult to tackle. The estimation of RLs in a non-stationary context has been the subject of a growing number of papers in recent years. While some authors consider that taking non-stationarity into account implies too large uncertainties to be justified [17,18], others suggest methods for doing so together with new definitions for an RL in such a context [19][20][21][22][23][24].
The objective of the present work was not to derive an RL for some design purpose, but to compare two ways of dealing with possible trends in rainfall when estimating near-future extremes. The focus was therefore not on the definition of a non-stationary RL, but on the type of trends to be considered and the consequences of the choice for both the BM and the POT approaches. Thus, two different approaches were used to compute future RLs: (1) The first consists of identifying trends in the parameters of the asymptotic distributions of extremes: location and scale parameters of the GEV for the BM and threshold and scale parameters of the GPD for the POT. The trend identification is based on likelihood ratio tests at a 5% significance level between models with linear or constant parameters [15], and the confidence intervals (CI) are computed by bootstrapping in order to take the uncertainty in the trends into account (see the Appendix in [12] for details). (2) A residual process that has been explained extensively elsewhere [12] is used in order to calculate the 20-year RLs for 2020. The idea is to use trends in the main characteristics of the whole distribution rather than those in the extreme values only. If the non-stationarity of extremes in a statistical framework is explained by means and variances, the daily mean and standard deviation in 2020 are estimated by linear extrapolation of the linear trends that have been estimated from observations so as to be able to compute the RLs for 2020.
Although linear trends are poorly suited to represent future evolutions in rainfall correctly, they were chosen here as a simple way to compare two different approaches under similar conditions. The second approach is appropriate for the use of climate simulation results to compute future means and standard deviations and is preferable if the aim is to estimate future RLs as design values for an installation or project.

Results
This section presents the main results of the 20-year RL (Z20) estimations. The present day was estimated using both all-days (Section 4.1) and rainy-days-only (Section 4.2) data, and the near future was estimated for 2020 (Section 4.3).

Present 20-Year RLs Estimated from All-Days Rainfall Data
Considering all days (rainy and non-rainy), the 20-year RLs for the two approaches were estimated for the present day by using the daily rainfall time series. Figure 2 shows the results for the three seasons considered. For autumn and spring, all the 20-year RLs estimated using the POT approach lay inside the CI of the 20-year RLs using GEV. In winter, only one gauge (vlm) showed no overlapping CIs, the RL estimated with the POT (55.09 [47.96; 62.22]) being greater than that corresponding to the BM (50. 16 [45.51; 54.81]). The difference between the two values is small, but since the CI with BM is narrower than with POT, the RL estimated with the POT does not fall inside the GEV CI. since the CI with BM is narrower than with POT, the RL estimated with the POT does not fall inside the GEV CI. In order to further analyze this exception, Figure 3 shows the GEV/POT comparison for this gauge (vlm) in winter. Looking at the top figure, one observes that the frequency of threshold exceedances has decreased over the last decade, and the last values are lower than the rest. Considering the red points for GEV in greater detail, one sees that, as one maximum is chosen for each season, high and low values are mixed, some of them being lower than the POT threshold. The shape parameter is different: ξ = 0 for POT and ξ = −0.2845 for GEV. In this latter case, the distribution is strongly bounded, which is surprising for precipitation. Since sigma measures the variability, the scale parameter is greater in GEV because the maximum is more variable than in POT when choosing such low values as the maximum (a season with less than 10 mm of rainfall is not really a maximum). The estimates were made as if the time series were stationary, which is probably not the case and may further explain the differences. In order to further analyze this exception, Figure 3 shows the GEV/POT comparison for this gauge (vlm) in winter. Looking at the top figure, one observes that the frequency of threshold exceedances has decreased over the last decade, and the last values are lower than the rest. Considering the red points for GEV in greater detail, one sees that, as one maximum is chosen for each season, high and low values are mixed, some of them being lower than the POT threshold. The shape parameter is different: ξ = 0 for POT and ξ = −0.2845 for GEV. In this latter case, the distribution is strongly bounded, which is surprising for precipitation. Since sigma measures the variability, the scale parameter is greater in GEV because the maximum is more variable than in POT when choosing such low values as the maximum (a season with less than 10 mm of rainfall is not really a maximum). The estimates were made as if the time series were stationary, which is probably not the case and may further explain the differences.
Water 2018, 10, x FOR PEER REVIEW 6 of 16 since the CI with BM is narrower than with POT, the RL estimated with the POT does not fall inside the GEV CI. In order to further analyze this exception, Figure 3 shows the GEV/POT comparison for this gauge (vlm) in winter. Looking at the top figure, one observes that the frequency of threshold exceedances has decreased over the last decade, and the last values are lower than the rest. Considering the red points for GEV in greater detail, one sees that, as one maximum is chosen for each season, high and low values are mixed, some of them being lower than the POT threshold. The shape parameter is different: ξ = 0 for POT and ξ = −0.2845 for GEV. In this latter case, the distribution is strongly bounded, which is surprising for precipitation. Since sigma measures the variability, the scale parameter is greater in GEV because the maximum is more variable than in POT when choosing such low values as the maximum (a season with less than 10 mm of rainfall is not really a maximum). The estimates were made as if the time series were stationary, which is probably not the case and may further explain the differences.

Present 20-Year RLs Estimated from Rainy Days
Considering rainy-days-only and studying when Z20-POT is inside the CI of Z20-GEV, we show in Figure 4 the results for the three seasons considered. For both autumn and winter, all the gauges except one in each case (pcr and cal, respectively) have Z20-POT inside the Z20-GEV CI, and for spring, there are three exceptions (cal, gra, ssp). In all cases, Z20-POT is greater than Z20-GEV. The CI obtained using GEV is narrower than that corresponding to POT. The gauges with non-overlapping CIs and different values for the RL present the same behavior. Figure 5 illustrates the GEV/POT comparison for the gauge "pcr" in autumn. There are many maxima in the second half of the study period with the BM approach, with low rainfall leading to smaller RL values.

Present 20-Year RLs Estimated from Rainy Days
Considering rainy-days-only and studying when Z20-POT is inside the CI of Z20-GEV, we show in Figure 4 the results for the three seasons considered. For both autumn and winter, all the gauges except one in each case (pcr and cal, respectively) have Z20-POT inside the Z20-GEV CI, and for spring, there are three exceptions (cal, gra, ssp). In all cases, Z20-POT is greater than Z20-GEV. The CI obtained using GEV is narrower than that corresponding to POT.

Present 20-Year RLs Estimated from Rainy Days
Considering rainy-days-only and studying when Z20-POT is inside the CI of Z20-GEV, we show in Figure 4 the results for the three seasons considered. For both autumn and winter, all the gauges except one in each case (pcr and cal, respectively) have Z20-POT inside the Z20-GEV CI, and for spring, there are three exceptions (cal, gra, ssp). In all cases, Z20-POT is greater than Z20-GEV. The CI obtained using GEV is narrower than that corresponding to POT. The gauges with non-overlapping CIs and different values for the RL present the same behavior. Figure 5 illustrates the GEV/POT comparison for the gauge "pcr" in autumn. There are many maxima in the second half of the study period with the BM approach, with low rainfall leading to smaller RL values. The gauges with non-overlapping CIs and different values for the RL present the same behavior. Figure 5 illustrates the GEV/POT comparison for the gauge "pcr" in autumn. There are many maxima in the second half of the study period with the BM approach, with low rainfall leading to smaller RL values.

Future 20-Year RLs
In this section, we shall present the main results of calculating the 20-year RLs for 2020. It is remarkable that, in a stationary context, according to the likelihood ratio test at a 95% confidence level, the value of the shape parameter is zero for most of the stations for all three seasons considered.

Using Trends in the Parameters
We first studied the future 20-year RLs obtained through extrapolating the trends in the parameters (scale and location for GEV and scale and threshold for POT), checking whether the RLs obtained with POT lay within the CI of the RLs obtained with GEV. Figure 6 shows the results. As can be seen, there are several gauges with non-overlapping CIs: 10 for autumn, 11 for winter and 5 for spring. For these gauges, Tables 2-4 list the degree of each parameter for GEV and POT, with zero meaning no trend and one a linear trend. Most of these gauges show discrepancies between the degree obtained for the parameters using the two methods; when there is a linear trend in the scale or location parameters using GEV, there is no trend in the scale parameter using POT, and vice versa. This leads to different values of the 20-year RL for 2020 as estimated by extrapolating the parameters.  Our consideration of a varying threshold for POT was roughly equivalent to considering a linear trend in the location parameter (μ) for GEV. Thus, those situations in which μ shows a linear trend, and σ for GEV and GPD shows no trend should be similar. The gauges that disagree (such as "agr") may be a result of the sampling, because lower values are considered using the BM approach. By way of example, Figure 7 shows the behavior with the two approaches for the gauge "agr" in autumn.

Future 20-Year RLs
In this section, we shall present the main results of calculating the 20-year RLs for 2020. It is remarkable that, in a stationary context, according to the likelihood ratio test at a 95% confidence level, the value of the shape parameter is zero for most of the stations for all three seasons considered.

Using Trends in the Parameters
We first studied the future 20-year RLs obtained through extrapolating the trends in the parameters (scale and location for GEV and scale and threshold for POT), checking whether the RLs obtained with POT lay within the CI of the RLs obtained with GEV. Figure 6 shows the results. As can be seen, there are several gauges with non-overlapping CIs: 10 for autumn, 11 for winter and 5 for spring. For these gauges, Tables 2-4 list the degree of each parameter for GEV and POT, with zero meaning no trend and one a linear trend. Most of these gauges show discrepancies between the degree obtained for the parameters using the two methods; when there is a linear trend in the scale or location parameters using GEV, there is no trend in the scale parameter using POT, and vice versa. This leads to different values of the 20-year RL for 2020 as estimated by extrapolating the parameters.  Our consideration of a varying threshold for POT was roughly equivalent to considering a linear trend in the location parameter (µ) for GEV. Thus, those situations in which µ shows a linear trend, and σ for GEV and GPD shows no trend should be similar. The gauges that disagree (such as "agr") may be a result of the sampling, because lower values are considered using the BM approach. By way of example, Figure 7 shows the behavior with the two approaches for the gauge "agr" in autumn.  Figure 8 shows the comparison of the two approaches for the gauge "dto" in winter. Again, the decrease in the location parameter is greater than that in the threshold due to the lower maximum values, especially in the second half of the time series. This leads to a lower value of Z20 for 2020 when using GEV than when using POT. The same is the case for the gauge "cal".   Figure 8 shows the comparison of the two approaches for the gauge "dto" in winter. Again, the decrease in the location parameter is greater than that in the threshold due to the lower maximum values, especially in the second half of the time series. This leads to a lower value of Z20 for 2020 when using GEV than when using POT. The same is the case for the gauge "cal".  Figure 8 shows the comparison of the two approaches for the gauge "dto" in winter. Again, the decrease in the location parameter is greater than that in the threshold due to the lower maximum values, especially in the second half of the time series. This leads to a lower value of Z20 for 2020 when using GEV than when using POT. The same is the case for the gauge "cal".

Using the Stationarity Test
The next step is to estimate the 20-year RLs for 2020 obtained by extrapolating the linear trends in the daily mean and standard deviation of the amount of rain for rainy days and in the number of rainy days [12]. The hypothesis that the non-parametric temporal evolution is essentially linked to the evolutions of the mean and the variance (the methodological approach mentioned in [12]) had previously been used to test the stationarity of the extremes of standardized residuals of the time series. Table 5 shows the percentage of gauges that verified this stationarity at a 90% confidence level. The stationarity test was reasonably well satisfied by both methods: for the BM approach, 87% of the gauges satisfied the test in autumn and spring and a lower proportion, 78%, in winter; for the POT approach, greater proportions of gauges satisfied the test: 95% in autumn and spring and 92% in winter.

Using the Stationarity Test
The next step is to estimate the 20-year RLs for 2020 obtained by extrapolating the linear trends in the daily mean and standard deviation of the amount of rain for rainy days and in the number of rainy days [12]. The hypothesis that the non-parametric temporal evolution is essentially linked to the evolutions of the mean and the variance (the methodological approach mentioned in [12]) had previously been used to test the stationarity of the extremes of standardized residuals of the time series. Table 5 shows the percentage of gauges that verified this stationarity at a 90% confidence level. The stationarity test was reasonably well satisfied by both methods: for the BM approach, 87% of the gauges satisfied the test in autumn and spring and a lower proportion, 78%, in winter; for the POT approach, greater proportions of gauges satisfied the test: 95% in autumn and spring and 92% in winter. Table 5. Percentage of the stations that completely satisfied the stationarity of the extremes of the residuals for the three seasons considered, using the two methods.

GEV POT
Autumn 87% 95% Winter 78% 92% Spring 87% 95% Having checked that the stationarity test is fairly well satisfied by both methods, we then estimated the 20-year RLs for 2020. Figure 9 shows the spatial distribution of the stations with an RL obtained with POT that was within the CI of the RL obtained with GEV. The comparison showed similar results for the two methods, giving overlapping CIs for all the gauges in winter, and all except one in both autumn and spring.
Water 2018, 10, x FOR PEER REVIEW 11 of 16 Table 5. Percentage of the stations that completely satisfied the stationarity of the extremes of the residuals for the three seasons considered, using the two methods. Having checked that the stationarity test is fairly well satisfied by both methods, we then estimated the 20-year RLs for 2020. Figure 9 shows the spatial distribution of the stations with an RL obtained with POT that was within the CI of the RL obtained with GEV. The comparison showed similar results for the two methods, giving overlapping CIs for all the gauges in winter, and all except one in both autumn and spring. It should be noted that we did not extrapolate the trends in parameters, but the trends in mean and variance. These trends are the same for the means in both approaches, but not for the variances, as shall be explained below. The differences in the Z20 results may thus come from the variance trends and/or from the selected high values of Y(t): the annual maxima or excesses of the 95th percentile, which may lead to different extremes for Y(t). We took into account the proportion of rainy days with the POT, but not with GEV, so that the computation of the extremes of X(t) is different from that of Y(t).

GEV POT
To go into greater depth of the exceptions (red circles in Figure 9), let us now analyze these two gauges explicitly. The gauge that is the exception in autumn is "mfr" with Z20-POT=81.04 (64.41;86.03) outside the CI of Z20-BM=73.39 (63.96;80.43).
This gauge shows a negative (positive) trend for variance using the BM (POT), leading to different results for Z20 when extrapolating these trends. There is also a positive trend in the number of rainy days. However, the variance with the two methods is different because all-days were used It should be noted that we did not extrapolate the trends in parameters, but the trends in mean and variance. These trends are the same for the means in both approaches, but not for the variances, as shall be explained below. The differences in the Z20 results may thus come from the variance trends and/or from the selected high values of Y(t): the annual maxima or excesses of the 95th percentile, which may lead to different extremes for Y(t). We took into account the proportion of rainy days with the POT, but not with GEV, so that the computation of the extremes of X(t) is different from that of Y(t).
To go into greater depth of the exceptions (red circles in Figure 9), let us now analyze these two gauges explicitly. The gauge that is the exception in autumn is "mfr" with Z20-POT=81.04 (64.41;86.03) outside the CI of Z20-BM=73.39 (63.96;80.43).
This gauge shows a negative (positive) trend for variance using the BM (POT), leading to different results for Z20 when extrapolating these trends. There is also a positive trend in the number of rainy days. However, the variance with the two methods is different because all-days were used for the BM, but rainy-days-only for the POT. This is another drawback of the BM approach: mixing rainy and non-rainy days mixes two different processes, so that the variance mixes them, as well. Figure 10 shows the comparison between the two approaches using this gauge "mfr" in autumn in the two top panels and the evolution of the mean, standard deviation and number of rainy days in the three bottom panels (seasonal values in black, linear trends in red). This highlights the importance of considering rainy-days-only because, when all days are considered, a positive trend in variance of the rain of rainy days and in the number of rainy days can result in a negative trend in variance.
Water 2018, 10, x FOR PEER REVIEW 12 of 16 for the BM, but rainy-days-only for the POT. This is another drawback of the BM approach: mixing rainy and non-rainy days mixes two different processes, so that the variance mixes them, as well. Figure 10 shows the comparison between the two approaches using this gauge "mfr" in autumn in the two top panels and the evolution of the mean, standard deviation and number of rainy days in the three bottom panels (seasonal values in black, linear trends in red). This highlights the importance of considering rainy-days-only because, when all days are considered, a positive trend in variance of the rain of rainy days and in the number of rainy days can result in a negative trend in variance. The gauge that is the exception in spring is "reu" with Z20-POT=51.58 (36.46; 59.10) outside the CI of Z20-GEV=44.80 (33.64; 51.11). It shows a positive trend for variance with both methods, with a more positive trend in the number of rainy days producing a higher Z20 when using the POT technique.

Changes in Future 20-Year RLs
Although extrapolating recent trends may not be the best approach to computing future RLs, the comparison here was aimed at analyzing the differences induced by changes similar to those already experienced. For both techniques, BM and POT, the spatial distribution of the differences between the 20-year RLs obtained for 2020 using the latter method (stationarity test) and those The gauge that is the exception in spring is "reu" with Z20-POT=51.58 (36.46; 59.10) outside the CI of Z20-GEV=44.80 (33.64; 51.11). It shows a positive trend for variance with both methods, with a more positive trend in the number of rainy days producing a higher Z20 when using the POT technique.

Changes in Future 20-Year RLs
Although extrapolating recent trends may not be the best approach to computing future RLs, the comparison here was aimed at analyzing the differences induced by changes similar to those already experienced. For both techniques, BM and POT, the spatial distribution of the differences between the 20-year RLs obtained for 2020 using the latter method (stationarity test) and those corresponding to the present day are shown in Figure 11. The changes for the BM are on the left, and those for the POT on the right, with blue (red) meaning decreasing (increasing) values of the 20-year RLs for 2020. The scale used is the same for both approaches. For the three seasons considered, the two approaches usually give similar results, although there are some remarks that need to be made. In autumn, there is an increase in the western IP and a decrease in the east, although there is also an increase in the southeast. The increase in the northeast with the BM (left) is greater than with the POT. In winter, most of the IP presents a decrease in the 20-year RLs for 2020, with only a small area in the southeast having increasing RLs, this being more notable with the POT approach. corresponding to the present day are shown in Figure 11. The changes for the BM are on the left, and those for the POT on the right, with blue (red) meaning decreasing (increasing) values of the 20-year RLs for 2020. The scale used is the same for both approaches. For the three seasons considered, the two approaches usually give similar results, although there are some remarks that need to be made. In autumn, there is an increase in the western IP and a decrease in the east, although there is also an increase in the southeast. The increase in the northeast with the BM (left) is greater than with the POT.
In winter, most of the IP presents a decrease in the 20-year RLs for 2020, with only a small area in the southeast having increasing RLs, this being more notable with the POT approach.

Discussion
The estimation of the 20-year RLs for the present day with both the BM and the POT approaches led to similar results regardless of whether rainy-days-only or all-days (rainy and non-rainy) datasets were used. There were only a few gauges showing non-overlapping CIs (one for autumn and winter and three for winter), with the RL obtained from the POT approach being greater than that corresponding to the BM. It is important to note that, because of the large number of zeros, the asymptotic assumption in the case of the BM is weakened, i.e., the maximum in the rainfall time series is not the maximum of 90 values (which is of course already quite far from infinity), but of far fewer.
The study of the 20-year RLs for a near-future time obtained through the extrapolation of the trends in the parameters of the extreme value distribution led to similar results for a great part of the gauges using both the BM and the POT approaches, although there were some exceptions. These were mainly related to the decrease of extreme rainfall at the end of the study period. The decrease in location was greater than that in the threshold due to the lower maxima during the last part of the time series. Even if one considers the maxima obtained with the BM approach that are lower than the mobile threshold (red circles in the upper panel of Figure 7), most of them occur at the end of the observation period and also correspond to years that have no maximum precipitation value with the POT approach (black dots in the same figure). This leads to a lower 20-year RL value in 2020 using GEV than using POT when extrapolating the parameters of the extreme value distribution.
The results of the study of the changes between the present day and future 20-year RLs by using the stationarity test were, in general, similar for the BM and POT. They pointed to an important consequence: autumn, instead of winter, would be the season with most extreme rainfall in the western IP due to the sharp increase in the RLs in the future. This is coherent with previous results [9,10] indicating a significant positive trend for extreme rainfall over the western half of the IP and with the findings of a regional study [12] of an area in the southwestern IP, which showed the same behavior. An increase in the frequency of extreme precipitation in autumn for the same area of the IP has been described as being due to increased northwesterly flows [25]. However, most of these changes are not significant because, for most of the gauges, the CIs overlap. This is important because the estimation of RLs is highly uncertain, especially in the context of climate change. In the comparison of the two approaches used to compute the RLs, the stationarity test led to more robust results than using trends in the parameters of the extreme value distribution. The disadvantages of the latter were its reduced sample size, sensitivity to the threshold selection and the fact that there were negative trends in extreme rainfall events at the end of the study period when using the BM. These problems are related to the conclusions drawn in a previous work [12].
The present work has estimated future RLs by extrapolating the observed trends. This does not allow the different signals involved-climate change and inter-annual variability-to be separated. In order to better understand the impact of climate change on extreme rainfall over the IP, complementary analyses with the aid of climate simulations might be necessary. The generally linear trends identified from observed time series and then extrapolated into the future involve both climate change and inter-annual variability signals present in the period of observation. In this sense, there has been a study of how to discriminate the sources of variation in changes of several precipitation characteristics, including heavy precipitation, for the end of the Twenty-First Century over the Rhine Basin [26]. One of the conclusions was that, except for the summer season, natural variability explains a considerable proportion of the variance of the changes in precipitation and is the dominant factor in the winter season. It is noteworthy that extreme rainfall changes for the Iberian Peninsula are different from those in other regions. For instance, some studies have found evidence of increasing trends in intensity for various European regions, with the Czech Republic being an example [27].

Conclusions
We have described an EVT study to estimate non-stationary RLs of extreme rainfall for the present day and the near future over the IP. We used a complete dataset of 76 gauges for the period 1961-2010.
Two EVT approaches, BM and POT, were applied to estimate the 20-year RLs for the present day and those expected for 2020.
Two methods were used to compute future rainfall RLs with the POT approach. In the first, trends in the extremes were identified considering all-days rainfall data, taking into account a time-varying threshold based on a linear quantile regression and, when appropriate, a trend in the GPD scale parameter. In the second, the RLs were calculated considering rainy-days-only data, examining the impact of evolutions of the mean and variance and of the number of rainy days. In this second case, a novel adaptation of a stationarity test that had been designed and used for temperature time series was applied to rainfall, finding that it was indeed satisfied for most of the gauges for all three seasons considered. As in [12], this second procedure was found to have fewer disadvantages than the first and could therefore be applied in using the evolution of mean and variance projected by climate models for the estimation of future RLs.
The estimation of the present-day 20-year RLs considering all-days (rainy and non-rainy) and rainy-days-only datasets led to similar results for the two approaches, with only a few gauges giving different RLs with BM and POT, as was described in detail in the previous section. When using the BM, the asymptotic assumption is weakened because, in some areas of the IP, the rainfall time series include very many days with zero rainfall, and the BM is sometimes not a real maximum.
This study has estimated the 20-year RLs expect for 2020, with the main objective having been to compare two EVT methods, BM and POT. There were some exceptions to the generally equivalent behavior, which confirmed that POT is a better method for estimating RLs, with BM being less reliable because fixed blocks lead to the selection of non-extreme values. Furthermore, the study of the changes in the near-future RLs over the IP showed that, for some regions, autumn, not winter, is the season with the heaviest rainfall over recent decades due to the increase in the RLs for the western IP.