Uncertainty of Rate of Change in Korean Future Rainfall Extremes Using Non-Stationary GEV Model

Interest in future rainfall extremes is increasing, but the lack of consistency in the future rainfall extremes outputs simulated in climate models increases the difficulty of establishing climate change adaptation measures for floods. In this study, a methodology is proposed to investigate future rainfall extremes using future surface air temperature (SAT) or dew point temperature (DPT). The non-stationarity of rainfall extremes is reflected through non-stationary frequency analysis using SAT or DPT as a co-variate. Among the parameters of generalized extreme value (GEV) distribution, the scale parameter is applied as a function of co-variate. Future daily rainfall extremes are projected from 16 future SAT and DPT ensembles obtained from two global climate models, four regional climate models, and two representative concentration pathway climate change scenarios. Compared with using only future rainfall data, it turns out that the proposed method using future temperature data can reduce the uncertainty of future rainfall extremes outputs if the value of the reference co-variate is properly set. In addition, the confidence interval of the rate of change of future rainfall extremes is quantified using the posterior distribution of the parameters of the GEV distribution sampled using Bayesian inference.


Introduction
In most studies investigating future rainfall extremes, it is assumed that future rainfall data simulated from global climate models (GCMs) or regional climate models (RCMs) are "observed" data in the future [1,2]. However, even recently developed climate models simulate temperature stably [3][4][5] but still expose many limitations to the simulation of rainfall extremes [6][7][8][9][10]. As can be found in the results of Kim et al. [11], very different future projection outputs are obtained depending on which RCM-simulated future rainfall data are used to analyze rainfall extremes. Recently, a series of attempts have been made to investigate future rainfall extremes based on the relationship between temperature and rainfall extremes [12][13][14]. The basic concept commonly included in these attempts is based on the fact that investigating the behavior of rainfall extremes under global warming conditions from the relationship between observed temperature and observed rainfall extremes is an approach to obtain more consistent future outputs [15,16]. co-variate-based non-stationary frequency analysis model, and it identifies the uncertainty of the rate of change of future rainfall extremes.
This study proceeds as follows: (1) A non-stationary frequency analysis of the annual maximum daily rainfall time series (AMR) using surface air temperature (SAT) or dew point temperature (DPT) as a co-variate is performed. (2) Using the Bayesian approach, parameters of the stationary and non-stationary generalized extreme value (GEV) distributions for the AMR are sampled from the posterior distribution. (3) Using this, we evaluate the performance of the stationary frequency analysis and the non-stationary frequency analysis in terms of uncertainty. (4) We investigate how uncertainty in the non-stationary frequency analysis can be reduced by determining the appropriate co-variate reference value corresponding to the rainfall quantile. (5) We present the rate of change in rainfall quantile under climate change scenario conditions using future information of SAT or DPT produced from a combination of various climate models. (6) The uncertainty of the rate of change of future rainfall extremes is analyzed using the information obtained from the non-stationary frequency analysis by the Bayesian approach.

Data
We used the daily mean SAT, daily mean DPT, and daily precipitation of two Automated Surface Observing System (ASOS) meteorological observation sites (Chuncheon and Cheonan) operated by the Korea Meteorological Agency (KMA) [40]. The Chuncheon site is located in a mountainous area of the central part of the Korean Peninsula, and the Cheonan site is located in a plain area of the western part of the Korean Peninsula. The data period was from 1973 to 2019. Figure 1 shows the spatial location of the sites used.
Atmosphere 2021, 12, x FOR PEER REVIEW 3 of 17 ity generated from various climate model combinations and the uncertainty of the co-variate-based non-stationary frequency analysis model, and it identifies the uncertainty of the rate of change of future rainfall extremes. This study proceeds as follows: (1) A non-stationary frequency analysis of the annual maximum daily rainfall time series (AMR) using surface air temperature (SAT) or dew point temperature (DPT) as a co-variate is performed. (2) Using the Bayesian approach, parameters of the stationary and non-stationary generalized extreme value (GEV) distributions for the AMR are sampled from the posterior distribution. (3) Using this, we evaluate the performance of the stationary frequency analysis and the non-stationary frequency analysis in terms of uncertainty. (4) We investigate how uncertainty in the nonstationary frequency analysis can be reduced by determining the appropriate co-variate reference value corresponding to the rainfall quantile. (5) We present the rate of change in rainfall quantile under climate change scenario conditions using future information of SAT or DPT produced from a combination of various climate models. (6) The uncertainty of the rate of change of future rainfall extremes is analyzed using the information obtained from the non-stationary frequency analysis by the Bayesian approach.

