Temporal and Spatial Characteristics of Multidimensional Extreme Precipitation Indicators: A Case Study in the Loess Plateau, China

Abstract

Extreme precipitation can seriously affect the ecological environment, agriculture, human safety, and property resilience. A full-scale and scientific assessment in extreme precipitation characteristics is necessary for water resources management and providing decision-making support to mitigate the potential losses brought by extreme precipitation. In the present study, a multidimensional risk assessment framework is developed to investigate the spatial–temporal changes in different extreme precipitation indicators. The Gaussian mixture model (GMM) is applied to fit the distribution for each indicator and carry out single index risk assessment. The joint probabilistic features of multiple extreme indicators can be explored through coupling the GMM distributions into copulas. In addition, the moving window approach and the Mann–Kendall test are integrated to examine non-stationary risks (evaluated by “AND”, “OR”, and Kendall return periods) of multidimensional indicators along with their changing trends and significance. The proposed assessment framework is applied to the Loess Plateau, China. Four extreme precipitation indicators are characterized: the amount (P95), the number of days (D95), the intensity (I95), and the proportion (R95) of extreme precipitation. The spatial–temporal changes of these indicators and their multidimensional combinations (including six two-dimensional and three three-dimensional combinations) are fully identified and quantitatively evaluated.

1. Introduction

In recent years, the increasingly catastrophic impact of hydrometeorological disasters has aroused public concern over extreme events due to climate change [1,2,3]. The drastic changes in the global water cycle may lead to increasing extremes such as extreme precipitation and drought in local areas [4,5]. The extreme precipitation can seriously affect eco-environment, agriculture, human safety, and property resilience [6,7,8]. Risk analysis for extreme precipitation is important for hydrologic designs and environmental systems management [9,10,11]. In the process of climate change, the spatial–temporal distribution of precipitation has also changed significantly in many regions around the world. The characteristics of extreme precipitation have also presented dynamic features [12], in which the changes of magnitude or/and frequency of extreme precipitation may have significant impacts on local flooding [10,13]. However, the commonly used risk assessment techniques with stationary assumption (constant probability) make the resulting hydrologic design questionable under changing environmental conditions [10]. Therefore, it is necessary to conduct an in-depth study on the non-stationary risk characteristics of extreme precipitation under climate change. Such research is remarkably valuable to enhance scientific understanding of the possible changes in floods and other secondary disasters, and provide key information for generating appropriate adaptation measures.

A variety of indicators are employed to characterize the duration, intensity, and amount of extreme precipitation [14,15,16,17]. Frequency characteristics estimated from probability distribution functions are often used to assess risks of those extreme precipitation indicators [18,19]. For example, based on the generalized extreme value theory [20] and the moving window approach, Ghosh et al. [21] have explored the non-stationary frequency characteristics of annual maximum rainfall in India. It is worthy to notice that the extreme precipitation indexes are usually correlated among each other, and univariate estimation of the extreme precipitation indicators may underestimate or overestimate the severity of the extreme event [11,15]. The joint probabilistic characteristics of correlated indicators of extreme precipitation are of great importance in hydrological risk management, and have also become essential factors under consideration for water resources planning, design, and management [11].

As a widely accepted and applied multidimensional joint analysis method, the copula method has been successfully applied in hydrological risk assessment [22,23,24]. Two prominent advantages of this method are that the marginal and joint distributions for indicators can be built separately, giving flexibility in choosing different marginal and dependence models, and there are a number of marginal distributions and copula functions available. Particularly, the bivariate copula functions have already been introduced into two-dimensional extreme precipitation studies, and different kinds of return periods (e.g., “AND”, “OR” return periods, see Section 3) are applied to analyze the joint risk characteristics of extreme precipitation indicators (see, e.g., in [15,25]). For instance, Liu et al. [26] have used the copula method to study the joint probabilistic features of six extreme precipitation indicators in southwestern China. Based on projections from climate models, Goswami et al. [15] have applied the copula method to assess the bivariate behaviors of seven extreme indicators in the North Sikkim in the future. Li et al. [11] have incorporated non-stationary analysis into a bivariate analysis of extreme precipitation by using the Archimedean copulas combined with a moving window scheme. The Akaike’s information criteria (AIC, [27]) and Kolomogorov–Smirnov test (KS) were usually employed to do the goodness-of-fit assessment for the copulas [15,16].

Although joint risk analysis of extreme precipitation indicators has been carried out by numerous studies, some challenges still need to be further addressed: (i) Different marginal distributions have been used in previous studies, which may produce noticeable uncertainty, and thus it is desired to adopt a widely applicable marginal distribution to eliminate or reduce model uncertainty. (ii) Bivariate behaviors of extreme precipitation indicators have been extensively characterized, but there are few studies on high-dimensional extremes larger than two, and also no comprehensive analysis on combinations of multiple-dimensional indicators. (iii) The joint return periods are frequently characterized by “AND” and “OR” types, which can hardly meet the needs of multivariate risk assessments in practical application. Moreover, no research in this field (extreme precipitation) is found to address the high-dimensional Kendall return period. Further, the vine copula, which is built by decomposed bivariate copulas, has been used to construct multivariate joint distributions in recent years [28]. However, for the three-dimensional cases, the parametric copula functions, such as Archimedean and elliptical copula families, are easy to calculate and are effective to capture the dependence structure among the variables [29,30].

