Improvement of Extreme Value Modeling for Extreme Rainfall Using Large-Scale Climate Modes and Considering Model Uncertainty

Abstract

Extreme value modeling for extreme rainfall is one of the most important processes in the field of hydrology. For the improvement of extreme value modeling and its physical meaning, large-scale climate modes have been widely used as covariates of distribution parameters, as they can physically account for climate variability. This study proposes a novel procedure for extreme value modeling of rainfall based on the significant relationship between the long-term trend of the annual maximum (AM) daily rainfall and large-scale climate indices. This procedure is characterized by two main steps: (a) identifying significant seasonal climate indices (SCIs), which impact the long-term trend of AM daily rainfall using statistical approaches, such as ensemble empirical mode decomposition, and (b) selecting an appropriate generalized extreme value (GEV) distribution among the stationary GEV and nonstationary GEV (NS-GEV) using time and SCIs as covariates by comparing their model fit and uncertainty. Our findings showed significant relationships between the long-term trend of AM daily rainfall over South Korea and SCIs (i.e., the Atlantic Meridional Mode, Atlantic Multidecadal Oscillation in the fall season, and North Atlantic Oscillation in the summer season). In addition, we proposed a model selection procedure considering both the Akaike information criterion and residual bootstrap method to select an appropriate GEV distribution among a total of 59 GEV candidates. As a result, the NS-GEV with SCI covariates generally showed the best performance over South Korea. We expect that this study can contribute to estimating more reliable extreme rainfall quantiles using climate covariates.

1. Introduction

Extreme value modeling for hydrological extreme events (e.g., heavy rainfall and floods) is an essential task for the design of hydraulic structures. For extreme value modeling, the block maxima approach has been widely used, and it uses a sequence of maximum observations extracted from equal periods, such as annual maximum (AM) daily rainfall. The set of block maxima is assumed to be independent and identically distributed (iid) and is fitted with a probability distribution model, such as a generalized extreme value (GEV) distribution, to estimate the design quantiles of extreme hydrological events.

As the impact of climate change has become a significant issue, many efforts have been made to consider nonstationarity in hydrologic applications. One popular approach applies a variety of candidate nonstationary models to the nonstationary data and selects an appropriate model based on model diagnostics. It generally employs the maximum likelihood estimation (MLE) method to estimate nonstationary model parameters due to its adaptability to changes in model structures [1]. So far, this approach has been studied extensively and can be referred to as a “user-friendly” method. However, there is an issue regarding ergodicity in nonstationary extreme value modeling from a statistical point of view [2,3,4]. More specifically, statistical properties from temporal statistics would involve the assumption of ergodicity since AM time series can be theoretically considered as a stochastic process. However, if the observed data has a trend or is affected by external variables (i.e., the process is nonstationary), the ergodicity cannot hold and inductive inference from the data would not be possible [3,4,5]. In short, it is necessary to determine the relationships between model parameters and covariates prior to model parameter estimation to deal with the ergodicity issue. This is an important issue that needs to be addressed in terms of statistical hydrology, but hydrological applications mostly accept nonstationary models, as we mentioned earlier. Hence, we apply the standardized approach in this study, and the limitations regarding ergodicity will be discussed in the discussion section later.

For nonstationary extreme value modeling, the nonstationary GEV (NS-GEV) is a representative model, and it employs time or exogenous covariates (e.g., large-scale climate modes and hydrological variables) [6,7,8,9,10,11,12,13,14,15]. Here, time t is commonly used as a covariate for modeling the linear or polynomial trend of GEV parameters. It is a simple and straightforward way to apply nonstationarity to the probability distribution model. However, considering only time t as a covariate could lead to problems such as an increase in uncertainty of quantile estimates and distortion of the probability distribution model according to extrapolation [16,17]. In addition, the trends of hydrological data can vary in the short or long term because of climate variability and other external forces [18,19,20]. Thus, large-scale climate modes have been employed as potential climate covariates since the mid-2000s as they can consider natural climate variability.

In general, extreme value modeling for rainfall with climate covariates has been performed as follows: (1) identifying significant climate indices related to extreme rainfall; and (2) selecting an appropriate probability distribution model among the various GEV models based on the statistical evaluation of the model fit. As such, various forms of NS-GEV can be modeled to estimate the magnitude and occurrence frequency of extreme rainfall events using climate indices as covariates [8,9,11,12,14,21]. Then, information criteria such as Akaike information criterion (AIC) and Bayesian information criterion (BIC) are generally used to select the most appropriate model. For example, Vasiliades et al. [8] considered the four climate indices (i.e., the Southern Oscillation Index (SOI), North Atlantic Oscillation (NAO), Pacific Decadal Oscillation (PDO), and Mediterranean Oscillation Index (MOI)) reported to illustrate the dependence of the Mediterranean climate on large-scale low-frequency climate patterns. These indices were employed as covariates for the GEV model with a conditional density network, and an appropriate GEV model was finally selected based on the corrected AIC (AICc) and BIC. Agilan and Umamahesh [9] conducted extreme value modeling for developing a nonstationary rainfall intensity-duration-frequency curve. They used five covariates, including El-Nino-Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD), which are known to impact extreme rainfall over India; then, they selected an appropriate model using the AICc. They concluded that the global processes (i.e., global warming, ENSO cycle, and IOD) are the best covariates for the long-duration extreme rainfall of Hyderabad city, India.

As shown in the abovementioned studies, many previous works considered several climate indices reported in the literature for modeling NS-GEV of rainfall extremes. The climate indices were directly used as covariates of location and scale parameters, which are closely linked to the mean and variance of the data, respectively. It could be possible to model a more reasonable NS-GEV by considering the connection between climate indices and the trend inherent in extreme rainfall. For this, a preliminary analysis is necessary to determine the trends in the mean and variance of AM rainfall and to identify the impact of climate indices on these trends. Recently, many researchers have employed modern signal processing known as ensemble empirical mode decomposition (EEMD) to extract long-term trends from a given data series. The EEMD has the advantage of detecting long-term trends after extracting oscillatory patterns from an original data series [22,23]. It is also successfully applied to hydro-climatology variables [24,25,26,27]. For example, Kim et al. [27] identified that the behavior of the Atlantic Multidecadal Oscillation (AMO) index indicates a long-term trend in the monthly precipitation series of South Korea through a statistical analysis procedure along with EEMD. Chen et al. [26] employed the EEMD to explore the trend of the AM daily precipitation data. They extracted long-term trends in the mean and variance of the AM daily precipitation data and successfully identified that the extracted trends can be applied to extreme value modeling.

