Future Projections and Uncertainty Assessment of Precipitation Extremes in the Korean Peninsula from the CMIP6 Ensemble with a Statistical Framework

Abstract

Scientists occasionally predict projected changes in extreme climate using multi-model ensemble methods that combine predictions from individual simulation models. To predict future changes in precipitation extremes in the Korean peninsula, we examined the observed data and 21 models of the Coupled Model Inter-Comparison Project Phase 6 (CMIP6) over East Asia. We applied generalized extreme value distribution (GEVD) to a series of annual maximum daily precipitation (AMP1) data. Multivariate bias-corrected simulation data under three shared socioeconomic pathway (SSP) scenarios—namely, SSP2-4.5, SSP3-7.0, and SSP5-8.5—were used. We employed a model weighting method that accounts for both performance and independence (PI-weighting). In calculating the PI-weights, two shape parameters should be determined, but usually, a perfect model test method requires a considerable amount of computing time. To address this problem, we suggest simple ways for selecting two shape parameters based on the chi-square statistic and entropy. Variance decomposition was applied to quantify the uncertainty of projecting the future AMP1. Return levels spanning over 20 and 50 years, as well as the return periods relative to the reference years (1973–2010), were estimated for three overlapping periods in the future, namely, period 1 (2021–2050), period 2 (2046–2075), and period 3 (2071–2100). From these analyses, we estimated that the relative increases in the observations for the spatial median 20-year return level will be approximately 18.4% in the SSP2-4.5, 25.9% in the SSP3-7.0, and 41.7% in the SSP5-8.5 scenarios, respectively, by the end of the 21st century. We predict that severe rainfall will be more prominent in the southern and central parts of the Korean peninsula.

1. Introduction

Heavy precipitation can have cascading effects on communities, infrastructure, agriculture, and livestock, as well as on economically and culturally important natural ecosystems. For example, extreme rain events can result in costly damage to wastewater treatment plants, culverts, and roads. Extreme precipitation can result in landslides and floods, accompanied with a loss of life and the deterioration of infrastructure. Thus, understanding and projecting heavy rainfall is of significant importance to climate change impact, adaptation, and vulnerability assessments.

Numerous studies have reported that extreme precipitation events have become more frequent during the last century, and are occurring even more often over the 21st century ([1,2,3,4,5,6,7], for example). A simplified and major reason for more frequent extreme rainfall is the following: Warming conditions mean more evaporation, which leads to more water vapor in the air. When rain-triggering conditions are favorable, more saturated air leads to heavier precipitation [8,9]. This has been the case across some areas of the world during the last century [10]. For example, heavy daily rainfall has accounted for an increased proportion of total annual rainfall over an increasing fraction of the Australian continent since the 1970s ([11], Section 4.2.3). The average annual rainfall has already increased by nearly 50% over parts (including the Netherlands, Belgium, and Luxembourg) of northern Europe [12]. This trend is likely to be accelerated with increased global warming over the 21st century [10,12]. Some studies have projected that global warming leads to a higher intensity of precipitation and longer dry periods, for example, in Europe and Asia [13,14,15,16]. Seemingly paradoxically, as written by Mann and Kump [12], “While many regions are likely to become drier, scientists predict that even in those regions individual rainfall or snowfall events will become more intense, although longer dry spells will separate them.”

Extreme rainfall occurs frequently over the Korean peninsula during the warm season from June through September, in association with synoptic disturbances, typhoons, or convective changes within the air masses over the region. Previous studies [17,18,19] have reported an increase in observed extreme precipitation in Korea. Lee et al. [20] predicted that the increasing changes in the future heavy rainfall across East Asia appear more distinctly in Korea at a local scale, which indicates a higher sensitivity of the Korean peninsula to global warming. It is thus crucial to project and assess the changes in extreme precipitation events in Korea under different scenarios.

Some authors [21,22,23,24,25,26,27] have predicted increasing changes in future extreme rainfall over Korea by using a single model, ensembles of regional climate models, or multiple Coupled Model Inter-Comparison Project Phase 5 (CMIP5) models. In this study, we update the previous studies based on the multiple CMIP6 models under the three shared socioeconomic pathway (SSP) scenarios—namely, SSP2-4.5, SSP3-7.0, and SSP5-8.5 [28].

Studies on the projection of future climate change have used ensembles of multiple climate simulations. Multi-model ensemble (MME) methods of climatic projection have been proven to improve upon the systematic bias and to have fewer of the general limitations that are typically associated with single simulation models. Among the many ensemble methods, model weighting or averaging is typically employed ([29,30,31], for example). Model averaging is a statistical method in which unequal or equal weights are assigned to those models. Despite some arguments, the equal weighting or “model democracy” [30] has been criticized because it does not take into account the performance, uncertainty, and independency of each model in constructing an MME ([32,33,34], for example).

One typical unequal weighting scheme involves giving more weights to those models that are more skillful and realistic for a specific process or application. This performance-based weighting method and its variants, including Bayesian model averaging (BMA), have been employed in many different studies. It has improved the accuracy of the projections and reduced the prediction uncertainty. However, it has been reported that only a few models often exhibit extremely high weights, and most others have very low weights [27,35]. This phenomenon may be because some models are more fit to the observations for given applications than others, and thus, they receive extremely high weights in a multi-model estimate of change [36]. Such an aggressive weighting based on performance will only be dangerous in the sense of an overfitting when observational uncertainty is large, and thus, it is not robust in quantifying the uncertainty.

In addition to the performance, some researchers have considered other criteria, such as model convergence [37], model independency [32,38,39,40], and a semi-performance measure [41]. A weighting scheme that accounts for both the independence and performance simultaneously is called the PI-weighting. In this study, we employ PI-weighting to robustly quantify uncertainty in an MME. In calculating the PI-weights, considering only one or two climate variables over a relatively small area can lead to an overfitting problem [36,40]. To avoid this problem, we thus consider five climate variables over East Asia, while our focus is the annual maximum daily precipitation (AMP1) over the Korean peninsula.

In applying the PI-weighting, we have to determine two shape parameters that control the strength of the weights. One way to select the shape parameters is a leave-one-out perfect model test [36,40], but it requires a huge amount of computing time. To overcome this computational problem, we suggest simple ways to determine these parameters based on the entropy and p-values of the chi-square statistic.