The challenges mentioned above prompt us to propose a comprehensive risk assessment framework, which can extensively evaluate the spatial–temporal changes of multidimensional extreme precipitation indicators. Therefore, the main purpose of this study is to achieve the following objectives. (i) Adopting a marginal distribution which can be applied uniformly for all indicators. (ii) Assessing the risks of multidimensional indicator combinations by different types of return periods, including two- and three-dimensional Kendall return periods. (iii) Establishing a non-stationary assessment model and evaluate the spatial–temporal changing features of risks for multidimensional indicators.

2. Overview of the study area and data

With an area of ~640,000 km², the Loess Plateau, which is located in the north of central China, is selected as the study area (Figure 1). The altitude of the Loess Plateau mainly ranges between 900 to 2200 m, and the terrain is characterized by high-west and low-east patterns. In the Loess Plateau, there are some mountains up to 3500 m above sea level, as well as plains and basins at an altitude of only 400 to 500 m. The Loess Plateau is also the largest area covered by loess in the world, presenting loose soil structure and fragile eco-environmental systems.

The average annual precipitation in the Loess Plateau is only 413.6 mm. Due to the influence of monsoon climate and complex geomorphic types, the precipitation in the Loess Plateau varies greatly in time and space. The spatial distribution shows more precipitation in the southeast and less in the northwest (Figure 1b). From southeast to northwest, the climate successively presents semi-humid, semi-arid, and arid types. The temporal distribution of precipitation in the Loess Plateau is also uneven. This mainly contains two aspects: First, the inter-annual variability of precipitation is noticeable. In the years with excessive rainfall, the precipitation may be several times more than that in the water-deficient years. Second, seasonal differences in precipitation are also large, with precipitation in Summer and Autumn accounting for 60–80% of the annual precipitation. The precipitation in rainy season may be concentrated in few intensive rainfall events. It is customary that a single rainstorm accounts for 30% or more of the annual rainfall, leading to coexistence of water shortage and soil erosion in the Loess Plateau.

The daily precipitation data with time period from 1959 to 2018 of 63 stations (Figure 1b, two of them are located near the Loess Plateau) are collected from the China Meteorological Administration. The ratio of missing values over all stations is from 1.23% to 3.15%, with an average of 1.83%, and the proportion of missing values in rainy season (July–September) is 0.44% on average. The missing values were replaced by average precipitation of neighboring stations, and the calculation of extreme indicators would be affected very little [31]. In order to facilitate the comparative study of the multidimensional risks of extreme precipitation indicators in different climatic regions, the study area is divided into two sub-regions based on the rainfall pattern (400 mm precipitation contour), taking into account the integrity of the basin and the similarity of other natural geographic factors. The first sub-region is the arid and semi-arid area (Region 1), including 25 stations, and the second area (Region 2) contains 38 stations, experiencing a semi-humid climate (Figure 1b).

3. Method

3.1. Gaussian Mixture Model

To quantitatively evaluate the probabilistic characteristics of each extreme precipitation indicator and meet the needs of fitting joint distributions for different indicator combinations, suitable marginal distributions should be constructed in advance. In order to reduce the uncertainty brought by differences in marginal distributions, the Gaussian mixture model (GMM) [22] is applied to fit the marginal distributions for all extreme precipitation indicators in this study. The probability density function (PDF) of GMM is expressed by a weighted sum of n-component Gaussian PDFs:

$$p(x)={\displaystyle \sum _{i=1}^n{\alpha }_i{G}_i(x;{\mu }_i,{\sigma }_i)}$$

(1)

where x denotes one extreme precipitation indicator, α_i (i = 1, 2,…, n) are the non-negative mixture weights which satisfy ${\displaystyle \sum _{i=1}^n{\alpha }_i=1}$, G_i(x; μ_i, σ_i) represents the ith Gaussian density component, and can be expressed as

$$G_i(x;{\mu }_i,{\sigma }_i)=\frac{1}{\sqrt{2\pi }{\sigma }_i}\exp (\frac{-{(x-{\mu }_i)}^2}{2{\sigma }_j^2})$$

(2)

where μ_i, and σ_i represent the mean and standard deviation of the ith Gaussian model, respectively. Given an appropriate number of components, the GMM can fully approximate a wide variety of distributions encountered in practical applications [22,32].

The maximum likelihood estimation (MLE) method is used to estimate the unknown parameters in the GMM model by maximizing the likelihood function through the Expectation–Maximization algorithm (EM) [28]. Let θ_i = (α_i, μ_i, σ_i), Θ = (θ₁, θ₂, …, θ_n), and the indicator x has a sample size of m, the likelihood function Լ of the model with n components can be expressed as