The increase of uncertainty in design levels is one of the most important issues in nonstationary extreme value modeling. The use of a more complex model (i.e., nonstationary distribution model) generally allows a good fit for the given samples, but this could provide quantile estimates with a large amount of uncertainty [3]. In general, however, the best model is selected using only the information criteria that assess the model performance based on the maximized log-likelihood (i.e., model fit) [26,28,29,30,31,32,33], so that the selected nonstationary model usually yields unreliable design quantiles [16]. Cooley [34] demonstrated that relying solely on the information criteria is not enough to select an appropriate model because they cannot take into account uncertainty. Nevertheless, most of the previous studies assessed uncertainties in model parameters and design quantiles after selecting the best model based on the information criteria [13,16,17,21,35,36,37,38,39,40,41]. Therefore, to provide more reliable design quantiles, the uncertainty of the candidate models should be assessed prior to selecting an appropriate model. The bootstrap method is generally recommended to measure the uncertainty of quantile estimates [36,37,38,42]. It is computationally efficient and provides realistic asymmetric confidence intervals (CIs) without asymptotic assumptions [3,10].

Many studies have attempted extreme value modeling by applying the NS-GEV with climate covariates to estimate rainfall quantiles considering climate variability. However, most of these are limited to modeling for several rainfall gauging stations [8,9,11,12,14,36], and few studies propose an overall procedure from the identification of appropriate climate indices impacting regional rainfall extremes to the selection of an appropriate nonstationary distribution model with climate covariates (e.g., India [21], southern France [38], and southern U.S. [43]). To the best of our knowledge, there has been less attention to considering the uncertainty in the model selection procedure. The main objective of this study is to propose a procedure for extreme value modeling with large-scale climate modes and to consider the uncertainty of model selection. It can be divided into two main subsections: (a) identifying significant seasonal climate indices (SCIs) that impact the long-term trends of AM daily rainfall using the EEMD, and (b) selecting an appropriate GEV distribution considering both model fit and uncertainty among the stationary GEV (ST-GEV) and various NS-GEVs using time and SCIs as covariates. The EEMD was applied to AM daily rainfall observed at 61 stations over South Korea, and the residue was extracted that represent the long-term trend of AM daily rainfall. By conducting correlation analysis between the extracted residues and the various climate indices, significant SCIs that impact the trend of AM daily rainfall were selected. Extreme value modeling was then performed using the significant climate indices as covariates of GEV parameters. Considering both model fit and uncertainty, an appropriate GEV distribution was finally selected at each station. We also discussed the physical meaning of the significant climate indices selected over South Korea and the feasibility of the procedure for nonstationary extreme value modeling.

2. Data

2.1. Annual Maximum Daily Rainfall Data

The spatial domain of the study is South Korea, East Asia. Rainfall over South Korea is mainly affected by the subtropical East Asian summer monsoon and typhoons. Therefore, approximately 70% of the annual rainfall is concentrated in the rainy season from June to September [44,45,46]. We extracted annual maximum (AM) daily rainfall data observed at 61 stations over South Korea. The data recorded until 2018 cover a minimum of 40 years with no missing year. Figure 1 illustrates the study area and the location of the rainfall stations, and associated information is presented in

Table 1. Locations and record lengths of the rainfall gauging stations in this study.

Code	Name	Longitude	Latitude	Record Length (year)	Code	Name	Longitude	Latitude	Record Length (year)
90	Sokcho	128.6	38.3	50	201	Ganghwa	126.4	37.7	45
100	Daegwallyeong	128.7	37.7	46	202	Yangpyeong	127.5	37.5	45
101	Chuncheon	127.7	37.9	52	203	Icheon	127.5	37.3	45
105	Gangneung	128.9	37.8	57	211	Inje	128.2	38.1	45
108	Seoul	127.0	37.6	57	212	Hongcheon	127.9	37.7	45
112	Incheon	126.6	37.5	57	221	Jecheon	128.2	37.2	45
114	Wonju	127.9	37.3	45	226	Boeun	127.7	36.5	45
115	Ulleungdo	130.9	37.5	57	232	Cheonan	127.1	36.8	45
119	Suwon	127.0	37.3	54	235	Boryeong	126.6	36.3	45
127	Chungju	128.0	37.0	45	236	Buyeo	126.9	36.3	45
129	Seosan	126.5	36.8	50	238	Geumsan	127.5	36.1	45
130	Uljin	129.4	37.0	46	243	Buan	126.7	35.7	45
131	Cheongju	127.4	36.6	51	244	Imsil	127.3	35.6	45
133	Daejeon	127.4	36.4	49	245	Jeongeup	126.9	35.6	45
135	Chupungnyeong	128.0	36.2	57	247	Namwon	127.3	35.4	45
138	Pohang	129.4	36.0	57	260	Jangheung	126.9	34.7	45
140	Gunsan	126.8	36.0	50	261	Haenam	126.6	34.6	45
143	Daegu	128.6	35.9	57	262	Goheung	127.3	34.6	45
146	Jeonju	127.2	35.8	57	272	Yeongju	128.5	36.9	45
152	Ulsan	129.3	35.6	57	273	Mungyeong	128.1	36.6	45
156	Gwangju	126.9	35.2	57	277	Yeongdeok	129.4	36.5	45
159	Busan	129.0	35.1	57	278	Uiseong	128.7	36.4	45
162	Tongyeong	128.4	34.8	50	279	Gumi	128.3	36.1	45
165	Mokpo	126.4	34.8	57	281	Yeongcheon	129.0	36.0	45
168	Yoesu	127.7	34.7	57	284	Geochang	127.9	35.7	45
170	Wando	126.7	34.4	45	285	Hapcheon	128.2	35.6	45
174	Suncheon	127.4	35.0	45	288	Miryang	128.7	35.5	45
184	Jeju	126.5	33.5	57	289	Sancheong	127.9	35.4	45
188	Seongsan	126.9	33.4	45	294	Geoje	128.6	34.9	45
189	Seogwipo	126.6	33.2	57	295	Namhae	127.9	34.8	45
192	Jinju	128.0	35.2	49