The remainder of this paper is structured as follows. The data construction and numerical models are described in Section 2. The statistical methods are briefly mentioned in Section 3. Section 4 describes the PI-weighting with the computational details, including simple ways for determining the shape parameters. The results of the model weights and projected future changes are presented in Section 5 and Section 6, respectively. Section 7 describes the results of an uncertainty assessment and projection by latitude based on an analysis of variance. In Section 8, relative improvement of the PI-weighted method over the simple average is quantified using the skill score and prediction variance. Discussions are then given in Section 9, followed by a summary of the paper in Section 10. Details including technical specifics, tables, and figures are provided in the accompanying Supplementary Material (hereafter referred to as the Supplementary Materials) file.

2. Data and Simulation Models

We consider five climate variables over East Asia to avoid overfitting in calculating the PI-weights, while our focus is the annual maximum daily rainfall (AMP1) over the Korean peninsula in this study.

Table 1. The five climate variables considered in this study.

Variable Acronym	Description
AMP1	Annual Maximum Daily Precipitation
AMP5	Annual Maximum Five-Day Precipitation
ATP	Annual Total Precipitation
AMCWD	Annual Maximum Consecutive Wet Days
AMCDD	Annual Maximum Consecutive Dry Days

lists five climate variables.

Consecutive wet days and dry days are defined as consecutive days with daily precipitation of ≥1 and <1 mm, respectively [42,43]. Table S1 in the accompanying Supplementary Materials lists the 21 CMIP6 climate models used in this study. The considered scenarios are shared socioeconomic pathways SSP2-4.5, SSP3-7.0, and SSP5-8.5 [28]. Hereafter, we shorten the scenario names to SSP2, SSP3, and SSP5, although we use both versions interchangeably. Three overlapping periods are considered for future data, namely, period 1 (2021–2050), period 2 (2046–2075), and period 3 (2071–2100), abbreviated by P1, P2, and P3 in this study.

To re-grid the common grid points of $1^{\circ }\times 1^{\circ }$, the iterative Barnes interpolation scheme [44] was employed for the observations and simulation data from the 21 models for each of five climate variables. The Barnes technique produces a rainfall field on a regular grid from irregularly distributed rainfall observation stations.

East Asia (EA) is considered by many authors (e.g., [4,14,15,45,46]) to be one of the regions most vulnerable to future increases in extreme weather and climate under global warming. For example, projected changes in extreme precipitation indices over the Asian monsoon domain are larger than over other monsoon domains, indicating the strong sensitivity of Asian monsoons to global warming [47].

Figure 1 depicts maps of EA and the Korean peninsula (KP) showing the observations sites and 1180 points over $1^{\circ }\times 1^{\circ }$ grids. The observations for the 38-year reference period (1973–2010) were obtained from the Meteorological Administrations of China, Japan, and Korea. In China, observations from 726 stations were used [48] to construct data for 1056 grids. In Japan, observations from 1325 sites were used [49] to construct data for 78 grids.

There are 64 and 27 observation stations in South and North Korea, respectively [50]. Because the 27 stations (black circles in Figure 1) in North Korea are too sparse, re-gridding using Barnes interpolation may not capture the locality of severe rainfall there. For a better re-gridding in North Korea, we used the Asian Precipitation Highly Resolved Observational Data Integration Towards Evaluation (APHRODITE) reanalysis data [51] as auxiliary information. However, the AMP1 in APHRODITE has serious bias in its mean and variance, as shown in Figure S1. We thus correct it using nearby observations using the quantile mapping technique [52]. Examples of time series plots of the observations, APHRODITE data, and the bias-corrected APHRODITE data near the observational stations are shown in Figure S1. Then, the bias-corrected APHRODITE data and the observations at stations were used together for better interpolation in re-gridding. On the right side of Figure 1, the map of the KP shows the spatial distribution of 91 rainfall observation stations, the 80 APHRODITE grid points in North Korea, and the final 46 grids used in this study.

3. Methods

3.1. Generalized Extreme Value Distribution

The generalized extreme value distribution (GEVD) is widely used to analyze extreme univariate values. The three types of extreme value distributions are sub-classes of GEVD. The cumulative distribution function of the GEVD is as follows:

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

(1)

when $1+\xi (x-\mu )/\sigma >0$, where $\mu ,\sigma$, and $\xi$ are the location, scale, and shape parameters, respectively. The particular case for $\xi =0$ in Equation (1) is the Gumbel distribution, whereas the cases for $\xi >0$ and $\xi <0$ are known as the Fréchet and the negative Weibull distributions, respectively [53].

It can be helpful to describe the changes in extremes in terms of the changes in extreme quantiles. These are obtained by inverting the following (1): $z_p=\mu -\frac{\sigma }{\xi }[1-{\{-log(1-p)\}}^{-\xi }]$, where $G\left(z_p\right)=1-p$. Here, $z_p$ is known as the return level associated with the return period $1/p$ because level $z_p$ is expected to be exceeded, on average, once every $1/p$ years [53]. For example, a 20-year return level is computed as the 95th percentile of the fitted GEVD and a 50-year return level as the 98th percentile. Conversely to the above, the return period $T\left(z\right)=1/p$ for the given value z is obtained by calculating $p=1-G\left(z\right)$. For the given value z, $T\left(z\right)$ is sometimes called the expected waiting time, and the value $p=1-G\left(z\right)$ is referred to as the exceedance probability of z.

Assuming the data approximately follow a GEVD, the parameters can be estimated by the maximum likelihood method [53,54] or the method of L-moment estimation. The L-moment estimator is more efficient than the maximum likelihood estimator in small samples for typical shape parameter values [55]. The L-moment method is employed in this study using the “lmom” package in R [56] because a relatively small number of samples are analyzed for each comparison period. Moreover, the formulae used to obtain the L-moment estimator are simple compared to those for obtaining the maximum likelihood estimator, which needs an iterative optimization until convergence.

3.2. Bias Correction

Although simulations from climate or meteorological models provide significant information, the simulated data are associated with potential biases in that their statistical distribution differs from the distribution of the observations. This is partly because of unpredictable internal variability that differs from the observations, and because global climate models (GCMs) have a very low spatial resolution for being employed directly in most impact models [57,58]. For example, in GCM precipitation fields, the bias may be due to errors in convective parameterizations and unresolved subgrid-scale orography [59]. Bias-correction (BC) methods are commonly applied to transform the simulated data into new data with no or fewer statistical biases with respect to an observed time series. In this study, future simulation outputs are bias corrected to compute the intensity of extreme rainfall.