Data
We used the daily mean SAT, daily mean DPT, and daily precipitation of two Automated Surface Observing System (ASOS) meteorological observation sites (Chuncheon and Cheonan) operated by the Korea Meteorological Agency (KMA) [40]. The Chuncheon site is located in a mountainous area of the central part of the Korean Peninsula, and the Cheonan site is located in a plain area of the western part of the Korean Peninsula. The data period was from 1973 to 2019. Figure 1 shows the spatial location of the sites used. The Coordinated Regional Climate Downscaling Experiment East-Asia (CORDEX-EA) team recently produced future climate information for the East Asian region from various GCMs and RCMs. There are two main types of future climate outputs developed by CORDEX-EA: 50-km spatial resolution outputs for East Asia and 12.5-km spatial resolution outputs for Korea. Future climate outputs with a spatial resolution of 12.5 km, produced from a total of 8 GCM-RCM combinations from 2 GCMs, including the Hadley Center Global Environment Model version 2-coupled Atmosphere-Ocean model (HadGEM2-AO, H2) and the Max Planck Institute Earth System Model-Low Resolution (MPI-ESM-LR, E6), and 4 RCMs, including the Fifth-Generation Penn State/NCAR Mesoscale Model (MM5), the Regional Spectral Model (RSM), RegCM4 (the latest version of the Regional Climate Model system), and the Weather Research and Forecasting model The Coordinated Regional Climate Downscaling Experiment East-Asia (CORDEX-EA) team recently produced future climate information for the East Asian region from various GCMs and RCMs. There are two main types of future climate outputs developed by CORDEX-EA: 50-km spatial resolution outputs for East Asia and 12.5-km spatial resolution outputs for Korea. Future climate outputs with a spatial resolution of 12.5 km, produced from a total of 8 GCM-RCM combinations from 2 GCMs, including the Hadley Center Global Environment Model version 2-coupled Atmosphere-Ocean model (HadGEM2-AO, H2) and the Max Planck Institute Earth System Model-Low Resolution (MPI-ESM-LR, E6), and 4 RCMs, including the Fifth-Generation Penn State/NCAR Mesoscale Model (MM5), the Regional Spectral Model (RSM), RegCM4 (the latest version of the Regional Climate Model system), and the Weather Research and Forecasting model (WRF), were used in this study. In addition, Representative Concentration Pathway (RCP) 4.5 and RCP 8.5 were applied as future climate change scenarios. That is, a total of 16 future climate ensembles were used. For reference, the climate information outputs produced from a combination of climate modes consist of present data from 1981 to 2010 and future data from 2021 to 2050 [41]. The climate model outputs' information is summarized in Table 1. All data of daily outputs produced by the combination of climate models were applied after bias correction by the quantile mapping process [42] linked to the Boe method [43].

Non-Stationary GEV Distribution
In this study, GEV distribution was applied as a probability distribution for frequency analysis of the AMR. The probability density function f(x) of the GEV distribution is as follows: where α is the scale parameter, β is the shape parameter, and x o is the location parameter. From Equation (1), one can calculate the quantile x T corresponding to the return level of T-year in the GEV distribution as follows: Traditionally, the GEV distribution assumes that observations are independent and identically distributed. However, in order to model non-stationarity in the AMR time series, many researchers allowed the parameters of the GEV distribution to change depending on the co-variate [10,22,[27][28][29]31,32,36,39]. In theory, all parameters of the GEV distribution can be applied as a function of various co-variates, but in this study, a co-variate was applied only to the scale parameter alpha, as shown below, in order to intuitively understand the effect of the co-variate. Non-stationary model n 1 : Non-stationary model n 2 : where SAT(y) is the SAT observed on the day on which the AMR occurred in year y, and DPT(y) is the DPT observed on the day on which the AMR occurred in year y. Therefore, one can calculate the quantile x T corresponding to the return level of the T-year of the non-stationary model as follows: Atmosphere 2021, 12, 227

of 16
That is, the quantile corresponding to the return level of the T-year of the nonstationary model has various values depending on the co-variate. Since the current hydraulic structure design practice is based on the fact that the quantile corresponding to the return level of the T-year in the stationary model has a single value, the dependence of the quantile on the co-variate in non-stationary models acts as a major cause that makes practical application of non-stationary models difficult.