$$L(x|\Theta )=\log \prod _{j=1}^m{\displaystyle \sum _{i=1}^n{G}_i(x;{\mu }_i,{\sigma }_i)}={\displaystyle \sum _{j=1}^m\log {\displaystyle \sum _{i=1}^n{G}_i(x;{\mu }_i,{\sigma }_i)}}$$

(3)

The EM algorithm uses an iterative approach to maximize the likelihood function, and the parameters in GMM are gradually updated. The EM algorithm consists of two steps: E-step and M-step. A brief description of the two steps can be shown as follows. E-step: Calculate the posterior probability of the ith component based on the current estimated parameters or the initial given parameters:

$${\beta }_{i,j}=E({\alpha }_i|x_j;\Theta )=\frac{{\alpha }_i{G}_i(x_j;\Theta )}{{\displaystyle \sum {}_{i=1}^n{\alpha }_i{G}_i(x_j;\Theta )}},1\le i\le n,1\le j\le m$$

(4)

M-step: Update all the parameters in the model according to the weighted averages over all elements of the indicator:

$${\alpha }_i^{\prime }=\frac{{\displaystyle {\sum }_j^m{\beta }_{i,j}}}{m}$$

(5)

$${\mu }_i^{\prime }=\frac{{\displaystyle {\sum }_j^m{\beta }_{i,j}{x}_i}}{{\displaystyle {\sum }_j^m{\beta }_{i,j}}}$$

(6)

$${\sigma }_i^{\prime }=\frac{{\displaystyle {\sum }_j^m{\beta }_{i,j}{(x_i-{\mu }_i^{\prime })}^2}}{{\displaystyle {\sum }_j^m{\beta }_{i,j}}}$$

(7)

More detailed information regarding the EM algorithm can be seen in [33]. In this study, the upper limit of the iteration times is set to 100,000, and the convergence condition is the difference between the update parameters and the last parameters, which is set to be less than 10⁻⁵. The goodness-of-fit for GMM is examined using the root mean square error (RMSE), Akaike’s information criteria (AIC [27]), and the Kolomogorov–Smirnov (KS) test (at the 5% significance level of the p-value) [30].

3.2. Joint Probability Distribution and Return Period

3.2.1. Joint Probability Distribution for Extreme Precipitation Indicators

Let X = (X₁, X₂, …, X_d) be a d-dimensional random vector formed by d extreme precipitation indicators. Based on the theory of copula proposed by Sklar [34], the d-dimensional distribution function F_X can be expressed by the Copula function as

$$F(X)=C(F_{X_1}(x_1),F_{X_2}(x_2),\cdot \cdot \cdot ,F_{X_d}(x_d))$$

(8)

where C represents the copula function that satisfying [0,1]^d → [0,1]; $F_{X_1}(x_1),F_{X_2}(x_2),\cdot \cdot \cdot ,F_{X_d}(x_d)$ denote continuous marginal distribution functions of the d-dimensional random vector X.

Various copula functions, which can be classified into different copula families (e.g., elliptical, extreme value, and Archimedean), are available in multivariate analysis. In this study, four candidate copulas at two- and three-dimensional forms are applied to analyze the joint probability distributions for two- and three-dimensional indicator combinations. These copulas are the Student t and Gaussian copulas from elliptical family, and Frank and Gumbel copulas from Arachmedean family. These copulas are adopted because these parametric copula functions can be easily constructed for two- or higher-dimensional situations, and be effective to capture a broad class of dependence structures [29,30]. The basic characteristics of those copulas can be seen in Xu et al. [29]. The unknown parameters in the chosen copulas are estimated based on the inversion of Kendall’s τ [35,36], by using the copula package in R [29]. The validation passes when the p-value in the KS test is greater than 0.05, and smaller AIC and RMSE values imply a better goodness-of-fit. Detailed descriptions on these testing methods can be found in Zhang [37].

3.2.2. Return Periods for Extreme Precipitation Indicators and Their Combinations

The concept of return period is defined as the average time interval between successive occurrence of a given event, which can be estimated from the exceedance probability. The return period $T_X$ for a single index (represented by X) beyond a certain value x can be formulated as

$$T_X=\frac{1}{P(X>x)}=\frac{1}{1-F(x)}$$

(9)

Three types of return periods have been proposed the two- and three-dimensional cases in hydrology, which are “AND”, “OR”, and Kendall return periods. Figure 2 illustrates the difference between the three return periods in a two-dimensional case.

In Figure 2, u and v represent the marginal distributions of two indicators X₁ and X₂ (F(X₁) = u, F(X₂) = v), respectively. C(u_i, v_i) = t is the joint probability for a given (u_i, v_i). The three return periods are defined based on three different types of probability expressions. The space formed by the dotted lines, u = 1 and v = 1 in the upper right corner indicates the probability that both the two indicators exceed their threshold values P{u> u_i ∩ v> v_i}. The return period $T_{X_1,X_2}^{and}$ corresponding to the simultaneous occurrence of two events can be expressed as