. The monthly rate of the AM daily rainfall occurrences over South Korea was identified, as depicted in Figure 2. Most of the AM daily rainfall occurred between June and September, particularly in July and August.

2.2. Climate Indices

Climate indices are quantitative indicators representing large-scale climate modes in the ocean-atmospheric system. They can be characterized by variations in the sea surface temperature, sea level pressure, atmospheric circulation, and so on. Each climate index is derived through various measurements and is widely used to represent regional climatic characteristics [27]. In the current study, we used a total of 12 climate indices that can cover the observed rainfall data periods: Atlantic Meridional Mode (AMM), AMO, Arctic Oscillation (AO), NAO, Niño 1 + 2 (NINO12), Niño 3.4 (NINO34), Niño 4 (NINO4), North Pacific pattern (NP), PDO, Pacific Meridional Mode (PMM), Pacific-North American pattern (PNA), and SOI. These climate indices were obtained from the Physical Sciences Division (PSD) of the National Oceanic and Atmospheric Administration Earth System Research Laboratory (NOAA/ESRL) (www.esrl.noaa.gov/psd/data/climateindices/list (accessed on 7 June 2019)).

3. Methodology

This study is grouped largely into two main steps to conduct extreme value modeling with large-scale climate modes: Step 1. identifying significant climate indices that impact the trend of AM daily rainfall based on statistical approaches; Step 2. selecting an appropriate GEV distribution by comparing the performance of various ST- and NS-GEV with climate covariates. The procedure can be subdivided into five steps. Figure 3 presents a brief overview of this study, and the details of the methodology are described in the following subsections.

3.1. Seasonal Climate Indices

All climate indices are provided on a monthly scale. To employ the SCIs as covariates in nonstationary extreme value modeling, we converted the monthly climate indices to seasonal ones by seasonal averaging (i.e., Spring: March-April-May (MAM); Summer: June-July-August (JJA); Fall: September-October-November (SON); and Winter: December-January-February (DJF)). To consider the time lag for using SCIs as a predictor, the SCIs observed prior to the occurrence of AM daily rainfall should be used. As most AM daily rainfall events over South Korea occur in the JJA season (see Figure 2), the four seasons of lagged SCIs before the JJA season are considered as described in

Table 2. Four seasons of lagged SCIs were considered in this study.

Abbreviation	JJA(−1)	SON(−1)	DJF(−1)	MAM
Description	Averaged from Jun. to Aug. (summer season) in the previous year of AM rainfall occurrence	Averaged from Sep. to Nov. (fall season) in the previous year of AM rainfall occurrence	Averaged from Dec. to Feb. (winter season) in the previous year of AM rainfall occurrence	Averaged from Mar. to May (spring season) in the year of AM rainfall occurrence

.

3.2. Ensemble Empirical Mode Decomposition

EEMD is a decomposition method that has been used widely in hydrological time series. By decomposing an original time series into a set of intrinsic mode functions (IMFs), the EEMD can subsequently extract a long-term trend inherent in the time series. The IMFs indicate physically meaningful information, and they should satisfy two conditions: (1) the number of extrema and zero crossings must either be equal to each other or differ at most by 1 in the original time series, and (2) the mean of the upper and lower envelopes, which are defined by connecting all local minima and maxima, should be zero at any point [47]. A monotonic function that remains after the decomposition is a residue, which indicates a long-term trend. Here, we briefly describe the EEMD procedure with an example of a time series, as shown in Figure 4. The red-dotted and blue-dotted lines are the upper and lower envelopes of the original time series, respectively, and the gray-dotted line is the mean value of the local maxima and local minima. Assuming that there is an original time series, we can obtain N number of IMFs as follows: (1) identify the upper and lower envelopes by connecting the local maxima and minima $(y_u\left(t\right), y_l\left(t\right))$ using a cubic spline; (2) calculate the mean value of the local maxima and local minima as $y_m\left(t\right)=\left(y_u\left(t\right)+y_l\left(t\right)\right)/2$; (3) extract the mean value from the original time series as $h\left(t\right)=y\left(t\right)-y_m\left(t\right)$; (4) repeat steps (1) to (3) until $h\left(t\right)$ satisfies the condition of the IMF; (5) define a new time series by extracting $h\left(t\right)$ from $y\left(t\right)$ and repeat steps (1) to (5) until no more IMFs can be extracted. Finally, $y\left(t\right)$ is composed of the sum of the IMFs and a residue $\left(res\left(t\right)\right)$, as presented in Equation (1):

$$y\left(t\right)={{\displaystyle \sum }}_{i=1}^{N}IM{F}_i\left(t\right)+res\left(t\right)$$

(1)

where $N$ is the number of IMFs.

In this study, the EEMD is applied to AM daily rainfall to extract the long-term trend of AM daily rainfall as follows (for additional details, see Chen et al. [26]): (1) define a residue $re{s}_m\left(t\right)$, extracted from the original time series $y\left(t\right)$ using the EEMD procedure, as a long-term trend in the mean of AM daily rainfall (hereafter referred to as the mean trend); (2) calculate the time series of the variance as $var\left(t\right)={\left(y\left(t\right)-re{s}_m\left(t\right)\right)}^2$; (3) define a residue $re{s}_v\left(t\right)$, extracted from the time series of the variance $var\left(t\right)$ using the EEMD procedure, as a long-term trend in the variance of AM daily rainfall (hereafter referred to as the variance trend).

3.3. Spearman’s Rank Correlation Analysis

Spearman’s rank correlation is one of the most widely used statistical estimators to measure the statistical dependence between two different hydrometeorological series [14,38,48,49]. It is based on the nonparametric rank correlation and describes the statistical relationship using a monotonic function. Spearman’s rank correlation coefficient ($r_s$) between two time series $X_1\left(t\right)$ and $X_2\left(t\right)$ is defined by Equation (2):