Parameter Estimation and Uncertainty
The parameters of the GEV distribution were estimated using the Metropolis-Hastings (MH) algorithm. Table 2 shows the parameters and negative log likelihood (nllh) estimated by fitting each of the three models to the observed AMR. The bold numbers denote the model representing the smallest nllh. The AMR of the Chuncheon site showed the best fit when applying the non-stationary model n 1 (i.e., using SAT as a co-variate of the scale parameter). The AMR of the Cheonan site showed the best fit when applying the nonstationary model n 2 (that is, using DPT as a co-variate of the scale parameter). However, the AIC (where AIC = 2 × (# of parameters) + 2 × nllh) reflecting the parsimony point of view recommends that it is more appropriate to apply the stationary model s for both AMRs of both sites. That is, as applied in Lee et al. [10] and Ouarda et al. [36], if an optimal model is selected using the AIC, which mainly considers the aspect of fitness of the observed AMR, the stationary model s is selected as the optimal model for the AMR at Chuncheon and Cheonan sites. The fundamental purpose of performing a frequency analysis is to estimate the rainfall quantile. However, since the parameters of the probability distribution necessary for estimation of the quantile are estimated from limited samples, uncertainty is bound to be included in the parameter of the probability distribution and the corresponding rainfall quantile. Therefore, if the sample size changes (i.e., when new data become available or when sampling is performed through other techniques), the selection result of the optimal frequency analysis model centered on fitness with a limited sample such as AIC is likely to change. Our argument is that the uncertainty of the probability distribution parameter and rainfall quantile gives important information in selecting the optimal frequency analysis model. In other words, identifying in which model the uncertainty of the parameter and rainfall quantile is the lowest can be an important determinant in selecting an optimal model. This is based on the fact that a model with less uncertainty in the parameter and rainfall quantile can be recognized as a more reliable model. We quantified the uncertainty of the model as follows: Atmosphere 2021, 12, 227 6 of 16 where 95PPU (Percent Prediction Uncertainty) means the 95% confidence interval of the corresponding variable [44]. A total of 10,000 parameters were sampled from the posterior distribution of parameters using Markov chain Monte Carlo (MCMC), and 10,000 rainfall quantile ensembles corresponding to the return level of 50 years were simulated. Equations (6) and (7) were used to quantify the uncertainty about the parameter and the uncertainty about the rainfall quantile. Table 3 shows the results. The parameters of the stationary model s are α, β, and x o , while the parameters of the non-stationary models (n 1 and n 2 ) are α 1 , α 2 , β, and x o . For direct comparison between models, the p-factor converted from α 1 and α 2 of non-stationary models to α is also expressed. The α of the non-stationary model was calculated using is the average of the SAT or DPT corresponding to each AMR, respectively. The bold numbers denote the parameter representing the smallest uncertainty. If a reliable model is selected based on the uncertainty of parameters, it can be recognized that the optimal model for the AMR at the Chuncheon site is the non-stationary model n 2 , and the optimal model for the AMR at the Cheonan site is the stationary model s. For the AMR at the Chuncheon site, the p-factor of the scale parameter and the location parameter was the smallest in the non-stationary model n 2 . For the AMR at the Cheonan site, the p-factors of the scale parameter and the shape parameter were the smallest in the stationary model s. However, the p-factor of all parameters in the optimal model is not smaller than the p-factor of parameters in other candidate models. Selecting the optimal model based on the p-factor of parameters may involve a lot of subjectivity.
The q-factor of the rainfall quantile corresponding to the return level of 50 years was calculated in two ways. The first method was to calculate the q-factor using the rainfall quantile ensemble reflecting all observed co-variate values. The second method was to calculate the q-factor of the rainfall quantile ensemble under the condition of E[SAT] or E[DPT] (numbers in parentheses in Table 3). Since the q-factor by the first calculation method adds co-variate change to the uncertainty caused by parameter estimation, the q-factor from the first calculation method becomes larger than the q-factor of the second calculation method. This fact was also discussed in Ouarda et al. [37]. However, it can be found that the difference between the values of the q-factor calculated by the first method and the second method is not very large. This suggests that in the range of values of SAT or DPT observed routinely, the effect of the change in co-variate on the uncertainty of the rainfall quantile is limited.
If the optimal model is selected using the uncertainty of the rainfall quantile, that is, the q-factor, it can be recognized that the optimal model for AMR at both sites is the non-stationary model n 2 . Table 4 shows the results of selecting the optimal model using nllh, AIC, p-factor, and q-factor. The penalty for the number of parameters in AIC is to evaluate the model by reflecting the uncertainty of the model. However, evaluating the model using the p-factor or q-factor is a more direct reflection of the uncertainty. As described above, since the optimal model for the p-factor may vary for each parameter, it is more preferable to use the q-factor rather than the p-factor to evaluate the reliability of the model. Our proposal is to evaluate the fitness of the model for AMR using nllh and the reliability of the model using the q-factor. Therefore, the optimal model for AMR at the Chuncheon site can be determined as the non-stationary model n 1 or n 2 . Similarly, it can be said that the non-stationary model n 2 is recommended as the optimal model for AMR at the Cheonan site. Table 4. Results of determining optimal model.

Site nllh AIC p-Factor q-Factor
Chuncheon n 1 s n 2 n 2 Cheonan n 2 s s n 2 Figure 2 shows the response of the rainfall quantile to the change in co-variate and the probability distribution of the observed co-variate. The ensemble average of the rainfall quantile by model s represents a single value regardless of the co-variate, and the width of the 95% confidence interval is not related to the co-variate. However, it should be noted that the ensemble average and 95% confidence interval by models n 1 and n 2 respond to the values of the co-variates observed on the day of AMR. The range of the ensemble of rainfall quantiles by non-stationary models is likely to be wider than the range of the ensemble of rainfall quantiles by the stationary model since the co-variate change is additionally reflected in addition to parameter uncertainty. In other words, the uncertainty of the rainfall quantiles by non-stationary models needs to be recalculated by reflecting only the uncertainty due to parameters by designating an appropriate co-variate value. The question is how to specify the appropriate co-variate value. We pay attention to Equations (2) and (5) to find the appropriate co-variate value, that is, the reference co-variate value. Equation (5) returns the rainfall quantile corresponding to the return level of the T-year in the non-stationary GEV distribution. Equa- We pay attention to Equations (2) and (5) to find the appropriate co-variate value, that is, the reference co-variate value. Equation (5) returns the rainfall quantile corresponding to the return level of the T-year in the non-stationary GEV distribution. Equation (2) returns the rainfall quantile of the stationary GEV distribution. Equations (2) and (5) can be used to calculate the value of the reference co-variate of the non-stationary GEV distribution, which returns the same rainfall quantile as the stationary GEV distribution. Reference co-variates for various return levels of AMR at the Chuncheon site and the Cheonan site are shown in Figure 3. In this study, based on the observed AMR, a method of estimating reference co-variates corresponding to various return levels is proposed under the assumption that the rainfall quantiles for each return level of the stationary model and the non-stationary model are the same. However, the relationship between the reference co-variate and the return level is not yet clear. In addition, the uncertainty of the reference co-variate should be further revealed. Reference co-variates may be an important research topic in co-variate-based non-stationary frequency analysis in the future. The rainfall quantile 50 corresponding to the return level of 50 years is shown. In the upper plot, the solid black line is the ensemble average of 50 for the stationary model s, and the area painted in green is the 95% confidence interval. Blue lines were drawn from the non-stationary model 1 . The red lines are the result of the non-stationary model 2 . The empirical probability distribution of the observed co-variate is shown in the plot at the bottom. The distribution of surface air temperature (SAT) is indicated by the blue line and the distribution of dew point temperature (DPT) is indicated by the red line.
We pay attention to Equations (2) and (5) to find the appropriate co-variate value, that is, the reference co-variate value. Equation (5) returns the rainfall quantile corresponding to the return level of the T-year in the non-stationary GEV distribution. Equation (2) returns the rainfall quantile of the stationary GEV distribution. Equations (2) and (5) can be used to calculate the value of the reference co-variate of the non-stationary GEV distribution, which returns the same rainfall quantile as the stationary GEV distribution. Reference co-variates for various return levels of AMR at the Chuncheon site and the Cheonan site are shown in Figure 3. In this study, based on the observed AMR, a method of estimating reference co-variates corresponding to various return levels is proposed under the assumption that the rainfall quantiles for each return level of the stationary model and the non-stationary model are the same. However, the relationship between the reference co-variate and the return level is not yet clear. In addition, the uncertainty of the reference co-variate should be further revealed. Reference co-variates may be an important research topic in co-variate-based non-stationary frequency analysis in the future.  Although not covered in this study, the bivariate frequency analysis related to covariates could be one of the important research topics. Although not related to co-variates, the studies of Shiau [45], which applied a bivariate extreme value distribution to model extreme floods for flood volumes and flood peaks, Serinaldi [46], which emphasized a clear understanding of the return level in non-stationary frequency analysis, and De Luca et al. [47], which evaluated the return level of a design hyetograph using the concept of copula, can be said to provide good insight in the context of future expansion of this topic.

Future Rainfall Extremes
In most studies, including those of Hosseinzadehtalaei et al. [17] and Kim et al. [11], the future rainfall extremes taking into account climate change scenarios are estimated based on the ratio of rainfall quantiles under present climate conditions to rainfall quantiles under future climate conditions (i.e., rate of change). When estimating the rainfall quantile in present or future climate conditions, stationary conditions are generally assumed for each period (e.g., 1980-2010, 2021-2050, and 2071-2100). That is, the present rainfall quantile is estimated under the assumption that the AMR produced from the climate model under the present climate condition satisfies the stationary condition, and the future rainfall quantile is estimated under the assumption that the AMR produced under the future climate condition also satisfies another stationary condition. Figure 4a shows the results of projecting the rate of change of future rainfall extremes using this approach. The rate of change is defined as [future rainfall quantile]/[present rainfall quantile]. That is, a Atmosphere 2021, 12, 227 9 of 16 rate of change of 1.2 means that the future rainfall quantile will increase by 20% over the present rainfall quantile.
sumed for each period (e.g., 1980-2010, 2021-2050, and 2071-2100). That is, the present rainfall quantile is estimated under the assumption that the AMR produced from the climate model under the present climate condition satisfies the stationary condition, and the future rainfall quantile is estimated under the assumption that the AMR produced under the future climate condition also satisfies another stationary condition. Figure 4a shows the results of projecting the rate of change of future rainfall extremes using this approach. The rate of change is defined as [future rainfall quantile]/[present rainfall quantile]. That is, a rate of change of 1.2 means that the future rainfall quantile will increase by 20% over the present rainfall quantile. Based on the ensemble average of the rate of change by eight climate model combinations, the rainfall extremes in the future period (2021-2050) are likely to increase from 8% to 24% over the rainfall extremes in the present period (1981-2010). However, depending on which climate model combination is applied, the rate of change has very different values (0.67-1.81). It can even be found that cases where the rate of change is greater than 1 (i.e., rainfall extremes increase) and less than 1 (i.e., rainfall extremes decrease) are occurring simultaneously. Figure 4b,c show the rate of change in annual average SAT and DPT. Compared with rainfall extremes, it can be seen that the fluctuation range of the rate of change for each combination of climate models is significantly reduced. In addition, Based on the ensemble average of the rate of change by eight climate model combinations, the rainfall extremes in the future period (2021-2050) are likely to increase from 8% to 24% over the rainfall extremes in the present period (1981-2010). However, depending on which climate model combination is applied, the rate of change has very different values (0.67-1.81). It can even be found that cases where the rate of change is greater than 1 (i.e., rainfall extremes increase) and less than 1 (i.e., rainfall extremes decrease) are occurring simultaneously. Figure 4b,c show the rate of change in annual average SAT and DPT. Compared with rainfall extremes, it can be seen that the fluctuation range of the rate of change for each combination of climate models is significantly reduced. In addition, regardless of the combination of climate models applied, the rate of change shows a consistent increase pattern. These results, coupled with the results of Section 3.1 which state that temperature is the recommended co-variate for the estimation of rainfall extremes, can be said to suggest a new approach to producing more reliable future rainfall extremes outputs.
A trend analysis of SAT and DPT produced from various climate model combinations was performed to estimate future rainfall extremes reflecting climate change scenarios from co-variate-based non-stationary models derived from observations described in Section 3.1. The trend was assumed to be linear, as shown in Equations (8) or (9). The yearly SAT and DPT calculated from the daily SAT and DPT time series for the present and future periods produced for each climate model combination were used for the linear trend analysis.
where Y is the future year, and ∆SAT Y and ∆DPT Y are the annual mean SAT or DPT at future year Y minus the average of SAT or DPT for the present period (1981-2010), respectively. In Equations (8) and (9), 2010 means the last year of the present period.
From the trend analysis of SAT and DPT, one can set the scale parameter of the GEV distribution for a specific year in the future as follows: Non-stationary model n 1 : Non-stationary model n 2 : where SAT r and DPT r are the reference SAT and DPT for the return level obtained in Figure 3, respectively. Figure 5 shows the time evolution of the scale parameter of model n 2 at the Cheonan site. It can be found that the scale parameter increases gradually in the future.

Decision Making from Ensemble Average and Uncertainty
How to reflect the impact of climate change in the estimation of the design rainfall depth is a task in the new normal era facing hydrologists. The most common way to investigate changes in future rainfall extremes is to use future rainfall data derived from climate models with climate change scenarios applied. In the early days, future rainfall extremes were estimated using rainfall data predicted from one climate model [19], but it has been established that it is desirable to use ensembles simulated in various climate models [20]. To quantify uncertainty in future climate change with a probabilistic approach, an ensemble of various climate model combinations is needed [1]. The reason why many researchers use the ensemble approach and try to obtain a better ensemble is to determine what the rate of change of future rainfall extremes will be in the end. In terms of climate change adaptation policy, the rate of change of the design rainfall depth at a specific duration for a specific return level at a specific site should eventually be given as a single value. This single value will be the ensemble average no matter what method we

Decision Making from Ensemble Average and Uncertainty
How to reflect the impact of climate change in the estimation of the design rainfall depth is a task in the new normal era facing hydrologists. The most common way to investigate changes in future rainfall extremes is to use future rainfall data derived from climate models with climate change scenarios applied. In the early days, future rainfall extremes were estimated using rainfall data predicted from one climate model [19], but it has been established that it is desirable to use ensembles simulated in various climate models [20]. To quantify uncertainty in future climate change with a probabilistic approach, an ensemble of various climate model combinations is needed [1]. The reason why many researchers use the ensemble approach and try to obtain a better ensemble is to determine what the rate of change of future rainfall extremes will be in the end. In terms of climate change adaptation policy, the rate of change of the design rainfall depth at a specific duration for a specific return level at a specific site should eventually be given as a single value. This single value will be the ensemble average no matter what method we take [2,11]. However, as shown in previous studies [10,17] and Figure 4a of this study, the future rainfall extremes estimated by the stationary frequency analysis method using future rainfall data showed a lot of variation depending on which climate model combination was applied. In other words, the inter-model variability was very high. Even the future rainfall extremes of a specific site tend to increase in some combinations of climate models and decrease in other combinations of climate models. If the deviation between the ensemble members constituting the ensemble average is too large (i.e., the uncertainty is too large), policy decisions based on the ensemble average struggle to obtain credibility from the general public [15]. In fact, the large inter-model variability of future rainfall extremes has been proven through many studies [6,7]. Therefore, the reliability of future rainfall extremes projection at a specific site is generally low. Therefore, we need a new approach that can convince the general public about the uncertain future rainfall extremes.
We thought that a new approach to presenting convincing future rainfall extremes to the general public is to reduce the gap between future ensemble members. From the perspective of reducing the differences arising from various climate model combinations, the comparison of the results of Figures 4a and 6 suggests that using future SAT or DPT in conjunction with a non-stationary model is better than using future AMR values. Consistency of future projections will contribute to future rainfall extremes gaining public trust and being reflected in policy.
The total uncertainty of future rainfall extremes decomposed into major sources is plotted in Figure 7. The total uncertainty was expressed as the coefficient of variation of samples extracted from the posterior distribution of future rainfall extremes. When the stationary model was applied, the uncertainty due to parameter estimation and the uncertainty due to various climate model combinations were similar. In fact, the uncertainty resulting from the estimation of the parameter of probability distribution is directly related to the data period we have observed. It is very difficult to reduce the uncertainty arising from parameter estimation since it is actually difficult to artificially increase the data period. The uncertainty of these two sources, parameter estimates and climate model combinations, can be found to propagate, without loss, to estimates of future rainfall quantiles. On the other hand, it can be found that most of the uncertainty included in the future rainfall quantile estimated from the co-variate-based non-stationary model is due to the parameter estimation of the GEV distribution. This result is because various combinations of climate models simulate similar future SATs or DPTs. In other words, this means that even if an ensemble of various climate model combinations is applied, the uncertainty of rainfall extremes does not increase significantly.

Uncertainty of Rate of Change
In this section, as an application of the proposed approach, we tried to quantify the rate of change of rainfall extremes for future temperature increase and the degree of uncertainty in such a rate of change.
By combining the non-stationary model n 1 and Equation (10) or the non-stationary model n 2 and Equation (11), an ensemble of rainfall quantiles corresponding to the increase in SAT or DPT can be obtained. Figure 8 shows the probability that the rainfall quantile for the return level of 50 years exceeds a certain rate of change (i.e., likelihood of increase, LoI) given that SAT or DPT increase (∆SAT or ∆DPT). For example, in the scenario in which the SAT rises by 6 • C at the Chuncheon site, the probability that the rainfall quantile for the 50-year return level increases by 20% or more is about 40%. Similarly, in the scenario where the DPT rises by 4 • C at the Cheonan site, it can be found that the probability of a rate of change of 1.1 or higher is about 50%.

tremes.
We thought that a new approach to presenting convincing future rainfall extremes to the general public is to reduce the gap between future ensemble members. From the perspective of reducing the differences arising from various climate model combinations, the comparison of the results of Figure 4a and Figure 6 suggests that using future SAT or DPT in conjunction with a non-stationary model is better than using future AMR values. Consistency of future projections will contribute to future rainfall extremes gaining public trust and being reflected in policy. The total uncertainty of future rainfall extremes decomposed into major sources is plotted in Figure 7. The total uncertainty was expressed as the coefficient of variation of samples extracted from the posterior distribution of future rainfall extremes. When the stationary model was applied, the uncertainty due to parameter estimation and the uncertainty due to various climate model combinations were similar. In fact, the uncertainty resulting from the estimation of the parameter of probability distribution is directly related to the data period we have observed. It is very difficult to reduce the uncertainty arising from parameter estimation since it is actually difficult to artificially increase the data period. The uncertainty of these two sources, parameter estimates and climate model combinations, can be found to propagate, without loss, to estimates of future rainfall quantiles. On the other hand, it can be found that most of the uncertainty included in the future rainfall quantile estimated from the co-variate-based non-stationary model is due to the parameter estimation of the GEV distribution. This result is because various combinations of climate models simulate similar future SATs or DPTs. In other words, this means that even if an ensemble of various climate model combinations is applied, the uncertainty of rainfall extremes does not increase significantly.

Uncertainty of Rate of Change
In this section, as an application of the proposed approach, we tried to quantify the rate of change of rainfall extremes for future temperature increase and the degree of uncertainty in such a rate of change.
By combining the non-stationary model 1 and Equation (10) or the non-stationary model 2 and Equation (11), an ensemble of rainfall quantiles corresponding to the increase in SAT or DPT can be obtained. Figure 8 shows the probability that the rainfall If Figure 8 is substituted for a specific DPT rise scenario, the cumulative probability distribution for the rate of change of the rainfall quantile can be obtained as follows. The LoI(r) in Figure 8 is the probability that the rate of change of the rainfall quantile for a specific return level is more than r under a specific DPT rising condition, so the probability that the rate of change is less than r is 1 − LoI(r). That is, the cumulative probability distribution of the rate of change r is 1 − LoI(r). Figure 9 shows the cumulative distribution of the rate of change of the rainfall quantile for the 50-year return level in the scenario where DPT rises by 1, 3, and 6 °C respectively. At the Cheonan site, the ensemble average of the rate of change of the rainfall quantile for the 50-year return level was 1.11, and the 95% confidence interval of the rate of change was calculated from 0.66 to 1.76 when the DPT increased by 3 °C .

Conclusions
In most studies, the approach to investigate future rainfall extremes takes a method of considering future rainfall data obtained from GCMs or RCMs as future observation data. However, the very large deviations between the statistics of rainfall extremes obtained using future rainfall time series produced from various climate models create a lot If Figure 8 is substituted for a specific DPT rise scenario, the cumulative probability distribution for the rate of change of the rainfall quantile can be obtained as follows. The LoI(r) in Figure 8 is the probability that the rate of change of the rainfall quantile for a specific return level is more than r under a specific DPT rising condition, so the probability that the rate of change is less than r is 1 − LoI(r). That is, the cumulative probability distribution of the rate of change r is 1 − LoI(r). Figure 9 shows the cumulative distribution of the rate of change of the rainfall quantile for the 50-year return level in the scenario where DPT rises by 1, 3, and 6 • C respectively. At the Cheonan site, the ensemble average of the rate of change of the rainfall quantile for the 50-year return level was 1.11, and the 95% confidence interval of the rate of change was calculated from 0.66 to 1.76 when the DPT increased by 3 • C. If Figure 8 is substituted for a specific DPT rise scenario, the cumulative probability distribution for the rate of change of the rainfall quantile can be obtained as follows. The LoI(r) in Figure 8 is the probability that the rate of change of the rainfall quantile for a specific return level is more than r under a specific DPT rising condition, so the probability that the rate of change is less than r is 1 − LoI(r). That is, the cumulative probability distribution of the rate of change r is 1 − LoI(r). Figure 9 shows the cumulative distribution of the rate of change of the rainfall quantile for the 50-year return level in the scenario where DPT rises by 1, 3, and 6 °C respectively. At the Cheonan site, the ensemble average of the rate of change of the rainfall quantile for the 50-year return level was 1.11, and the 95% confidence interval of the rate of change was calculated from 0.66 to 1.76 when the DPT increased by 3 °C .

Conclusions
In most studies, the approach to investigate future rainfall extremes takes a method of considering future rainfall data obtained from GCMs or RCMs as future observation data. However, the very large deviations between the statistics of rainfall extremes obtained using future rainfall time series produced from various climate models create a lot

Conclusions
In most studies, the approach to investigate future rainfall extremes takes a method of considering future rainfall data obtained from GCMs or RCMs as future observation data. However, the very large deviations between the statistics of rainfall extremes obtained using future rainfall time series produced from various climate models create a lot of confusion in establishing climate change adaptation measures. The uncertainty of future rainfall extremes hinders the public's confidence in the outputs produced by many experts. In this study, we proposed an approach that integrates future SAT and DPT information and a non-stationary GEV model to obtain a future rainfall extremes ensemble. We also produced outputs of future rainfall extremes that reflect future SAT and DPT information simulated in various climate model combinations. The future rainfall extremes estimated from the future rainfall data differed greatly depending on the applied climate model combination, but the future rainfall extremes estimated using the proposed approach were relatively consistent regardless of the climate model combination.
The non-stationarity of AMR was implemented by expressing the scale parameter of the GEV distribution as a function of the SAT or DPT observed on the day of the AMR. We were able to find a non-stationary GEV distribution that gave a smaller value than the nllh of the stationary GEV distribution. However, since the difference in nllh is likely to change anyway during the parameter estimation process, it was recognized that the performances of the stationary GEV distribution and the non-stationary GEV distribution in terms of the goodness of fit of the data were almost similar. In this study, the rainfall quantile for various parameter combinations was generated using MCMC sampling from the posterior distribution of parameters derived by the MH algorithm. From the point of view that the model with less uncertainty inherent in the rainfall quantile, which is the result of frequency analysis, is more reliable, we found that performing an SAT-or DPT-based non-stationary frequency analysis instead of a stationary frequency analysis is advantageous in obtaining a more robust rainfall quantile.
In this study, a method of estimating the reference value of the co-variate corresponding to the return level based on the observed AMR was proposed, but the relationship between the reference value of the co-variate and the return level is not yet clear. Furthermore, the uncertainty of the reference co-variate needs to be further studied. The issue of reference co-variate will be an important research topic in non-stationary frequency analysis.
Finally, the uncertainty of the rate of change of the rainfall quantile in the future SAT or DPT conditions could be quantified by using the rainfall quantile ensemble in the present and future SAT or DPT conditions that can be obtained during the uncertainty analysis process.
The approach we propose in this study is not a new theory in the domain of nonstationary frequency analysis. In this study, even though there are still many limitations, a method to reduce the uncertainty of non-stationary frequency analysis was discussed by introducing a reference value to the co-variate of the non-stationary frequency analysis. In addition, when projecting the quantile of future rainfall extremes, it was suggested that using future co-variate data could improve the consistency between climate models rather than using future rainfall data simulated in the climate model. Finally, using the proposed method, the uncertainty of the rate of change of future rainfall extremes could be identified by dividing the uncertainty based on the parameter estimation of the nonstationary model and the uncertainty due to the inter-model variability of various climate model combinations. These points will be able to provide additional insights based on current knowledge.  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: [www.weather.go.kr] (last accessed on 6 February 2020).