To correct the model outputs more efficiently by taking account of the dependency among variables or nearby grids, several multivariate BC methods have recently been proposed. In this study, we chose the multivariate bias correction (MBC) method by Cannon [59] among the many available BC methods [52]; it is a multivariate extension of quantile delta mapping (QDM). QDM has an advantage of preserving the approximate trends of the model data. In this study, the MBC is applied to the five climate variables provided in Table 1 to take into account the dependency among these variables. The MBC is applied one time to the future simulation data for all periods (2021–2100) based on the historical information, not for each period separately. More details are provided in the Supplementary Materials.

4. Weighting Method for Ensembles

4.1. Performance and Independence Weighting

Knutti et al. [39] argued that the growing number of models with different characteristics and considerable interdependence finally justifies abandoning a strict model democracy. They provided at least five reasons for why PI-weighting is required. Brunner et al. [36] illustrated that PI-weighting leads to an increase in the investigated skill score for temperature and precipitation while minimizing the probability of overfitting (or overconfidence).

As the basic idea of PI-weighting, models that agree poorly with observations for a selected set of diagnostics receive less weight, as do models that largely duplicate existing models [39]. Weights are calculated for each model based on a combination of the distance $D_i$ (informing the performance) and the model similarity $S_{ij}$ (informing the dependence):

$$w_i=\frac{exp(-\frac{D_i}{{\sigma }_D})}{1+{\sum }_{j\ne i}^{M}exp(-\frac{S_{ij}}{{\sigma }_S})},$$

(2)

with the total number of model runs M and the shape parameters ${\sigma }_D$ and ${\sigma }_S$. The weights are normalized such that their sum equals 1. The details in computing the performance distance $D_i$ are given in the Supplementary Materials.

The numerator represents the modeling skill when using a Gaussian weighting, where the weight decreases exponentially the farther away a model is from the observations. The denominator is the “effective repetition of a model” [32] and is intended to account for the model interdependency [39]. The details in computing the distance $S_{ij}$ are given in the next subsection.

The shape parameters define the strength of the weighting and the relative importance of the performance and independence [36]. Large values will lead to an almost equal weighting, whereas small values will lead to aggressive (or one-sided) weighting, giving a few models most of the weight. The shape parameters are often determined through a perfect model test (or a model-as-truth experiment) using the continuous rank probability score [36,40]. The perfect model test picks each model from a multi-model ensemble in turn and treats it as the true representation of the climate system. This leave-one-out procedure requires a huge amount of computing time. To address this computational problem, we consider relatively simple ways to determine the shape parameters in the next subsections.

4.2. Computing Independence Weights

If a model has no close neighbors, then all $S_{ij}(i\ne j)$ are large, and the denominator of the PI-weight is approximately one and has no effect. If two models i and j are identical, then $S_{ij}=0$ and the denominator equals two, so each model gets half the weight.

To calculate the model similarity $S_{ij}$, we follow a technique among several methods proposed by Sanderson et al. [60]. A method employed in this study is based on the empirical orthogonal function (EOF) or principal component analysis. The following process is done for each grid: First, for each model, the historical data from 1850 to 2015 and the future simulation data from three scenarios are lined up as one time series data, as in Figure S2 in the Supplementary Materials. The bias correction is not applied to all data for this process. We can choose the historical data only, as was done by Brunner et al. [61], but we deploy all simulation data for a maximum use of available information. For each time series induced from each model, seven-year moving averages are obtained. Then, for each climate variable, a correlation matrix R among all M models is constructed by applying the Spearman correlation coefficients to those M numbers of the series of seven-year moving averages. That is, R is the correlation matrix of M models, with size $M\times M$.

A singular value decomposition (SVD) is performed on $R^{1/2}$ and truncated to t modes to obtain the dominant modes of multivariate ensemble variability such that

$$R^{1/2}=U\lambda V^T,$$

(3)

where U is an orthogonal matrix of model loadings (size M by t) whose columns are the eigenvectors of the model correlation matrix R, $\lambda$ (size t by t) are the eigenvalues of R, and V (size M by t) are the eigenvectors of R. The dimensions are sorted by decreasing eigenvalue, such that the basis set can be truncated to a smaller number of modes t [60]. Note that t is often determined by selecting a number of the eigenvalues greater than 1.

The model loadings U now define a t-dimensional space (where t is the truncation length of the SVD) in which intermodel and observation–model Euclidean distances may be defined. The intermodel distances can then be measured in a Euclidean sense in the loadings matrix, such that the distances $S_{ij}$ between two models i and j can be expressed as [60]

$$S_{ij}={\left\{\sum _{l=1}^t{[U(i,l)-U(j,l)]}^2\right\}}^{1/2}.$$

(4)

$U(i,l)$ is interpreted as a correlation or a dependency of the model i to the l-th principal component. Thus, a small $S_{ij}$ value means high dependency or similarity between models i and j.

The above procedures are done separately for each of five climate variables listed in Table 1. Thus, we have five distances between models i and j, corresponding to five climate variables. Then, the final distances are obtained by averaging those five distances for models i and j. An example of the final distances $S_{ij}$ between two models i and j averaged from the five distances is given in Table S2 in the Supplementary Materials. A small value indicates high dependency or similarity between two models. In Table S2, the first four models (UKESM2, CanESM5, EC-Earth3-Veg, and EC-Earth3) show the highest similarity, whereas the last four models (MPI-ESM1-2-LR, GFDL-ESM4, INM-CM5-0, and MPI-ESM1-2-HR) show the lowest similarity.

4.3. Selection of ${\sigma }_S$

To select an appropriate value of the shape parameter ${\sigma }_S$ for the I-weights, we consider an entropy-based approach. We denote $I_i\left({\sigma }_S\right)$ as a normalized I-weight for model i and for the given ${\sigma }_S$, as defined in the following:

$$I_i\left({\sigma }_S\right)=\frac{s_i\left({\sigma }_S\right)}{{\sum }_{l=1}^M{s}_i\left({\sigma }_S\right)},$$

(5)

where $s_i\left({\sigma }_S\right)=\frac{1}{1+{\sum }_{j\ne i}^{M}exp(-\frac{S_{ij}}{{\sigma }_S})}$.

The entropy of the I-weights as a measure of uncertainty [62] from these weights is defined by the following:

$$E\left({\sigma }_S\right)=-\sum _{i=1}^M{I}_i\left({\sigma }_S\right)log{I}_i\left({\sigma }_S\right)$$

(6)

as a function of ${\sigma }_S$. When all $I_i\left({\sigma }_S\right)$s are almost equal, the entropy has a high value. We thus expect the entropy to increase because ${\sigma }_S$ has a large value. Note that the calculation of $S_{ij}$ does not depend on ${\sigma }_S$, and thus, the $S_{ij}$ values obtained are fixed for the entropy computation. The entropy is computed as ${\sigma }_S$ changes from 0.1 to 1.0 in increments of 0.01.