$$r_s=1-\frac{6{{\displaystyle \sum }}_{t=1}^n{d_t}^2}{n\left(n^2-1\right)}$$

(2)

where $d_t=rank\left(X_1\left(t\right)\right)-rank\left(X_2\left(t\right)\right)$ is the difference between the two ranks of each observation at time $t$, and $n$ is the number of observations. If the two time series have a perfect monotonic function, the value of $r_s$ is close to −1 or 1. As a rule of thumb, a value of $r_s$ between 0.7 and 1.0 (−0.7 and −1.0) represents a strong positive (negative) correlation, a value between 0.3 and 0.7 (−0.3 and −0.7) represents a moderate positive (negative) correlation, and a value between 0 and 0.3 (0 and −0.3) represents a weak positive (negative) correlation [50,51,52].

To identify the significant SCIs impacting AM daily rainfall over South Korea, $r_s$ between the SCI and the mean and variance trends of AM daily rainfall are calculated for all employed stations. Then, the percentage of stations with a significant $r_s$ value ($P_{sig}$) at the 1% significance level is calculated for each SCI.

3.4. Generalized Extreme Value Distribution Modeling

The GEV distribution is widely used for estimating the magnitude and occurrence probability of hydrological extreme events. Let $X=x_1, x_2, \cdots , x_n$ be the AM time series data with sample size $n$, which are assumed to be independent and identically distributed (iid). The cumulative distribution function $\left(F\left(x\right)\right)$ of the ST-GEV for $X$ is expressed by Equation (3):

$$F\left(x\right)=\exp \left[-{\left\{1+\xi \frac{x-\mu }{\sigma }\right\}}^{-1/\xi }\right]$$

(3)

where $\mu$ is the location parameter related to the mean of the data, $\sigma$ (>0) is the scale parameter related to the variability of the data, and $\xi$ is the shape parameter related to the heaviness of the distribution tail.

When a trend or the impact of external variables in the observed data is considered in extreme value modeling, the distribution parameters can be modeled as functions of covariates, such as time or climate indices. In this study, the location and scale parameters of the GEV distribution were defined as functions of time-dependent covariates as follows:

$$\mu \left(t\right)={\mu }_0+{\mu }_{1}co{v}_1\left(t\right)+{\mu }_{2}co{v}_2\left(t\right)+\cdots +{\mu }_{p}co{v}_p\left(t\right)$$

(4)

$$\sigma \left(t\right)=\exp \left({\sigma }_0+{\sigma }_{1}co{v}_1\left(t\right)+{\sigma }_{2}co{v}_2\left(t\right)+\cdots +{\sigma }_{p}co{v}_p\left(t\right)\right)$$

(5)

where $p$ is the number of covariates and $co{v}_p\left(t\right)$ is the $p$-th covariate at time $t$. As the scale parameter should be positive for all $t$, an exponential function was employed as a link function for the scale parameter. The shape parameter was assumed to be stationary, as it is difficult to estimate the shape parameter with precision [1,10].

The combinations of significant SCIs are used as covariates of location and scale parameters. Depending on the nonstationarity in the GEV parameters, two types of NS-GEV are defined: NS-GEV(1,0,0), in which only a location parameter is assumed to be a function of covariates, and NS-GEV(1,1,0), in which both location and scale parameters are assumed to be functions of covariates. The parameters of GEV distributions are estimated using the MLE method.

3.5. Appropriate Model Selection Considering Uncertainty

The Akaike information criterion (AIC) is an information-theoretic model selection method based on Kullback–Leibler information loss [53]. Since the AIC evaluates model performance based on the maximized log-likelihood and the number of parameters for the fitted distribution, a distribution model with a good fit and parameter parsimony was selected as an appropriate model. The equation of the AIC is expressed as

$$AIC=-2\log \left(ML\right)+2k$$

(6)

where $\log \left(ML\right)$ is the maximized log-likelihood of the fitted distribution and $k$ is the number of distribution parameters. By comparing the AIC value of candidate distribution models, the model with the smallest AIC value is selected as an appropriate model, considering both the parameter parsimony and the goodness-of-fit. The AIC of each model can be rescaled as follows:

$$rAIC=AIC-\min \left(AIC\right)$$

(7)

where $rAIC$ is the rescaled $AIC$ and $\min \left(AIC\right)$ is the smallest $AIC$ among all candidate models. The models with $rAIC\le 2$ can also be reasonable choices [53].

The bootstrap method is generally recommended to measure the uncertainty of the quantile estimate in a nonstationary distribution model because it is computationally efficient and provides realistic asymmetric confidence intervals (CIs) without asymptotic assumptions [3,10]. For the NS-GEV, the residual bootstrap method has been employed to calculate the confidence intervals for the parameter and quantile estimates [9,14,54,55,56]. The residual bootstrap method can be conducted as follows:

• For the fitted GEV distribution, transform the AM time series data ($X=x_1, x_2, \cdots , x_n$) into the standardized residuals with no trend ($\tilde{X}={\tilde{x}}_1, {\tilde{x}}_2, \cdots , {\tilde{x}}_n$) as follows [1]:

  $$\widetilde{x_t}=\frac{1}{\widehat{\xi }}\ln \left[1+\widehat{\xi }\left(\frac{x_t-\widehat{\mu }\left(t\right)}{\widehat{\sigma }\left(t\right)}\right)\right],\hspace{1em}\hspace{1em}t=1, 2,\cdots ,n$$

  (8)

  where the $\widehat{\mu }\left(t\right)$, $\widehat{\sigma }\left(t\right)$, and $\widehat{\xi }$ are the location, scale, and shape parameters of the fitted GEV distribution, respectively.