$$\begin{array}{ll}{T}_{X_1,X_2}^{and}& =\frac{1}{P(X_1>x_1\cap X_2>x_2)}=\frac{1}{1-F(x_1)-F(x_2)+F(x_1,x_2)}\\ & =\frac{1}{1-u-v+C(u,v)}\end{array}$$

(10)

For a three-dimensional case, in which three indicators exceeding their thresholds simultaneously, the joint return period $T_{X_1,X_2,X_3}^{and}$ is shown as

$$\begin{array}{l}{T}_{X_1,X_2,X_3}^{and}=\frac{1}{P(X_1>x_1\cap X_2>x_2\cap X_3>x_3)}\\ =\frac{1}{1-F(x_1)-F(x_2)-F(x_3)+F(x_1,x_2)+F(x_1,x_3)+F(x_2,x_3)-F(x_1,x_2,x_3)}\\ =\frac{1}{1-u-v-w+C(u,v)+C(u,w)+C(v,w)-C(u,v,w)}\end{array}$$

(11)

where X₁, X₂, and X₃ represent extreme indicators; x₁, x₂, and x₃ are the certain thresholds of the indicators; F(x), F(x₁,x₂), …, F(x₁,x₂,x₃) denote the cumulative probability distributions (CDFs) or joint CDFs corresponding to the thresholds; u, v, and w represent the marginal distributions of the indicators; and C(u,v), …, C(u,v,w) stand for copula functions describing joint distributions of extreme precipitation indicators.

The remaining space excluding the shadow part in Figure 2 indicates the probability that at least one of the two indicators exceeds the threshold (P{u > u_i ∪ v > v_i}). The corresponding return period expressed by $T_{X_1,X_2}^{or}$ can be calculated as

$$T_{X_1,X_2}^{or}=\frac{1}{P(X_1>x_1\cup X_2>x_2)}=\frac{1}{1-F(x_1,x_2)}=\frac{1}{1-C(u,v)}$$

(12)

For the three-dimensional case, in which at least one of the three indicators exceeding the threshold, the return period $T_{X_1,X_2,X_3}^{or}$ can be expressed as

$$\begin{array}{l}{T}_{X_1,X_2,X_3}^{or}=\frac{1}{P(X_1>x_1\cup X_2>x_2\cup X_3>x_3)}\\ =\frac{1}{1-F(x_1,x_2,x_3)}=\frac{1}{1-C(u,v,w)}\end{array}$$

(13)

The space formed by the curve C = t, u = 1 and v = 1 (C = t is a surface in the three-dimensional case) represents the probability P(k) = P(C > t), and the corresponding Kendall return period is 1/P(k). The Kendall return period is calculated by introducing the Kendall measure Kc [38]:

$$Kc(t)=\mathrm{P}(C(u_1,u_2,\dots ,u_d)\le t)$$

(14)

For bivariate Archimedean copulas, Kc can be calculated from the generators of the copula functions. For the t copula, Gaussian copula and Archimedean copulas of three or higher dimensions, the Kc can hardly be calculated directly. In this study, the random sampling algorithm (10,000 sampling in this study) is used to estimate Kc [39], and then Kendall return periods of all two-dimensional and three-dimensional cases are obtained. Sampling algorithm steps to estimate Kc:

• Simulate a sample u₁, ..., u_m from the d-dimensional copula C;
• For i = 1 ,..., m calculate w_i = C(u_i);
• Estimate Kc: ${\widehat{K}}_C(t)=\frac{1}{m}{\displaystyle {\sum }_{i=1}^{m}1(w_i\le t)}$.

4. Case Study and Results Analysis

4.1. Case Study

In order to quantitatively assess the changes in extreme precipitation from different perspectives, a number indicators focusing on various aspects of extreme precipitation have been proposed. According to the recommendations of the World Meteorological Organization, other relevant studies [14,15], four extreme precipitation indicators are used in this study to describe extreme precipitation features in Loess Plateau, including, the amount (P95), the number of days (D95), the intensity (I95), and the proportion (R95) of extreme precipitation. Detailed definitions of these indicators are shown in

Table 1. Concept of four extreme precipitation indicators.

Indices	Abbreviations	Definitions	Unit
Number of extreme precipitation days	D95	Number of days with P > 95th percentile (daily precipitation exceeding the 95th percentile of precipitation series during 1971–2000).	days
The amount of extreme heavy precipitation	P95	Annual total amount of precipitation with P > 95th percentile	mm
The intensity of extreme precipitation	I95	Average daily precipitation intensity of extreme precipitation	mm/day
Ratio of extreme precipitation	R95	Ratio of annual total precipitation due to events exceeding the 95th percentile	-

.

Six two-dimensional and three three-dimensional combination schemes of the four indicators can be formulated. The indicator combination schemes are shown in

Table 2. The two-dimensional and three-dimensional indicator combination schemes.