Figure 2 presents the entropy function of ${\sigma }_S$ computed from the data used for this study, which indicates that it is at its minimum at ${\sigma }_S=0.4$. It is interesting to note that the entropy function increases as ${\sigma }_S$ decreases from 0.4 to zero. This is explained by looking into the similarity measure $1+{\sum }_{j\ne i}^{M}exp(-\frac{S_{ij}}{{\sigma }_S})$. As ${\sigma }_S$ moves toward zero, this measure converges at one for all i. Thus, $s_i$ moves toward one, and $I_i$ is close to $1/M$ for all i. Because we want to have a shape parameter ${\sigma }_S$ that can differentiate the I-weights most distinctly with minimum uncertainty, the value ${\sigma }_S=0.4$ minimizing the entropy is chosen in this study.

4.4. Selection of ${\sigma }_D$

To select an appropriate value of ${\sigma }_D$ for the P-weights, we attempted to use the entropy criteria again, but we were not fortunate enough to obtain the optimal result, as with ${\sigma }_S$. Thus, a technique based on the p-value of the chi-square statistic [54] is considered in this study.

We denote $P_i\left({\sigma }_D\right)$ as a normalized P-weight for model i and for the given ${\sigma }_D$, which is defined as follows:

$$P_i\left({\sigma }_D\right)=\frac{exp(-\frac{D_i}{{\sigma }_D})}{{\sum }_{l=1}^{M}exp(-\frac{D_l}{{\sigma }_D})},$$

(7)

where $D_i$ is the performance measure of the i-th model.

For testing the hypothesis frame, the null hypothesis and the alternative hypothesis for $i=1,\cdots ,M$ are as follows:

$$\begin{array}{c}{H}_0:\mathrm{all}\mathrm{weights}\mathrm{are}\mathrm{equal}\iff P_i=\frac{1}{M}\mathrm{for}\mathrm{all}i\\ H_1:\mathrm{some}\mathrm{weights}\mathrm{are}\mathrm{not}\mathrm{equal}\iff P_i\ne \frac{1}{M}\mathrm{for}\mathrm{some}i.\end{array}$$

(8)

For the given $P_i$, the chi-square statistic used to test the above hypothesis is as follows:

$${\chi }_0^2\left({\sigma }_D\right)=\sum _{i=1}^M\frac{{(\frac{1}{M}-P_i\left({\sigma }_D\right))}^2}{\frac{1}{M}}.$$

(9)

Because we do not want to accept equal weights, ${\sigma }_D$ should be selected to reject the null hypothesis. That is, the p-value [54] obtained from the chi-square statistic should be less than a preassigned value $\alpha$ (significance level), e.g., $\alpha =0.05$. In addition, because we also do not want aggressive weights, a ${\sigma }_D$ can be selected as the maximum value of ${\sigma }_D$ in which we still reject $H_0$ with $\alpha$ level. That is, our selection is

$${\sigma }_D^{\ast }=max\{{\sigma }_D:p\text{-}\mathrm{value}\left({\sigma }_D\right)<\alpha \},$$

(10)

where p-value$\left({\sigma }_D\right)=Pr[{\chi }^2>{\chi }_0^2\left({\sigma }_D\right)|H_0]$.

Here, ${\chi }^2$ indicates a random variable of (9) under the equal $P_i$ weights. Although this selection assures the use of the least aggressive weights, it is still statistically significantly different from the equal weights.

The p-values are computed by a Monte-Carlo simulation in which random numbers of weights are generated from the Dirichlet distribution [63]. When the parameters are all equal to 1, the Dirichlet distribution is the same as the multivariate uniform distribution with values between 0 and 1, which represents the null hypothesis. We used the “MCMCpack” package [64] in R to generate the random weights that satisfy $H_0$.

The detailed steps of computing the p-value for given ${\sigma }_D$ and ${\chi }_0^2\left({\sigma }_D\right)$ are:

• Step 1: Generate random weights $P_i^{\left(k\right)}$ from the Dirichlet distribution with all parameters equal to 1 (under $H_0$), for $i=1,\cdots ,M$;
• Step 2: Compute ${\chi }^2={\sum }_{i=1}^M\frac{{(\frac{1}{M}-P_i^{\left(k\right)})}^2}{1/M}$, and denote it ${\chi }_{\left(k\right)}^2$;
• Step 3: Iterate the above two steps K (=1000, for example) times;
• Step 4: Calculate $p\text{-}\mathrm{value}\left({\sigma }_D\right)={\sum }_{k=1}^K\frac{I[{\chi }_{\left(k\right)}^2>{\chi }_0^2\left({\sigma }_D\right)]}{K}$,

where $I\left[A\right]$ denotes the identity function, which takes 1 or 0, depending on if the condition A is satisfied or not. Note that $P_i^{\left(k\right)}$s generated in Step 1 do not depend on ${\sigma }_D$.

Figure 3 depicts the chi-square statistic values computed from AMP1 with some p-values as ${\sigma }_D$ changes from 0.10 to 0.30 in increments of 0.01. We calculated the ${\sigma }_D$ for each of the five climate variables, and then calculated the average from those five ${\sigma }_D$s. When $\alpha =0.05$, as is usually applied in testing a hypothesis in statistics, the averaged ${\sigma }_D^{\ast }$ from five different ${\sigma }_D$ is 0.21. When $\alpha =0.1$, ${\sigma }_D^{\ast }=0.25$. We use ${\sigma }_D^{\ast }=0.21$ in this study.

5. Results: Model Weights

5.1. Model Similarity

Table S2 provides the similarity values $S_{ij}$ for certain models. Figure 4 shows the intermodel distance matrix for the 21 CMIP6 models considered in this study. The distances are obtained from the five climate variables listed in Table 1 for East Asia. Each box represents a pairwise combination, where red indicates a greater distance. According to Figure 4, models MPI-ESM1-2-HR, INM-CM5-0, GFDL-ESM4, MPI-ESM1-2-LR, and INM-CM4-8 were found to be the most independent, whereas models UKESM, CanESM5, EC-Earth3-Veg, and EC-Earth3 were found to be the most dependent.

5.2. PI-Weights

The normalized PI-weights are obtained using Equation (2) with ${\sigma }_S=0.4$ and ${\sigma }_D=0.21$. Figure 5 demonstrates the distributions of the P-, I-, and PI-weights. The variability of the I-weights is smaller than that of the P-weights.