• Obtain a new sample of $\tilde{X}$ by resampling residuals with replacement and back-transforming the resampled residuals using Equation (8).
• For back-transformed samples, estimate the T-year quantile at each time $t$ (${\widehat{Q}}_T\left(t\right)$, $t=1, 2,\cdots ,n$) using the same GEV distribution.
• Repeat steps (2)–(3) $N$ times and calculate the time-averaged 95% CIs for the T-year quantiles ($C{I}_{Q_T}^{95\%}$) as follows:

  $$C{I}_{Q_T}^{95\%}=\frac{1}{n}{{\displaystyle \sum }}_{t=1}^n\left({\widehat{Q}}_T^{97.5\%}\left(t\right)-{\widehat{Q}}_T^{2.5\%}\left(t\right)\right)$$

  (9)

  where ${\widehat{Q}}_T^{97.5\%}\left(t\right)$ and ${\widehat{Q}}_T^{2.5\%}\left(t\right)$ are the 97.5th and 2.5th percentiles of $N$ ordered samples of ${\widehat{Q}}_T\left(t\right)$.

In this study, the model selection procedure combining the AIC and residual bootstrap method is performed to select an appropriate GEV model considering both model fit and uncertainty. First, the AIC values of all GEV candidates are calculated, and the GEV distributions with $rAIC\le 2$ are selected. Second, for each selected GEV distribution, the residual bootstrap method is repeated 1000 times to calculate the time-averaged 95% CIs for the 100-year quantiles ($C{I}_{Q_{100}}^{95\%}$). Finally, the GEV distribution with the smallest $C{I}_{Q_{100}}^{95\%}$ is selected as an appropriate distribution model.

4. Results

4.1. Mean and Variance Trends of AM Daily Rainfall Extracted Using EEMD

At first, the residue of AM daily rainfall is extracted using EEMD. Figure 5 presents the time series data and extracted residues of AM daily rainfall and its variance at two representative stations (Ulleungdo and Daegu) for illustration purposes. On the left side, the black line indicates the observed AM daily rainfall at the station, and the red line represents a residue $\left(re{s}_m\left(t\right)\right)$ indicating the long-term mean trend of AM daily rainfall. On the right side, the black line is the time series of the variance, and the red line is a residue $(re{s}_v\left(t\right))$ indicating the long-term variance trend of AM daily rainfall. Chen et al. [26] modeled these residues as linear or quadratic functions of time and directly used them as time-dependent parameters of nonstationary probability distribution models. In this study, the extracted residues were used to find the significant climate indices that impact the long-term trend of AM daily rainfall.

4.2. Significant Seasonal Climate Indices for AM Daily Rainfall over South Korea

Table 3. Percentage of stations with a significant $r_s$ value between the lagged SCI and $re{s}_m\left(t\right)$ at the 1% significance level.

Climate Index	Season
	JJA(−1)	SON(−1)	DJF(−1)	MAM
AMM	-	77%	67.2%	-
AMO	75.4%	77%	77%	77%
AO	-	-	-	-
NAO	57.4%	-	13.1%	-
NINO12	4.9%	-	-	-
NINO4	-	-	-	-
NINO34	19.7%	-	-	-
NP	-	-	-	6.6%
PDO	9.8%	-	-	9.8%
PMM	-	-	-	-
PNA	-	-	-	8.2%
SOI	-	-	-	-

presents the percentage of stations with a significant $r_s$ value ($P_{sig}$) between the lagged SCIs and $re{s}_m\left(t\right)$. A $P_{sig}$ greater than 50% is highlighted in bold. Among the 12 climate indices, the AMM, AMO, and NAO indices include SCIs with more than 50% $P_{sig}$.

Table 4. Percentage of stations with a significant $r_s$ value between the lagged SCI and $re{s}_v\left(t\right)$ at the 1% significance level.

Climate Index	Season
	JJA(−1)	SON(−1)	DJF(−1)	MAM
AMM	-	75.4%	68.9%	-
AMO	72.1%	82%	80.3%	73.8%
AO	-	-	-	-
NAO	68.9%	-	13.1%	-
NINO12	3.3%	-	-	-
NINO4	-	-	-	-
NINO34	14.8%	-	-	-
NP	-	-	-	3.3%
PDO	13.1%	-	-	11.5%
PMM	-	-	-	-
PNA	-	-	-	8.2%
SOI	-	-	-	-

presents the $P_{sig}$ between the lagged SCIs and $re{s}_v\left(t\right)$. Similar to the results presented in Table 3, the AMM, AMO, and NAO indices include SCIs with more than 50% $P_{sig}$.

The $r_s$ values between the mean trend of AM daily rainfall and (a) AMM, (b) AMO, (c) NAO indices are depicted in Figure 6. The seasonal AMM indices indicate positive $r_s$ values in the overall inland area, whereas they indicate negative $r_s$ values in some eastern and southern coastal areas of the Korean Peninsula. For each season, the average absolute $r_s$ value for all stations ($\overline{\left|r_s\right|}$) was calculated. Among the four seasonal AMM indices, the AMM index for the SON(−1) season, AMM_SON(−1), indicates the strongest correlation ($\overline{\left|r_s\right|}$ = 0.55) with the mean trend of AM daily rainfall over South Korea in Figure 6a. The results of the AMO index are very similar to those of the AMM_SON(−1). For all seasons, the AMO index indicates a strong correlation with the mean trend of AM daily rainfall in Figure 6b. Among them, AMO_SON(−1) indicates the strongest correlation ($\overline{\left|r_s\right|}$ = 0.59). In contrast to the results of the AMO indices, the NAO index demonstrates different patterns according to the seasons. While the NAO_JJA(−1) and NAO_SON(−1) present overall negative correlations with the mean trend of AM daily rainfall, the NAO_DJF(−1) and NAO_MAM indicate very weak positive correlations in Figure 6c. In particular, NAO_JJA(−1) presents the strongest correlation ($\overline{\left|r_s\right|}$ = 0.34), and its spatial pattern is opposite to those of AMM_SON(−1) and AMO_SON(−1).