ID	Combinations	ID	Combinations	ID	Combinations
1	{D95, P95}	4	{P95, I95}	7	{D95, P95, R95}
2	{D95, I95}	5	{P95, R95}	8	{D95, I95, R95}
3	{D95, R95}	6	{I95, R95}	9	{P95, I95, R95}

, and the joint return period of “AND”, “OR”, and Kendall for these combinations are characterized. Taking “AND” return periods of the combination events of {D95, P95} and {D95, P95, R95} as examples, the combination {D95, P95} indicates D95 and P95 simultaneously exceed their thresholds (Equation (10)), which means that an extreme precipitation event with long time and large amount would happen and severe flood may occur. The combination {D95, P95, R95} indicates that an extreme precipitation event with long duration and large amount occurs, and the proportion of this event to annual precipitation also exceeds the threshold (Equation (11)). For the reason that the annual rainfall is largely completed by extreme precipitation in this situation, adaptation strategies should be devised to achieve reasonable storage of water resources in the flood season and allocation of water resources throughout the year. The meaning of the other types of return periods for each joint event can be inferred from the definitions in Section 3.2.2.

The spatial risk characteristics for single indicators and multi-indicator combinations throughout the Loess Plateau can be obtained from the results over all stations in the region. In order to explore the trend in different extreme precipitation indicators and their joint risks, a 30-year moving window approach is applied. From the beginning of 1959–1988 to the end of 1989–2018, 31 groups of results for all stations were obtained. For each group, every indicator is fitted with marginal distribution (the GMM) and copulas are fitted for all two and three-dimensional combinations. The initial group (from 1959 to 1988) and the final group (from 1989 to 2018) of the 31 results belong to two independent time periods. The comparison of the results for these two groups can directly reflect the changes of extreme precipitation in the past two representative periods. By evaluating the outcomes of 31 groups in time series, the trend for both the univariate and multivariate risks can then be obtained. The Mann–Kendall (MK) test [40] is applied to evaluate whether the trend is significant, and the linear regression method is used to quantitative assess the degree of change. The annual extreme precipitation indicators are used in this study, and the Mann–Kendall test is basically not affected by serial correlation.

Taking Yulin station (Station 34 in Figure 1) as an example, the whole calculation process for the three indicators of P95, D95, and R95 during the period 1959–1988 is illustrated. GMM is used to construct the marginal distributions for the three indicators, and the test results of KS, RMSE, and AIC are shown in

Table 3. Statistical test results for Gaussian mixture model (GMM) and best-fitted copulas at Yulin Station.

Scheme	Distribution	KS Test		RMSE	AIC
		$\mathit{T}$	p-Value
D95	GMM	0.143	0.568	0.055	−170.089
P95	GMM	0.089	0.972	0.029	−208.096
R95	GMM	0.082	0.987	0.032	−202.164
{D95, P95}	Gaussian	0.022	0.615	0.027	−215.108
{D95, R95}	t	0.024	0.697	0.032	−203.671
{P95, R95}	Gaussian	0.028	0.790	0.030	−208.887
{D95, P95, R95}	t	0.032	0.648	0.021	−228.451

KS test: the Kolomogorov-Smirnov test; $T$: the value of the test statistic; p-Value: the significance level; RMSE: root mean square error; AIC: Akaike’s information criteria.

. Table S1 in the Supplementary Materials also shows the statistical test comparisons between GMM and the best fitted ones of five other commonly used parametric distributions (Gamma, Pearson type III (P III), lognormal (LN), Log Pearson type III (LP III), and generalized extreme value (GEV)) [23]. The p-values of KS tests are all greater and the AIC values are all smaller for GMM than other parametric distributions. Figure S1 shows the PDF fitting results of GMM with other parameter distributions in Table S1. It can be seen that even in the case of non-multimodal PDF, the performance of GMM can be similar to the optimal parameter distribution (Figure S1c) or even better (Figure S1b). Candidate copulas of Frank Gumbel, Gaussian, and t copulas (at two- or three-dimensional forms) are applied to build the joint probability distributions of multi-indicator combinations. The optimal copulas selected according to the KS, RMSE and AIC tests are shown in Table 3, and the corresponding test values are also displayed. It can be seen that the p-values of KS test for all indicators and multi-indicator combinations are all greater than 0.05, which means that the fitted distributions can effectively describe the probability characteristics of the indicators in various dimensions. The test results of the other copulas such as Frank and Gumbel are shown in Table S2 in the Supplementary Materials. Figure 3a–f further shows the goodness-of-fit by comparing the empirical PDFs and CDFs with those from GMM distributions of the three indicators, and the theoretical distributions show good agreements with the empirical distributions. Figure 3g–i shows the bivariate CDFs for three two-dimensional indicator combinations.

Given the 10-year univariate return period of three indicators, joint return periods of “AND”, “OR”, and Kendall for two- and three-dimensional cases are then evaluated. It can be seen from

Table 4. Joint returned periods of “AND” (T_and), “OR” (T_or), and Kendall (T_k) for two and three-dimensional cases.