The high P-weights of the CanESM5 and EC-Earth3-Veg models decrease in the PI-weights owing to the low I-weights. The PI-weights of the BCC-CSM2-MR, FGOALS-g3, and GFDL-ESM4 models increase owing to a relatively high independency. The PI-weight is not located in the middle of the P- and I-weights, but is close to the P-weight, except in a few cases. When the P-weight (I-weight) is almost the same as the equal weight, as in the GFDL-ESM4 (KACE-1-0-G and IPSL-CM6A-LR) model, it seems that the PI-weight is wholly influenced by the I-weight (P-weight). The performances for some of the models, such as ACCESS-ESM1-5, MIROC6, INM-CM5-0, and INM-CM4-8, are so low that their (even relatively high) I-weights do not affect the final weights. Based on this view, the performance is more influential to the PI-weights than the independency. Some of these observations may be changed if different ${\sigma }_S$ and ${\sigma }_D$ are used.

6. Results: Future Projection of Extreme Precipitation

Using the PI-weights obtained in the above section, the future extreme precipitations are projected by the MME. Note that the future climate data are used after the bias correction with the MBC method [59]. Figure 6 illustrates the time series plots of the nine-year moving averages of AMP1 in Seoul, for example, from the observations, from the PI-weighted ensemble of the historical data, and from the PI-weighted ensembles for the future data from the three SSP scenarios, with a 90% confidence band. In Figure 6, the line for the observations shows more variation than the lines of the historical data and the future projections by the MME. This means that the MME has smaller variance and can catch the signal more clearly than a single simulation model.

6.1. Return Levels

Figure 7 (S3) displays boxplots of the 20-year (50-year) return levels of the AMP1 in the KP. The boxplot of the historical data after the BC is similar to that from the observations, whereas the boxplot before the BC is much smaller than that from the observations. The increasing trends from P1 to P3 are evident in every scenario. Summary statistics of the corresponding values of these boxplots are provided in Table S3.

These values are lower in the 20-year return levels than the results by Lee et al. [27], who used the BMA method with the CMIP5 models. The median values for P3 in Table S3 are higher than the results by Shin et al. [41], who used the CMIP6 models and a hybrid weighting method between the BMA and equal weighting. It seems that the hybrid weighting by Shin et al. [41] and the PI-weighting in this study are similar. These last two weighting schemes assigned relatively lower weights to the CMIP6 models with heavy future precipitation than a weighting method based on the performance only would give to the same models.

Figure 8 (S4) shows isopluvial maps of the 20-year (50-year) return levels of the AMP1 for the three future periods under the three scenarios. The difference by latitude is more evident than that by longitude. Jeju island, the southern coast, and the central region receive more downpours than the northern parts.

6.2. Changes in Return Levels

Figure S5 in the Supplementary Materials exhibits relative changes (unit: %) of the 20- and 50-year return levels of extreme precipitation over the Korean peninsula relative to the years 1973–2010. The increasing changes are more evident in the northern part than in the southern area. Summary statistics of relative changes are presented in Table S4 in the Supplementary Materials. The relative increase in the observations for the spatially averaged 20-year (50-year) return level was approximately 18.4% (18.9%) in the SSP2-4.5, 25.9% (26.3%) in the SSP3-7.0, and 41.7% (44.0%) in the SSP5-8.5 scenario by the end of the 21st century.

These rates of change are lower than the results by Lee et al. [27], who used the BMA method with the CMIP5 models. This is perhaps due to the differences in the reference period and the research methods, as well as the difference between the CMIP5 and CMIP6. For the 20-year return level, our result in the KP is approximately 1.8 to 2.0 times faster than the changes in the globally averaged value (10% in the SSP2-4.5 and 20% in the SSP5-8.5) reported by Kharin et al. [65].

6.3. Change in Return Periods

Figure 9 displays boxplots for the 20-year and 50-year return periods, as compared to the reference years (1973–2010) for the three future periods under the three scenarios. The corresponding statistics are presented in Table S5 in the Supplementary Materials. We realize that a 1-in-20 year (1-in-50 year) AMP1 in the Korean peninsula will likely become 1-in-11 (1-in-26) year, 1-in-8 (1-in-20) year, and 1-in-7 (1-in-15) year events in the median by the end of the 21st century based on the SSP2, SSP3, and SSP5 scenarios, respectively, as compared to the observations from 1973 to 2010. These findings indicate that both 20-year and 50-year return periods are likely to decrease by approximately 32% under the SSP2 and 60% under the SSP5 by the end of the 21st century.

These projections of extreme rainfall are similar to or less frequent than the results obtained by previous studies [23,27] based on CMIP5 models, i.e., approximately 1-in-11 (1-in-21) year and 1-in-7 (1-in-13) year events under SSP2 and SSP5, respectively. Shin et al. [41] realized the occurrence as approximately 1-in-10 (1-in-30) and 1-in-8 (1-in-17) year events under SSP2 and SSP5, respectively, which is similar to our projection.

6.4. Exceedance Probability and Waiting Time

Because of computational issues and the defects of the return period [66], the exceedance probability is often used as an alternative to the return period [67]. This is defined as $Pr[Y\left(\widehat{\theta }\right)>z]$, where z is a specified precipitation value and $Y\left(\widehat{\theta }\right)$ is a random variable following a GEVD with a parameter estimate $\widehat{\theta }$. Here, $Y\left(\widehat{\theta }\right)$ depends on the models, periods, and scenarios.

The spatially averaged estimates of the exceedance probability are presented in Figure S6 and Table S6. There are relatively large differences in the exceedance probability of a downpour of $100$ to 250 mm compared with that for over $250$mm, as shown in Figure S6. The differences between the past and future scenarios are distinct during the period P3.

From Table S6, the return period or expected waiting time ($T\left(z\right)$) until the reoccurrence of a specific AMP1 value (z) is computed by $T\left(z\right)=1/p\left(z\right)$, where $p\left(z\right)$ is the exceedance probability of the AMP1 z. These values are listed in

Table 2. The expected waiting time (unit: year) until reoccurrence or the return period of specific annual maximum daily precipitation (AMP1) values from 100 to 500 mm in the Korean peninsula obtained from the observations (OBS) and the CMIP6 model ensemble under the three SSP scenarios for the three future periods (P1, P2, and P3).

		SSP2-4.5			SSP3-7.0			SSP5-8.5