The $r_s$ values between the variance trend of AM daily rainfall and the (a) AMM, (b) AMO, and (c) NAO indices are depicted in Figure 7. The spatial pattern of $r_s$ values of the variance trend is similar to that of the mean trend, except that the $r_s$ values in some inland areas (lying between latitude 35.5° N–37.5° N and longitude 126.5° E–129.5° E) have an opposite sign compared with those of the mean trend. As with the results for the mean trend, AMM_SON(−1), AMO_SON(−1), and NAO_JJA(−1) present the strongest correlations ($\overline{\left|r_s\right|}$ = 0.56, 0.60, and 0.36, respectively) with the variance trend of AM daily rainfall over South Korea. Therefore, these three SCIs were finally selected as the significant climate indices linked closely to the long-term trends of AM daily rainfall over South Korea.

4.3. Appropriate Model Selection Considering Model Fit and Uncertainty

The combinations of three significant SCIs selected in Section 4.2 (AMM_SON(−1), AMO_SON(−1), NAO_JJA(−1)) were used as covariates of location and scale parameters in the GEV distribution. As presented in

Table 5. Description of 59 GEV candidate models in this study.

Distribution Type	GEV Parameters			The Number of Models
	Location, $\mathit{\mu }\left(\mathit{t}\right)$	Scale, $\mathit{\sigma }\left(\mathit{t}\right)$	Shape, $\mathit{\xi }$
ST-GEV	$\mu$	$\sigma$	$\xi$	1
NS-GEV(1,0,0) with time covariate	${\mu }_0+{\mu }_{1}t$	$\sigma$	$\xi$	1
NS-GEV(1,1,0) with time covariate	${\mu }_0+{\mu }_{1}t$	$\exp \left({\sigma }_0+{\sigma }_{1}t\right)$	$\xi$	1
NS-GEV(1,0,0) with SCI covariates	(1) ${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)$ (2) ${\mu }_0+{\mu }_{1}SC{I}_2\left(t\right)$ (3) ${\mu }_0+{\mu }_{1}SC{I}_3\left(t\right)$ (4) ${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)+{\mu }_{2}SC{I}_2\left(t\right)$ (5) ${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)+{\mu }_{2}SC{I}_3\left(t\right)$ (6) ${\mu }_0+{\mu }_{1}SC{I}_2\left(t\right)+{\mu }_{2}SC{I}_3\left(t\right)$ (7) ${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)+{\mu }_{2}SC{I}_2\left(t\right)+$ (8) ${\mu }_{3}SC{I}_3\left(t\right)$	$\sigma$	$\xi$	7
NS-GEV(1,1,0) with SCI covariates	(1) ${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)$ (2) ${\mu }_0+{\mu }_{1}SC{I}_2\left(t\right)$ (3) ${\mu }_0+{\mu }_{1}SC{I}_3\left(t\right)$ (4) ${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)+{\mu }_{2}SC{I}_2\left(t\right)$ (5) ${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)+{\mu }_{2}SC{I}_3\left(t\right)$ (6) ${\mu }_0+{\mu }_{1}SC{I}_2\left(t\right)+{\mu }_{2}SC{I}_3\left(t\right)$ (7) ${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)+{\mu }_{2}SC{I}_2\left(t\right)+$ (8) ${\mu }_{3}SC{I}_3\left(t\right)$	(1) $\exp \left({\sigma }_0+{\sigma }_{1}SC{I}_1\left(t\right)\right)$ (2) $\exp \left({\sigma }_0+{\sigma }_{1}SC{I}_2\left(t\right)\right)$ (3) $\exp \left({\sigma }_0+{\sigma }_{1}SC{I}_3\left(t\right)\right)$ (4) $\exp \left({\sigma }_0+{\sigma }_{1}SC{I}_1\left(t\right)+{\sigma }_{2}SC{I}_2\left(t\right)\right)$ (5) $\exp \left({\sigma }_0+{\sigma }_{1}SC{I}_1\left(t\right)+{\sigma }_{2}SC{I}_3\left(t\right)\right)$ (6) $\exp \left({\sigma }_0+{\sigma }_{1}SC{I}_2\left(t\right)+{\sigma }_{2}SC{I}_3\left(t\right)\right)$ (7) $\exp ({\sigma }_0+{\sigma }_{1}SC{I}_1\left(t\right)+{\sigma }_{2}SC{I}_2\left(t\right)+$ (8) ${\sigma }_{3}SC{I}_3\left(t\right))$	$\xi$	49

Note. $SC{I}_1$: AMM_SON(−1), $SC{I}_2$: AMO_SON(−1), $SC{I}_3$: NAO_JJA(−1).

, 59 GEV distributions including ST-GEV and NS-GEV with time covariates were considered as candidate models. Following the model selection procedure described in Section 3.5, an appropriate GEV distribution was selected for each station.

Table 6. Model selection result at the Daegu station (No. 143).

Distribution Type	rAIC Value	GEV Parameters			Uncertainty,
		Location, $\mathit{\mu }\left(\mathit{t}\right)$	Scale, $\mathit{\sigma }\left(\mathit{t}\right)$	Shape, $\mathit{\xi }$	$\mathit{C}{\mathit{I}}_{{\mathit{Q}}_{100}}^{95\%}$
NS-GEV(1,0,0) with SCI covariates	0	${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)$	$\sigma$	$\xi$	102.2
NS-GEV(1,0,0) with SCI covariates	0.24	${\mu }_0+{\mu }_1\mathit{S}\mathit{C}{\mathit{I}}_2\left(\mathit{t}\right)$	$\sigma$	$\xi$	88.8
NS-GEV(1,1,0) with SCI covariates	1.69	${\mu }_0+{\mu }_{1}SC{I}_2\left(t\right)$	$\exp \left({\sigma }_0+{\sigma }_{1}SC{I}_3\left(t\right)\right)$	$\xi$	116.6
ST-GEV	1.74	$\mu$	$\sigma$	$\xi$	109.5
NS-GEV(1,1,0) with SCI covariates	1.80	${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)$	$\exp \left({\sigma }_0+{\sigma }_{1}SC{I}_2\left(t\right)\right)$	$\xi$	121.3
NS-GEV(1,0,0) with SCI covariates	1.81	${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)+{\mu }_{2}SC{I}_2\left(t\right)$	$\sigma$	$\xi$	103.1
NS-GEV(1,1,0) with SCI covariates	1.81	${\mu }_0+{\mu }_{1}SC{I}_2\left(t\right)$	$\exp \left({\sigma }_0+{\sigma }_{1}SC{I}_2\left(t\right)\right)$	$\xi$	109.7
NS-GEV(1,0,0) with SCI covariates	1.82	${\mu }_0+{\mu }_{1}SC{I}_2\left(t\right)+{\mu }_{2}SC{I}_3\left(t\right)$	$\sigma$	$\xi$	108.2
NS-GEV(1,0,0) with SCI covariates	1.95	${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)+{\mu }_{2}SC{I}_3\left(t\right)$	$\sigma$		115.6
NS-GEV(1,1,0) with SCI covariates	1.95	${\mu }_0+{\mu }_{1}SC{I}_1\left(t\right)$	$\exp \left({\sigma }_0+{\sigma }_{1}SC{I}_3\left(t\right)\right)$		124.7