Return Periods	{D95, P95}	{D95, R95}	{P95, R95}	{D95, P95, R95}
Tor	8.09	7.28	7.79	7.01
Tk	11.66	13.50	11.77	12.66
Tand	13.09	15.99	13.96	17.45

that, for the same indicator combination scheme, the ranking order of the three return periods is T_or < T_k (Kendall return period) < T_and. This can also be inferred from the corresponding probabilities of various return periods in Figure 2, that is, P(and) < P(k) < P(or). Compared to the two-dimensional cases, it can be found that the T_or of the three-dimensional return period is smaller. The maximum return period (17.45 years) can be found when the three indicators occur simultaneously.

4.2. Changes of Different Schemes in Two 30-Year Stages

Based on the fitted GMM distributions, for a given return period, the values of the corresponding indicators can be obtained through the univariate return period formula (Equation (9)). In order to conduct a more illustrative comparison of the 10-year return period values between different regions, the values of four indicators at all stations under two 30-year periods are interpolated to the whole region by using the Kriging method (Figure 4). Figure 4a–d provides the comparisons of four indicator values between the time period of 1959–1988 (T1 period, left ones) and time period of 1989–2018 (T2 period, right ones).

The results show that the univariate return period values of the 10-year stage for all groups are similar in spatial distribution between the two 30-year stages over the entire region. The minimum values of D95 and P95 appear in the northwest of Loess Plateau, and there is a stepwise increasing trend from northwest to southeast for these two indicators. In the southeastern area, it is more prone to heavy rainfall for a long time (D95), which brings greater amounts of precipitation (P95). The values of I95 exhibit obvious increasing direction from west to east, while R95 presents a clear North–South gap. The mean values (averaged over the stations) of D95, P95, and I95 in Region 1 are 0.67 days, 78.2 mm, and 15.6 mm/day smaller than those in Region 2 for T1 period. Also, these indicators in Region 1 are 0.63 days, 81.2 mm, and 16.7 mm/day smaller than those in Region 2 for T2 period. The mean value of R95 in Region 1 is slightly higher than that in Region 2, and the values of R95 could be above 0.5 in the north part of Region 1. This indicates that, with the least annual precipitation in the northwest, nearly half of the annual precipitation would fall by extreme heavy precipitation.

The spatial distributions of the indicator values behave differently between the two time periods. The 10-year return period values of the four indicators in the T2 period (1989–2018), when compared with those in T1 period (1959–1988), exhibit the following changes. The value of D95 decreases slightly throughout the region, with an average decrease of 0.26 days. The station with the largest decrease of D95 appears in the northeast with a decrease of 1.95 days, whereas the station with the largest growth of D95 appears in the southwest region with an increase of 1.2 days. Similar to the changes of D95, P95 decreases significantly in north and southeast, with an average decrease of 13.13 mm. I95 increases significantly in some stations in the southwest and southeast regions with the largest increase of 9.96 and 25.4 mm/day, respectively. The spatial difference of R95 in different time periods is not obvious, and the change of annual precipitation does not bring obvious change of R95.

Table 5. Statistics of changes in four indicators in different regions.

Region	Number of Stations	D95		P95		I95		R95
		A	N	A	N	A	N	A	N
1	25	−0.24	7	−11.29	11	0.36	14	−0.013	11
2	38	−0.27	14	−8.27	14	1.44	24	−0.004	16

A: The average change over all stations in the region; N: Number of stations with increased value.

displays the average changes of four indicators in different regions and the number of stations which show growth changes by comparing with that in T1 period.

Under the condition that the indicator values at a 10-year return period, the joint return periods of T_and for two-dimensional or three-dimensional indicator combinations can then be obtained. Figure 5 shows T_and at each station in T1 period (i.e., 1959–1988) (sub-graphs on the left) of six two-dimensional indicator combinations, and the spatial changes in T2 (i.e., 1989–2018) compared with that in the T1 period (right ones of each group). The red areas in the right sub-graphs indicate decreasing T_and in past 30 years, which signify that the regions are more prone to joint risk than the T1 period.

The averaged T_and of three two-dimensional combination events {D95, P95}, {D95, R95}, and {P95, R95} in the T1 period are 13.05, 14.68, and 13.35 years, respectively. The T_and values are only slightly longer than the return period of a single event (10 years), and the minimum values of the three combined events are 11.03, 12.25, and 11.31 years, respectively. It is proved that these three groups of combination events are extremely easy to happen, and the events of different indicators in each combination at a 10-year return period level usually occur together. There is no obvious regularity in the spatial distributions of T_and for these three combinations. The spatial distributions of the changes in T_and for the three groups in T2 period by comparing with that in T1 are obviously different from each other. It is worth noticing that in the central region of the Loess Plateau, the return periods of the three indicator combinations are all decreased, indicating that increased chances for occurrence of the three joint events in this region.

