Effect of warming climate on extreme daily rainfall depth using non-stationary Gumbel model with temperature co-variate

In this study, non-stationary frequency analysis was carried out to apply non-stationarity of extreme rainfall driven by climate change using the scale parameter of two parameters of the Gumbel distribution (GUM) as a co-variate function. The surface air temperature (SAT) or dew-point temperature (DPT) is applied as the co-variate. The optimal model was selected by comparing AICs, and 17 of 60 sites were found to be suitable for the non-stationary GUM model. In addition, SAT was chosen as the more appropriate co-variate among 13 of the 17 sites. As a result of estimating changes in design rainfall depth with future SAT rises at 13 sites, it is likely to increase by 10% in 2040 and 18% in 2070.


INTRODUCTION
The rise in temperature changes the distribution of extreme rainfall from the physical relationship between temperature and the saturated vapor pressure of water vapor in the atmosphere (O'Gorman & Schneider 2009). As a result, global warming in recent decades has produced unprecedented extreme rainfall events around the world (Min et al. 2011). Future climate simulation results presented by many GCMs in the IPCC AR5 climate change scenarios predict an increase in extreme rainfall and flooding in a wide range of regions. However, the spatial resolution of GCMs, which is about 100-km, is not sufficient to quantify accurately extreme rainfall events at the regional level (Kunkel et al. 1999). In addition, the current climate model is reported to have many problems in terms of stably simulating extreme rainfall, although large-scale patterns of temperature changes are reliably simulated (Romps 2011). Therefore, looking at the response of extreme rainfall to global warming from appropriate analysis of observed temperature and rainfall data would be one of the rational approaches (O'Gorman 2012).
Frequency analysis has been widely used to determine the design rainfall depth for occurrence probability by linking the magnitude of extreme rainfall with an appropriate probability distribution (Katz et al. 2002).
In these design practices, the probability distribution fits best to represent observations, assuming that the underlying processes of the observations are independent, identically distributed, and do not change over time. However, naturally occurring climate variability for decades has had a significant impact on the magnitude and frequency of extreme rainfall (Towler et al. 2010), and long-term trends such as global warming have made the stationary assumption no longer valid (Milly et al. 2008). In other words, from the point of view of analyzing the effects of climate change, the term 'stationarity' is no longer a valid concept in the frequency analysis, and a new method for probabilistic interpretation of extreme rainfall is needed.
Several methods have been proposed to deal with time series non-stationarity (Cunha et al. 2011;Yilmaz & Perera 2013;Jang et al. 2015;Moon et al. 2016), and many studies on changes in design rainfall depth or its return level under non-stationary conditions have been conducted (Salvadori & DeMichele 2010;Graler et al. 2013;Hassanzadeh et al. 2013;Salas & Obeysekera 2013;Shin et al. 2014;Choi et al. 2019). In order to explain the non-stationary in the frequency analysis and to examine the changes in the extreme rainfall probability distribution with temperature rise, studies have been conducted to express the parameters of the probability distribution using co-variates (Coles et al. 2001). Recently, studies on non-stationary frequency analysis using meteorological variables in the annual maximum rainfall time series have been attempted (El Adlouni et al. 2007). Aissaoui-Fqayeh et al. (2009) developed the non-stationary generalized extreme value model of annual maximum rainfall time series in California, USA by applying the Southern Oscillation Index as a climate variable. Villarini & Vecchi (2012) performed a frequency analysis of precipitation and flood volume by expressing the parameters of the Gumbel distribution as a function of the North Atlantic Oscillation Index and the Mediterranean Index. In particular, if climate variables reflecting large-scale atmospheric circulations simulated in GCMs are applied, they may be useful for examining the effects of local extreme rainfall on climate change (Tramblay et al. 2013). Also, studies have recently been conducted to analyze extreme rainfall using the link between rainfall and temperature (Towler et al. 2010;Tramblay et al. 2013).
The purpose of this study was to develop a non-stationary Gumbel model with temperature variables that can be applied to extreme rainfall events observed at major sites in Korea. The reason for applying the Gumbel distribution is because the Gumbel distribution has been the most widely used among many probability distributions for frequency analysis in Korea, and it is currently designated as the standard probability distribution for frequency analysis in the design of various hydraulic structures. In this study, daily surface air temperature (SAT) and dew-point temperature (DPT) were considered in the frequency analysis of annual maximum daily rainfall time series (AMR) for ease of application to climate change scenarios. The degree of change in design rainfall depth with respect to future temperature rises is projected using the non-stationary frequency model and projected temperature rise scenarios.