Note. $SC{I}_1$: AMM_SON(−1), $SC{I}_2$: AMO_SON(−1), $SC{I}_3$: NAO_JJA(−1).

presents the model selection result at the Daegu station (No. 143) as an example. Among the 59 GEV candidate models, 9 NS-GEVs and ST-GEV with $rAIC\le 2$ were selected to assess the uncertainty. Thereafter, the bootstrap time-averaged 95% CIs for the 100-year quantile ($C{I}_{Q_{100}}^{95\%}$) were calculated. As a result, the NS-GEV(1,0,0) that uses AMO_SON(−1) as covariates of location parameter and constant scale parameter was finally selected as an appropriate GEV model (as shown in bold in Table 5). To evaluate the quality of the model fit, the QQ-plots were also examined. Figure 8 presents the QQ-plots of the ST-GEV, NS-GEV with time covariate, and selected NS-GEV with SCI covariates at the Daegu station (No. 143). The selected NS-GEV with SCI covariate provides a better fit than the ST-GEV and NS-GEV with time covariate.

Figure 9 illustrates the model selection results at 61 stations over South Korea according to the distribution type classified in Table 5. The ST-GEV, NS-GEV with time covariate, and NS-GEV with SCI covariates were selected at 18 stations (29.5%), 3 stations (4.9%), and 40 stations (65.6%), respectively. It is worth noting that the NS-GEV with SCI covariates was selected at a majority of stations, whereas the NS-GEV with time covariate was selected at only 3 stations. In the case of the NS-GEV with SCI covariates, the number of stations where the NS-GEV(1,0,0) was selected was larger than that where the NS-GEV(1,1,0) was selected.

5. Discussions

5.1. Significant Climate Indices for Annual Maximum Rainfall over South Korea

Based on the correlation analysis with the long-term trend of AM daily rainfall over South Korea, the three significant SCIs (AMM_SON(−1), AMO_SON(−1), and NAO_JJA(−1)) were selected. In nonstationary extreme value modeling, it is necessary to verify that external variables such as climate indices are physically meaningful and justifiable to use as covariates [10]. Therefore, the physical mechanisms of teleconnections between the selected climate indices and the East Asian summer monsoons (EASMs), which mainly affect rainfall over South Korea, were investigated by previous literature as follows.

The AMO, which exhibits an oscillation of the sea surface temperatures (SSTs) in the Atlantic, is one of the most significant indices leading to global-scale climate variability, particularly in the Northern Hemisphere [57,58]. Although the influence of the AMO on the variability in EASMs is not fully understood, there have been several studies to explain the relationship between AMO and EASM variability [58,59,60,61]. Lu et al. [58] suggested that a positive AMO causes positive SST anomalies in the eastern Indian Ocean and maritime continent, and this atmospheric heating leads to an anticyclonic anomaly at the low troposphere over the western North Pacific, resulting in strong EASMs. Si and Ding [60] described that through the AMO-Northern Hemisphere teleconnection, the positive and negative AMO affect the SST anomalies in the midlatitude central Pacific and the geopotential height anomalies over the Huanghe–Huaihe River and Yangtze River, resulting in variability in summer precipitation over East Asia. Li et al. [61] explained that the tropical easterly jet (TEJ) is a significant component of the Afro-Asian monsoon system and directly affects rainfall near East Asia. They described that warm AMO phases cause intensification of the TEJ, which could increase rainfall over East Asia. On the other hand, there is a lack of studies on how the AMM physically affects the EASMs. However, the AMM, characterized by meridional variations in SST and sea level pressure (SLP) in the tropical Atlantic, is highly correlated with AMO and is affected by AMO on decadal time scales [62]. In addition, the AMM_SON(−1) and AMO_SON(−1) also present very similar patterns, as displayed in Figure 5 and Figure 6. These results imply that the AMM may affect the EASMs indirectly.

The NAO, which describes the variation of the sea level pressure between the Icelandic Low and the Azores High [63], is known to affect the weather and climate over East Asia in addition to the North Atlantic [64,65,66,67]. In particular, some studies have illustrated that the summer NAO affects the EASMs [68,69,70]. Linderholm et al. [68] suggested that the North Atlantic storm tracks and transient eddy activity related to the summer NAO can affect the EASMs. Sun and Wang [69] described that the summer NAO can cause a meridional dipole pattern over East Asia by changing the zonal wave activity over the Eurasian Continent, and it would finally contribute to summer rainfall over central and northern East Asia. Wang et al. [70] explained that the summer NAO affects the summer rainfall in East China through an intermediate bridge effect, resulting from the thermal forcing of the Tibetan Plateau.

Although there are some differences in the lag and time scale of teleconnections between the climate indices and the EASMs, the abovementioned studies support the acceptability of the selected climate indices by explaining the physical mechanism for teleconnections. In addition, this implies that the residue of AM daily rainfall extracted by the EEMD was properly used to detect the significant SCIs that impact extreme rainfall over South Korea. The selected SCIs could be employed as predictors in extreme value modeling for AM daily rainfall.

5.2. Consideration of Uncertainty in Model Selection