The mean T_and values of the other three two-dimensional combination events of {D95, I95}, {P95, I95}, and {I95, R95} are 32.65, 22.37, and 22.45, respectively. The occurrence probability is relatively small compared with the above three combinations, but the severity of the events may be far greater than the former ones. For example, {D95, I95} represents a long-time heavy rainfall event with high precipitation intensity. For these three combination events, their spatial distributions of T_and are basically the same (T1 period), and their spatial distributions of changing trend in the past 30 years are also similar to each other. The three combination events share a common feature that all of them contain the I95 index. The above results may be because the joint distributions are significantly affected by the I95 indicator, which means that the contribution of I95 to the joint return period is far greater than other indicators. The T_and of the three combined events in the past 30 years have decreased greatly in the northern and central regions, which implies the occurrence probability of extreme joint events has increased significantly. A comprehensive evaluation of T_and of the six two-dimensional indicator combinations shows that the probability of all the six joint extreme events in the central part of the region has increased.

Figure 6 shows the T_and values of three three-dimensional indicator combinations at each station in 1959–1988, and also the spatial changes in 1989–2018 compared with that in the T1 stage. The T_and values of {D95, P95, R95} are significantly smaller than those for the other two groups, with an average of only 15.53 years in 1959–1988. For the two three-dimensional events that contain the precipitation intensity indicator, the T_and values are generally greater than that of the rest one. The mean T_and of {D95, I95, R95} and {P95, I95, R95} from 1959 to 1988 is 31.2 and 24.8 years, respectively. The spatial distribution and the change direction of joint return periods of the two events are similar with each other and are also similar to those for the two-dimensional events in Figure 5b,d,e. In T1 stage, the T_and values of {D95, I95, R95} and {P95, I95, R95} in Region 1 are generally smaller than those in Region 2. Figure 6b,c shows that the change direction of T_and in T2 stage compared with T1 stage has decreased in the northern and central regions, and increased in other regions. This may also be because the T_and values of different combination events are highly affected by the precipitation intensity indicator. In the north part of Region 1, where extreme precipitation occurs frequently in T1 stage, the three-dimensional T_and get even smaller in T2 stage.

Figures S2–S5 (in the Supplementary Materials) provide the T_k and T_or values of two- and three-dimensional indicator combinations at each station in 1959–1988, and the spatial changes in 1989–2018 compared with those in the T1 stage. In two-dimensional cases at T1 stage, the T_k values are all smaller than T_and over all stations (see explanation in Section 4.1) and are similar to T_and in magnitude distribution and change direction. The average T_k values corresponding to Figure S2a–f are 1.55, 10.09, 2.37, 5.66, 1.72, and 5.73 years smaller than those of T_and, respectively. The magnitude distribution of T_or values in T1 stage are opposite to T_and, as well as the change direction in T2 stage. This is because, given the CDF values of two indicators, the value of T_and is only related to bivariate CDF like T_or (Equations (10) and (12)). With the change of one bivariate CDF, the change directions of the two return periods are opposite with each other. In three-dimensional cases at T1 stage, the average T_K values corresponding to Figure S3a–c are, respectively, 2.82, 8.80, and 6.83 years smaller than T_and. The average T_or values in Figure S5a–c are only 6.98, 5.50, and 5.82 years, respectively.

4.3. The Trend of Return Period of Multidimensional Moving Window Series

Figure 7 shows the changing trends in 10-year return period values for the four individual indicators calculated from the 30-year moving time series over all stations. The significance of changing trend and magnitude are also displayed in Figure 7.

Table 6. Statistical information on the increase and decrease trend of the four indicators in two regions: The proportion of stations with significant increasing trend (PI) and the total growth over the moving window series (TG); the proportion of stations with significant decreasing trend (PD) and total reduction over the moving window series (TD).

Indices	Region 1		Region 2
	PI/TG	PD/TD	PI/TG	PD/TD
D95	0.31/0.54	0.46/−0.63	0.16/0.53	0.42/−0.95
P95	0.19/22.04	0.38/−35.7	0.24/39.40	0.50/−37.01
I95	0.27/6.79	0.38/−6.80	0.42/12.17	0.29/−11.63
R95	0.38/0.05	0.46/−0.06	0.34/0.04	0.34/−0.05

provides statistics on the proportion of stations with significant increasing/decreasing trends, and average changes in the two regions. The results show that in the semi-arid area of the Loess Plateau (Region 1), the changing trends of D95 and P95 are similar in significance and direction. The 10-year return period values exhibit a significant decrease trend for D95 and P95 at 46% and 38% of the stations, respectively. However, some stations in the southwest area show a significant increasing trend. For these two indicators in Region 2, 42% and 50% of the stations display a significant decreasing trend, while only 16% and 24% of the stations show a significant increasing trend. Forty-two percent of the stations in Region 2 show a significant decreasing trend for I95, with a decrease ranging in magnitude from 2.2 to 31.6 mm/day. The southern region generally exhibits a significant increase in extreme precipitation intensity, which may increase the risk of flooding in that region. The distribution of stations with significant decreasing trends of R95 is consistent with that of P95. However, more stations in the southern part of Region 1 show significant increasing trends. The proportion of extreme precipitation in annual precipitation increases. This means that the precipitation concentration degree increases, which will aggravate the spatial–temporal distribution of precipitation