Meteorological data
In this study, daily weather data from 60 observation sites were used. In order to perform frequency analysis, AMR was extracted from daily rainfall data from 1973 to 2017. In addition, as a co-variate for non-stationary frequency analysis, SATs and DPTs up to 6-day prior to the day AMR occurred were extracted. The locations of the major observation sites are shown in Figure 1. Information and locations of the sites used in this study are provided in Figure S1 and Table S1. In Figure 1, the dot represents the linear trend of AMR (red dot : positive trend, blue dot : negative trend), and the square represents the linear trend of the mean DPT time series from June to September ( : positive trend, : negative trend). The diamond represents the linear trend of the mean SAT time series from June to September ( : positive trend, : negative trend). The AMR trend shows an average increase of 0.4 mm/year, and the SAT trend shows an average increase of 0.02°C/year. The DPT trend shows a mixture of sites that tend to increase and decrease.

Gumbel distribution
The cumulative distribution function F GUM of the Gumbel distribution is defined as follows: where x is the AMR, a is the scale parameter, and x o is the location parameter. Several methods can be used to estimate parameters of this distribution. In this study, the maximum likelihood method is applied. Parameters are estimated using optimization techniques to minimize negative log likelihood (nllh) as follows: where n is the number of data. After parameters are estimated, the daily rainfall depth x (1Àp) corresponding to a non-exceedance probability 1 À p can be calculated as follows: where 0 , p , 1. Traditionally, Gumbel distribution assume that observations are stationary, but this assumption can be mitigated by introducing a co-variate to account for the non-stationarity. For example, one of parameters in Gumbel distribution can be expressed as a function of the given co-variate. In theory, all parameters of Gumbel distribution can be applied as functions of various co-variates, but in this study, the scale parameter a was applied as a function of the co-variate z, as shown in Equation (4), to clarify the effect of co-variate. In general, various studies consider co-variates in the form of exponential functions as shown in the following equation (Tramblay et al. 2011(Tramblay et al. , 2013Razmi et al. 2017;Patra 2020): where A and B are estimated in the process of minimizing nllh in Equation (5), and quantiles of the resulting Gumbel distribution have various values depending on the co-variate z. i ¼ 1, 2, 3 Á Á Á before 6 days : In monsoon climatic regions such as Korea, because of the cooling effect of rainfall, SAT and DPT several days before the AMR day may affect the value of AMR (Ali & Mishra 2017). Therefore, the largest value of SAT (or DPT) from the day of AMR to the preceding t-day was considered to be related to AMR. In this study, t-day from 0-day to the previous 6-day was used.
In this study, the preceding t-day was applied from 0-day to 6-day. In other words, each site is composed of seven non-stationary Gumbel models related to the SAT, seven non-stationary Gumbel models related to DPT, and one stationary Gumbel model using Equation (4) and the preceding t-day. A model that minimizes Akaike Information Criterion (AIC) is selected as the optimal model among the 15 models, and AIC is calculated as follows (Akaike 1974): where nllh is determined from Equation (2) for the stationary model and from Equation (5) for the non-stationary models, and K is the number of parameters of the model (K ¼ 2 for the stationary model and K ¼ 3 for nonstationary models).

RESULTS AND DISCUSSION
Non-stationary frequency analysis The co-variate Gumbel distributions using SAT or DPT were fitted for AMR. Table 1 shows the parameters for 15 models at the Namhae site. As shown in Table 1, the model using a t ¼ exp(A þ B Á SAT tÀ2 ) showed the minimum AIC value. Therefore, it can be seen that the optimal model of AMR at the Namhae site is a model using the SAT two days before the AMR event as a co-variates were fitted to AMR. In this way, a model representing the minimum AIC value for each AMR of 60 sites was selected as the optimal model. The AMR of the 17 sites revealed that the nonstationary model was more appropriate than the stationary model, and among them, the number of sites for which DPT was selected as the co-variate was 5. The AMR of the remaining 43 sites did not indicate that the non-stationary model was superior to the stationary model (Table 2). Table 4 shows more detailed information about the selected sites. Figure 2 shows the stationary Gumbel cumulative distribution ('MLS' in Figure 2) and non-stationary Gumbel cumulative distributions ('MLNSAT' and 'MLNDPT' in Figure 2), each constructed using the average value of SAT or DPT at the Namhae site. It can be seen that as the return period increases, the deviation between the stationary model and the non-stationary model increases. The gray line is 95% confidence intervals.
In order to compare the results of the stationary and non-stationary models, the daily design rainfall depths corresponding to 90 and 99% non-exceedance probabilities are shown with the observations in Figure 3 (Namhae site). In Figure 3, 'obs' is the observed AMR. Not surprisingly, the design daily rainfall depths of the stationary model ('xps90' and 'xps99' in Figure 3) do not change over time and appear as horizontal lines (dotted lines in Figure 3). Conversely, non-stationary models ('xpn90' and 'xpn99' in Figure 3) can be seen to have various design daily rainfall depths depending on the corresponding co-variate SAT (solid lines in Figure 3).

Corrected Proof
By comparing the probability distributions for stationary and non-stationary models for the year with the smallest annual maximum daily rainfall depth (e.g., 139.8 mm at Namhae sites in 1980) and the year with the largest annual maximum daily rainfall depth (e.g., 231.5 mm at Namhae sites in 2011), the variation in design rainfall with respect to the SAT can be explained more clearly (see Figure 4). The non-stationary probability density function ('NS (1980)' in Figure 4) was compared with the stationary probability density function ('stationary' in Figure 4) for the year 1980 with a relatively small maximum daily rainfall depth. The positions of the modes of the two probability density functions are similar to each other, but the spread of the non-stationary probability density function is much smaller. By contrast, the non-stationary probability density function ('NS (2011)' in Figure 4) for the year in which the relatively  If the data follow the Gumbel distribution, the larger the standard deviation of the data, the larger the value of the scale parameter. In other words, in the non-stationary model applied in this study, the higher the SAT, the larger the scale parameter value. The concept of the stationary model serves to increase the standard deviation of the data. Changes in scale parameter due to SAT can be clearly seen when looking at design rainfall depths corresponding to non-exceedance probabilities. As shown in Figure 4, the design rainfall depth of the non-stationary model for the SAT of 1980, corresponding to a nonexceedance probability of 99%, is 358.9 mm, and the design rainfall depth of the non-stationary model for the SAT in 2011 is 491.1 mm. The corresponding design rainfall depth of the stationary model is 414.8 mm independent of the SAT. In other words, as the SAT increases, the value of the scale parameter of the Gumbel distribution increases, which leads to an increase in the design rainfall depth with the effect of increasing the standard deviation.  In this section, the proposed co-variate-based non-stationary model is used to analyze the effect of future air temperature rise on design rainfall depth. In this study, the HadGEM3-RA regional climate model (RCM) produced based on the HadGEM2-AO global climate model (GCM) developed by the United Kingdom Met Office Hadley Center was used. The spatial resolution of this model is 12.5 km. KMA (2017) divided the future period into 2011-2040 and 2041-2070 for the four RCP scenarios and presented the results of the projections of air temperature rise for each period as shown in Table 3.

Corrected Proof
To estimate the future design rainfall, the change in future SAT (or DPT) provided in the RCP scenario was applied to the model selected at 17 points (see Table 3). The scale parameter a was recalculated using the changed temperature, and the changed future probability rainfall value was applied. Figure 5 shows the change in daily design rainfall depth over the 100-year return period of 2040 and 2070 under various RCP climate change scenarios at Suwon site. Although the results vary by RCP climate change scenario, the design rainfall depth will increase by 6-13% in 2040, and it is likely to increase from 13% to 26% in 2070. Figure 6 shows the rate of change in design rainfall depth by RCP climate change scenarios and future periods at 12 sites with SAT as the co-variate. At all sites, the future design rainfall depth was found to increase. Overall, it is necessary to prepare for an increase in design rainfall depth of about 10% by 2040 and about 18% by 2070.
As a result of analyzing the rate of change in the design rainfall compared to the average present for each RCP scenario and period, the results shown in the table below were obtained. As of 2040, RCP 2.6 shows the largest increase in design rainfall, but as of 2070, the RCP 8.5 scenario result shows the largest increase compared to the present.

CONCLUSION
In this study, a non-stationary Gumbel model is proposed in which climate variables act for frequency analysis of the daily annual maximum rainfall depth time series. Candidate climate variables are SAT and DPT for the purpose of facilitating the application of climate change scenario data. Various candidate models were constructed to best reflect the variability of extreme rainfall with SAT or DPT, and the optimal model was selected using the AIC. In 17 of the 60 sites in Korea, the non-stationary Gumbel model was found to be more suitable than the stationary Gumbel model. In addition, among the 17 sites where the non-stationary model was more appropriate, SAT was chosen to be a better co-variate than DPT at 12 sites except for five sites. Applying future climate change scenarios to 13 sites where the SAT was adopted as a co-variate, the design rainfall depth is likely to increase to about 10% in 2040 and about 18% in 2070.
The presented results are significant in that the rainfall data were analyzed by using the relationship with the temperature (dew-point temperature) data and the estimation of the extreme rainfall through the non-stationary frequency analysis, not the stationary frequency analysis. In addition, it is meaningful to predict future changes in rainfall by applying future climate change scenarios based on these results.
However, since the reliability of rainfall simulated by climate models is not high, the uncertainty about the estimated future design rainfall depth is also very high. The results presented in this study are very advantageous for practical policy use because they maintain the consistency of change rates between sites. However, since the results presented in this study are derived from the Gumbel distribution alone, it is difficult to generalize them. In addition, there will be a need to consider regional frequency analysis. Nevertheless, the methodology of this study is expected to provide an important direction in establishing climate change adaptation policies against extreme rainfall.