In nonstationary extreme value modeling, the uncertainty caused by the increase in the complexity of the model is a crucial issue [3,4]. However, information criteria such as AIC or BIC (Bayesian information criterion) are solely used to determine the best model [8,9,14,26,32,33]. Although these methods attempt to select a less complex model considering model parsimony, they are not enough to consider model uncertainty. In this study, a model selection procedure combining the information criterion and residual bootstrap method was proposed to consider both the model fit and uncertainty. That is, instead of selecting the best model with the smallest AIC value, distribution models whose $rAIC$ is equal to or less than 2 were first selected as acceptable candidates. Among the selected candidate models, the model with the smallest uncertainty measured by the residual bootstrap method was finally selected as the appropriate model. This stepwise procedure allows us to select the most appropriate model with acceptable goodness-of-fit and low uncertainty and to estimate more reliable design quantiles in practice. In addition, it reduces the analytical and computational burden of uncertainty assessment by narrowing down the number of model candidates. Even if several exogenous variables are used as covariates for nonstationary extreme value modeling, a number of candidate models can be considered due to probability distribution types [19], combinations of covariates [8,9], and their link function types [26,31,33]. Therefore, our suggested model selection procedure can be an efficient and reasonable way to select the most appropriate model.

5.3. Suggestions for Extreme Value Modeling with SCI Covariates

Extreme value modeling was conducted using the three significant SCIs (AMM_SON(−1), AMO_SON(−1), and NAO_JJA(−1)) as covariates of GEV parameters. Consequently, the ratio of stations where the NS-GEV with SCI covariates was selected was 65.6%, whereas the ratio of stations where the NS-GEV with time covariate was selected was 4.9%. This result indicates a notable improvement in extreme value modeling in terms of both model fit and uncertainty by applying the significant SCIs as covariates. Moreover, the NS-GEV with SCI covariates could rationally explain the long-term trend and variability of AM daily rainfall compared to ST-GEV and NS-GEV with the time covariate. In nonstationary extreme value modeling, much attention has focused on modeling the relationship between extreme data and their attributions [71,72], instead of modeling the linear or polynomial trend in extreme data using time t. Our results also support the use of climate indices as covariates for modeling extreme rainfall. Therefore, here, we propose to identify the significant climate indices properly based on EEMD and perform extreme value modeling using the selected climate indices to estimate more reliable rainfall quantiles in practice.

5.4. Limitations of this Study

One limitation of this study is that we only consider the GEV distribution to conduct nonstationary extreme value modeling, as the GEV distribution is recommended for the estimation of the design quantile of extreme rainfall in South Korea [73]. Other distribution models, such as generalized logistic, generalized normal, and Pearson type III distributions, can also be the candidate distributions for fitting AM rainfall series [74,75,76]. To the best of our knowledge, there are a limited number of studies that have considered nonstationarity in other distributions so far [35,77]. Therefore, we believe that this study can be improved by applying nonstationarity to various probability distributions and comparing their performance with that of NS-GEVs for future studies.

Another limitation is the sensitivity of data length to the result. In this study, we used observed rainfall and climate indices available until 2018. Due to the chaotic nature of the climate system, there is a possibility of providing varying results by adding the latest observations. For example, the SCI may be differently selected because there is no significant difference in $P_{sig}$ between seasons (see Table 3 and Table 4). This can lead to another source of uncertainty in nonstationary extreme value modeling. Thus, we remain this limitation for future studies that should be assessed for the selection of significant SCIs.

Last but not least, we employ the nonstationary extreme value modeling based on the “user-friendly” method, not considering the ergodicity issue, as we mentioned earlier in the Introduction section. Thus, this study estimated parameters of model candidates with all possible combinations of SCIs by directly applying the MLE method to the nonstationary data. Then, we selected the best model by comparing all candidate models. To deal with the ergodicity issue, it is necessary to define the deterministic relationships between the selected SCI and extreme rainfall based on the physical mechanism prior to parameter estimation. Since the significant SCIs that affect the long-term trend of extreme rainfall were identified throughout this study, we believe further study can be extended to adopt the ergodicity issue together with the current study.

6. Conclusions

In nonstationary extreme value modeling, it is essential to reflect the trend in statistical characteristics, such as the mean and variance of the observations, because nonstationarity is generally considered by the time-dependent location and/or scale parameters of the probability distribution model. Using physically meaningful information, such as climate indices as covariates, a probability distribution model can be improved in terms of not only considering nonstationarity but also reducing uncertainty. In this study, the procedure of extreme value modeling with large-scale climate modes was proposed, from the identification of appropriate SCIs impacting regional AM rainfall based on EEMD to the selection of appropriate NS-GEV, considering both model fit and uncertainty. Using the EEMD, residues that indicate the long-term trends in the mean and variance were extracted from AM daily rainfall. Subsequently, the correlation coefficient was calculated between the extracted residues and various lagged SCIs. After the identification of appropriate SCIs, the AMM_SON(−1), AMO_SON(−1), and NAO_JJA(−1) were finally selected as significant SCIs with an impact on the long-term trends in both the mean and variance of AM daily rainfall over South Korea. The combinations of these significant SCIs were employed as covariates of location and scale parameters for constructing various NS-GEVs.

As the uncertainty increases with the complexity of the probability distribution model, it is necessary to consider the uncertainty to select an appropriate probability distribution model. However, there have been few studies that consider uncertainty in the model selection procedure. In this study, the model selection procedure combining the AIC and residual bootstrap method was proposed to select an appropriate GEV model among ST-GEV, NS-GEV with time covariate, and NS-GEV with SCI covariates. Thus, the NS-GEV with SCI covariates was selected as an appropriate model for more than 65% of the applied stations in South Korea.

Nowadays, uncertainty is the main issue in nonstationary extreme value modeling. Rather than simply using the time covariate, employing appropriate climate indices as covariates can reduce uncertainty. The EEMD could be employed to detect the significant climate indices for using the covariates of the nonstationary extreme value model because it can extract the long-term trend inherent in the time series. Further, by selecting an appropriate probability distribution model considering both model fit and uncertainty, more accurate and reliable quantiles could be estimated. Although this study focused on the determination of significant climate indices and their application to extreme value modeling of extreme rainfall over South Korea, the method and discussion presented are expected to be extended to various hydrological variables, as well as to other regions.