Figure 8 displays the significance and magnitude of the changing trends for T_and of the six two-dimensional indicator combinations based on the 30-year moving time series over all stations. Figure 8b,d,f shows that the joint return periods of three two-dimensional events {D95, I95}, {P95, I95}, and {I95, R95} in the whole region present a relatively consistent trend. The joint return periods of these three events in the entire region show significant increasing trends at 30, 25, and 29 stations, respectively. Moreover, 19 of them are shared stations, and these stations are mainly located in the south of Region 1. There are 17, 16, and 19 stations showing significant decreasing trends, 14 of which are shared stations, which are mainly located in the middle and north parts of the Loess Plateau. The average reductions in the T_and return period are 26.8, 9.25, and 9.39 years for the three events, respectively. These results suggest that the possibility of flood events with long time and high intensity in the middle of the Loess Plateau would increase. The decrease in T_and indicates that extreme precipitation events with the same degree occur more frequently. For the other three two-dimensional combined events, it can be seen that the return period of T_and for {P95, R95} is most significantly reduced in the entire region, with a total number of 24 stations showing a significant decreasing trend.

Figure 9 shows the significance and magnitude of the trends in joint return periods of T_and for three-dimensional indicator combinations over all stations. It can be seen that some stations in the central region of Loess Plateau exhibit significant downward trends for all three-dimensional events. The changing trend of T_and for {D95, P95, R95} in the middle of Region 2 shows a considerable differentiation state: The changing trend for stations in the western region decreased significantly, whereas those in the eastern region increased significantly. The stations with significant decreasing trend in T_and for {D95, P95, R95} are mainly located in the sub-humid regions, with an average reduction of 20.5 years. However, the stations with significantly decreasing return periods for {P95, I95, R95} are located near the boundary of 400 mm rainfall, with an average decrease of 13.5 years. This means that the combination event of {P95, I95, R95} in the semi-arid and semi-humid junction area tends to be more frequent.

Figures S6–S9 in the Supplementary Materials show the significance and magnitude of the trend changes for T_k and T_or for the two- and three-dimensional situations over all stations. In two-dimensional cases, the T_or and T_and basically show opposite changing trends, whereas T_k and T_and show similar changing trends and directions, with only differences in magnitude. For example, for {D95, I95}, {P95, I95} and {I95, R95}, the number of stations showing an increasing trend in T_k are 30, 23, and 25, respectively, in which 29, 22, and 24 stations are shared with the T_and cases. For the events of {D95, I95, R95} and {P95, I95, R95}, the stations with significant decreasing trends in T_k (mainly located in the central Loess Plateau) decreased by an average of 10.61 and 5.05 years, respectively. The sites with significant decrease trends in T_or are mainly located in the southwest and northeast parts of Loess Plateau.

5. Discussions and Conclusions

Compared with the previous studies, this study presents several outstanding features in multivariate risk inferences: (i) The joint risks of multidimensional indicators are evaluated in this study. (ii) Three types of joint return periods, including “AND”, “OR”, and Kendall, are evaluated for both two- and three-dimensional indicator combinations. (iii) By using the moving window approach, the proposed approach can quantitatively evaluate the spatial–temporal changes in risks of multidimensional cases. (iv) The GMM is used in fitting the marginal distributions for all indicators, which helps reduce the uncertainty brought by marginal model selection. The Kendall return period for high-dimensional copulas and two-dimensional elliptical copulas based on the random sampling algorithm is another innovation in this field.

The developed extreme precipitation assessment framework is applied in the Loess Plateau, China. The results indicate that the mean values of P95, D95, and I95 at 10-year return period in Region 1 (arid and semi-arid climate) are all smaller than those in Region 2 (semi-humid climate). However, the mean value of R95 in Region 1 is slightly higher than that in Region 2, and nearly half of the annual precipitation would occur by extreme heavy precipitation in the north of Region 1. Given the values of the indicators at 10-year return period, the combined events of {D95, P95}, {D95, R95}, {P95, R95}, and {D95, P95, R95} are extremely easy to occur in the Loess Plateau. The occurrence probability for combined events with the I95 indicator is relatively small compared with other combinations. The spatial distributions of changing trends in joint return periods for the multidimensional events containing the I95 indicator are generally the same. In the central region of Loess Plateau, the T_and values of all joint extremes are decreased. T_k values are smaller than T_and over all stations and are similar to T_and in magnitude distribution and change direction. The magnitude distribution of T_or values in T1 stage are opposite to T_and, as well as the change direction in T2 stage.

The provided assessment framework is repeatable and can be effectively transferred to other research areas to explore the spatial–temporal changes in multidimensional characteristics of extreme precipitation under climate change. Thanks to the development of climate models, this framework can also be used to project multidimensional probability variation characteristics of future extreme precipitation under different scenarios [41]. A variety of return period assessment methods are available to meet the needs of different hydrologic designs, and can help provide countermeasures for possible changes in floods and other secondary disasters.