AMP1	OBS	P1	P2	P3	P1	P2	P3	P1	P2	P3
100 mm	1.8	1.4	1.3	1.3	1.4	1.3	1.2	1.4	1.3	1.2
150 mm	4.4	4.4	3.9	3.7	4.0	3.2	2.5	3.6	3.0	2.5
200 mm	11.1	14.5	11.1	10.6	13.2	10.2	7.1	13.0	8.0	5.9
250 mm	30.7	41.7	31.7	27.1	42.9	28.7	16.9	34.4	24.9	15.0
300 mm	78	93	76	55	104	59	42	76	50	29
400 mm	485	400	297	202	469	254	201	227	166	107
500 mm	1970	1167	814	577	1148	649	563	538	493	307

. For $z=200$mm of rainfall, for example, the expected waiting times until a reoccurrence are 11.1 years in the past, 13.0 years in the future period P1, 8.0 years in P2, and 5.9 years in P3 based on the SSP5 scenario. For the case of a $z=300$ mm downpour, the expected waiting times are 78 years in the past, 76 years in P1, 50 years in P2, and 29 years in P3 based on the SSP5 scenario.

6.5. Expected Number of Reoccurring Years

Another quantity that we can obtain is the expected number of reoccurrences during a certain period. By multiplying 30 years by the exceedance probability $p\left(z\right)$, we can estimate the expected frequency of such years over a period of 30 years in which we have more than z amount of AMP1 for a year. These values are given in Table S7. For $z=150$ mm (200 mm), as a specific example, during the last 30 years, we have experienced 6.8 (2.7) years in which AMP1 was greater than 150 mm (200 mm). In addition, we are likely to have an expected number of years of 8.2 (2.3) for the future period P1, 10.1 (3.8) for P2, and 12.1 (5.1) for P3 under the SSP5 scenario.

From this comparison, particularly for AMP1 of 100 to 300 mm, the expected number of years of occurrence for the future periods under the SSP5 scenario increases by approximately 1.25 times that over the last 30 years for P2 and 1.91 times that for P3 by the end of 21st century. These results are based on the spatially averaged values. When the exceedance probability is considered for each grid, different local results will be obtained.

7. Results: Projection by Latitude and Quantifying Uncertainty

Three or four major sources of climate projection uncertainty might be considered when trying to understand uncertainties in projected metrics: (1) climate model, (2) emission scenario, and (3) internal variation or randomness unexplained by other sources [68]. We apply the analysis of variance (ANOVA) technique [69] in this study. Other methods considered in [67] can also be applied.

7.1. Variance Decomposition with Interaction

Hawkins and Sutton [68] applied a method of variance decomposition to quantify uncertainty when assuming no interactions between the sources. In addition to their assumption, we added the interaction term between the model and scenario. Other factors, i.e., (4) future period and (5) location, such as the latitude applied in this study, are also considered in the ANOVA model. Although the period and latitude may not be real sources of uncertainty, they can be included as independent variables in the ANOVA model because they may affect the response variable. Here, the response variable is a 20-year return level, for example.

Figure 10 presents the proportions of variances contributed by each variable for each period. In P1, the variance by the model is the largest, whereas that by the scenario is the smallest. The proportion of internal variation (residual) decreases from P1 to P3. The variance proportion by model decreases from P1 to P3, whereas the variation due to the scenario increases significantly from P1 to P3.

It is notable that the variation contributed by the interaction between the model and scenario is relatively large during period 1, but becomes smaller during period 3. The details of this variation are presented in Figure 11, where the return levels from the models generally increase in mean and variance as the scenario changes from SSP2 to SSP5. However, the return levels of some models, such as CESM2-WACCM, ACCESS-ESM1-5, NorESM2-MM, MIROC6, MRI-ESM2-0, GFDL-ESM4, and INM-CM5-0, decrease from SSP2 to SSP3, contrary to the expectation. Moreover, the return values of some models, such as MPI-ESM1-2-HR, MPI-ESM1-2-LR, NorESM2-LM, and FGOALS-g3, decrease from SSP3 to SSP5. The latter two “unexpected” situations make the interaction variance large in P1, as shown in Figure 12. Despite such “unexpected” situations still appearing in P3, the variations owing to the scenario and model themselves are relatively too large, and thus, the proportion of the interaction variance becomes relatively small in P3.

7.2. Return Levels by Latitude

Figure 13 depicts boxplots of the 20-year return levels (unit: mm) from the 21 CMIP6 models as the latitude changes from south to north. It is evident that the return levels decrease initially from ${33}^{\circ }$ to ${36}^{\circ }$, but increase from ${36}^{\circ }$ to ${38}^{\circ }$. However, they decrease from ${38}^{\circ }$ to ${42}^{\circ }$ more rapidly. In Figure 13, the solid lines around the boxes across the latitude are the spatial medians of the 20-year return levels for each latitude and the four periods (P0, P1, P2, and P3). Here, P0 indicates the reference period, and thus, the result is obtained from the observations. The dashed lines depict the values for the three scenarios (SSP2-4.5, SSP3-7.0, and SSP5-8.5), in which three dashed lines locate the inside range of the solid lines. The spatial medians for the scenario (period) are obtained over all values across the periods (scenarios) for the 21 models. It is notable and interesting that the variation owing to the scenario is smaller than that from the period. Future projections fluctuate less than the observations.

The detailed values of the 20-year return levels of AMP1 and annual maximum five-day precipitation (AMP5) by latitude are presented in

Table 3. The 20-year return levels (unit: mm) of the annual maximum daily precipitation (AMP1) and of the annual maximum five-day precipitation (AMP5) averaged over each latitude for the observations (OBS) during 1973–2010 and for future periods—namely, P1 (2021–2050), P2 (2046–2075), and P3 (2071–2100). For each future period, the values were averaged over the three scenarios.

	20-Year Return Level (AMP1)				20-Year Return Level (AMP5)
Latitude	OBS	$\mathbf{P}\mathbf{1}$	$\mathbf{P}\mathbf{2}$	$\mathbf{P}\mathbf{3}$	OBS	$\mathbf{P}\mathbf{1}$	$\mathbf{P}\mathbf{2}$	$\mathbf{P}\mathbf{3}$
33	288	303	319	345	443	440	471	519
34	275	288	303	329	416	416	446	487
35	253	264	280	305	384	386	412	445
36	235	246	261	282	360	368	389	414
37	243	252	270	296	396	409	429	461
38	249	262	281	302	402	419	438	466
39	227	245	267	287	377	398	417	442
40	186	209	227	247	299	319	337	357
41	147	166	181	201	236	249	266	285
42	136	149	164	181	207	213	235	250
43	143	151	167	184	223	227	254	271

. These localized values provide different information from the previous results based on the spatial medians and averages over the Korean peninsula (Table S3). Table 3 indicates that a much higher intensity occurs in the southern part of the Korean peninsula than in the northern part. When the areas below and above ${38}^{\circ }$ are considered the south and north, respectively, the south is likely to receive approximately 1.54 (1.72) times heavier rainfall in the 20-year return levels of AMP1 on average than the north (the spatial median) based on a rough calculation using the above tables. Although we did not obtain these values, the localized values of the return periods, the relative changes in the return levels, the exceedance probabilities, and the expected number of occurring years can be obtained for each latitude. Those values in the south would show higher intensity during extreme precipitation with shortened expected waiting times compared to the spatial medians throughout the Korean peninsula.

8. Comparison of PI-Weighted Ensembles to the Simple Average

For comparison between the PI-weighting scheme and the simple average for the multimodel ensemble, we consider two measures. The first one is the error index $I^2$ based on Baker and Taylor [70] for the reference (historical) period. The second is the weighted variance of return level prediction [71,72] for the future periods. In computing these measures, non-bias-corrected data were used.

8.1. Error Index for the Historical Period

The error index can measure combined errors, that is, over multiple variables, compared to the observed climate [40]. However, we only evaluate our target variable AMP1 over the reference period 1973–2010. In a first step, the normalized error is calculated for the T-year return level as the difference between the PI-weighted multimodel and the observations:

$$e_w^2=\sum _{n=1}^N\frac{{(S_n-o_n)}^2}{{\sigma }_n^2},$$

(11)

where $S_n$ is the weighted multimodel T-year return level per grid point n. $S_n$ is calculated from the historical data that are not bias corrected. $o_n$ is the T-year return level obtained from the observations, and ${\sigma }_n$ is the standard deviation of $o_n$. The corresponding $e_{eq}^2$ is calculated for the nonweighted multimodel T-year return level. Note that, for each n, $o_n$ and ${\sigma }_n$ are unchanged between $e_w^2$ and $e_{eq}^2$ because those are computed from the observations on grid point n. Then, the error index is obtained by:

$$I^2=\frac{e_w^2}{e_{eq}^2}.$$

(12)

This index can be useful to evaluate if the weighted mean is improved compared to the simple average [40]. The smaller the $I^2$, the larger the improvement in the weighted average compared to the simple average for the target variable (AMP1) for the reference period.

The computed values of $I^2$ in this study are 0.622 for PI-weights and 0.627 for P-weights for the 20-year return level. Thus, there are 37.8% and 37.3% improvements in the PI-weighted and the P-weighted averages, respectively, compared to the nonweighted average. The decrease of the error index by adding I-weights from P-weights is 0.005, which is a 0.80% decrease relative to that of the P-weights. The results for the 50-year return level are similar to these.

8.2. Prediction Variance for the Future

The second measure to compare is from the weighted variance formula [71,72]:

$$Var\left(r\left(T\right)\right)=\sum _{k=1}^K{[r_k\left(T\right)-\overline{r}\left(T\right)]}^2{w}_k+\sum _{k=1}^{K}Var\left(r_k\left(T\right)\right)w_k,$$

(13)

where $r\left(T\right)$ is the T-year return level, $r_k\left(T\right)$ is the T-year return level from the k-th model, and $\overline{r}\left(T\right)=\frac{1}{K}{\sum }_{k=1}^K{r}_k\left(T\right)w_k$. The first term of (13) is the among-model variance, and the second one is within-model variance. $Var\left(r_k\left(T\right)\right)$ is estimated using a bootstrap technique in this study. Based on this variance, we quantify the skill of the PI-weighted method with the Brier skill score [54]:

$$BSS=100\times \{1-\frac{Va{r}_{PI}}{Va{r}_{Eq}}\},$$

(14)

where $Va{r}_{PI}$ is the variance of the T-year return level calculated by (13) from the PI-weights, and $Va{r}_{Eq}$ is the variance calculated by (13) from the equal weights for each grid point. BSS is a relative quantity and shows a percentage improvement of the PI-weighted method over the simple average.

Figure 14 shows the differences of variances of 20-year return levels calculated from the PI-weights and the equal weights plotted for 46 grid points for each future period and for each scenario. $Va{r}_{PI}$ for the Y-axes is the variance calculated by Equation (13) from the PI-weights. $Va{r}_{Eq}$ for the X-axes is the corresponding variance calculated from the equal weights. In the top panel, more points are located below the diagonal lines, which means that the variances calculated from the PI-weights are smaller more often than those from the simple average. The middle (bottom) panel shows the differences of between-model (within-model) variances. $Va{r}_{BM,PI}$ and $Va{r}_{WM,PI}$ for the Y-axis are between-model and within-model variances from the PI-weights, respectively, and $Va{r}_{BM,Eq}$ and $Va{r}_{WM,Eq}$ for the X-axis are the corresponding variances calculated from the equal weights. The within-model variances from both weights are similar, as seen in the bottom panel. Thus, we know that the differences in total variances in the top panel are mainly because of the differences of between-model variances.

Table 4. The averaged values of the Brier skill scores (BSSs), averaged differences of variances (unit: mm${}^2$) of the 20-year return level over 46 grid points in the Korean peninsula, and relative improvements (RIs) of the PI-weighted method over the simple average (last column) for the three periods and for the three scenarios. RI$=\frac{V_{Eq}-V_{PI}}{V_{PI}}\times 100$. $DV=Va{r}_{Eq}-Va{r}_{PI}$, $D{V}_{BM}=Va{r}_{BM,Eq}-Va{r}_{BM,PI}$, and $D{V}_{WM}=Va{r}_{WM,Eq}-Va{r}_{WM,PI}$. Other acronyms are the same as in Figure 14.

ssp	p	BSS(%)	${\mathit{BSS}}_{\mathit{BM}}$	${\mathit{BSS}}_{\mathit{WM}}$	DV	${\mathit{DV}}_{\mathit{BM}}$	${\mathit{DV}}_{\mathit{WM}}$	RI(%)
	P1	−2.2	−12.3	2.0	−50	−77	27	−3.3
SSP2-4.5	P2	7.7	19.7	3.2	129	82	47	8.5
	P3	13.4	39.0	2.7	309	267	42	17.1
	P1	17.9	42.9	5.4	300	238	62	24.0
SSP3-7.0	P2	11.0	22.5	4.7	228	153	75.0	13.7
	P3	16.7	36.2	1.4	426	378	48	19.5
	P1	6.4	24.7	−4.5	113	138	−25	8.0
SSP5-8.5	P2	10.4	33.9	−6.9	178	291	−113	8.0
	P3	22.5	42.2	−1.1	826	834	−8	26.0

reads the averaged values of BSS over 46 grid points, the averaged differences of variances (unit: mm${}^2$) of the 20-year return level, and relative improvements (RIs) for three periods and for three scenarios. The BSS values range from −2.2% to 22.5%. The averages of BSS over all periods are 6.29% for SSP2, 15.2% for SSP3, and 13.1% for SSP5 (11.5% overall). These averages of BSS also increase as the period moves from P1 to P3. The BSS due to between-model variance is bigger than that due to within-model variance. Watching the sixth to the eighth columns of this table, we know that the positive differences of two variances ($Va{r}_{Eq}-Va{r}_{PI}$) are mainly because of the differences of two between-model variances ($Va{r}_{BM,PI}-Va{r}_{BM,Eq}$). These patterns are prominent for the SSP5-8.5 scenario with increasing values from P1 to P3. The last column shows the RI of the PI-weighted method over the simple average in the variances. The RI values range from −3.3% to 26.0%. The averages of RI over all periods are 7.45% for SSP2, 19.1% for SSP3, and 14.0% for SSP5 (13.5% overall). These averages of RI also increase as the period moves from P1 to P3.

9. Discussion

It is generally accepted that increasing greenhouse gases induce atmospheric temperature warming, which results in increasing equivalent potential temperatures and specific humidity, according to the Clausius–Clapeyron relationship [73]. The increase in atmospheric water vapor is the main factor in generating convective instability. In view of this insight, Kim et al. [26] obtained an indication that increasing the extreme precipitation over South Korea in the past and future scenarios is related more to a change in convective instability rather than to synoptic conditions. Another plausible explanation for the increase in the maximum precipitation over the Korean peninsula is the increase in both the frequency and strength of typhoons in the region. Typhoons have a greater influence on the southern area than the northern part. The occurrence of extreme rainfall in the southern part is related to more sources or variables than in the northern area [25,26].

We observed a scale discrepancy between the model grids and the observation stations, which may compromise the credibility of our results. Regional projections require a fine resolution of the simulation models, whereas some of the CMIP6 models have a coarse resolution. It is known that the model spread is one of the major sources of uncertainty in regional predictions (Hawkins and Sutton, 2009). However, large numbers of simulation models can reduce the uncertainty and provide a robust projection [74,75]. A total of twenty-one CMIP6 models with different scales used in this study can cover the study region with a finer spatial structure to increase the reliability of the projection.

The daily precipitation data consist of measurements from 00:00 to 24:00 throughout the day. In daily observations, the rainfall does not accumulate between 22:00 and 02:00. In the data used in this study, such precipitation is divided and recorded in two separate days. The actual serious daily risk due to heavy rainfall does not exactly depend on the precipitation over the time duration from 00:00 to 24:00 exclusively. It is therefore recommended to consider the AMP1 data based on the maximum precipitation during the 24 h movement. In this sense, the results presented in this study underestimate the actual intensity and frequency of AMP1. More realistic daily data, such those as obtained after moving for 24 h and the annual maximum of two (and several) days of precipitation, should be used in a future study for assessment of risk owing to extreme rainfall.

Brunner et al. [36,61] considered multiple observational (or reanalysis) datasets to include an estimate of the observational uncertainty. They proposed a novel approach to account for the observational spread and uncertainty in a multi-model weighting study, which can lead to robust results and a more precise uncertainty quantification. In addition, considering multiple observational datasets may address the problem in which the BC and performance-based weighting scheme utilize an excessive number of observations. We believe that using the observations twice in the BC and weight calculation is unadvisable. Xu et al. [35] considered a Bayesian weighting method that removes observations during the initial phase of the downscaling and adds them in the estimation of the posterior distribution. However, if the series of observations is sufficiently long to divide into two parts, we may use one part for the BC and the other part for the weight calculation. Although we did not apply these methods, this would be a good approach in a future study.

In determining the shape parameter ${\sigma }_D$ for the performance weights, we used the standardized return levels. However, as a different approach, the BMA method described in [27,41,72] can be employed. The BMA weighting does not need to standardize the return levels or determine the shape parameter ${\sigma }_D$. It requires a bootstrap to estimate the variances of the return levels, which is straightforward in current computing facilities. Thus, calculating the P-weights based on the BMA, with an effort to minimize the probability of overconfidence, is recommendable.

10. Summary

We estimated the future changes in precipitation extremes within the Korean peninsula using observations, 21 multiple CMIP6 models, generalized extreme value distribution, the multivariate bias correction technique, and the model weighting method (PI-weighting), which account for both the performance and independence of the models. To avoid overfitting in the PI-weighting, we considered five climate variables for East Asia.

In applying the PI-weighting method, we suggest two ways of selecting two shape parameters, based on the p-value of the chi-square statistic and entropy. The suggested methods are simple and intuitively appealing, although they may need more justification for use in other studies.

From the analysis described this study, we realized that a 1-in-20 year (1-in-50 year) annual maximum daily precipitation within the Korean peninsula will likely become a 1-in-11 (1-in-26) year, a 1-in-8 (1-in-20) year, and a 1-in-7 (1-in-15) year event in terms of the median by the end of the 21st century under the SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios, respectively, as compared to the observations from 1973 through 2010. These results are similar to or less frequent than those obtained by previous studies [23,27,41], but still, they predict more frequent and intensified extreme precipitation events by the end of the 21st century as compared to 1973 through 2010.

The expected frequency of the reoccurring years, particularly for AMP1 from 100 to 300 mm under the SSP5 scenario, is projected to increase by approximately 1.0 times that of the past 30 years for period 1 (2021–2050), approximately 1.25 times that for period 2 (2046–2075), and approximately 1.91 times that for period 3 (2071–2100).

From the analysis based on latitude, we found that extreme rainfall is more prominent in the southern and central parts of the peninsula. The downpour in the southern part is approximately 1.54 times heavier than that of the northern part and approximately 1.72 times that of the spatial median of the Korean peninsula. For example, for 150 mm of AMP1, the expected waiting times until reoccurrence in the southern part (spatial median of the peninsula) are 2.6 (4.4) years during the reference period, 2.1 (3.6) years during P1, 1.7 (3.0) years during P2, and 1.5 (2.5) years during P3 based on the SSP5-8.5 scenario.

Heavy rainfall can have a significant effect on human life, communities, infrastructure, agriculture and livestock, and natural ecosystems. Thus, in addressing the impact of climate change due to more frequent extreme precipitation events, governments and communities should prepare the proper infrastructure and systems more carefully and securely to prevent critical damage, such as loss of life from landslides and flooding.