Construction of a Time-Variant Integrated Drought Index Based on the GAMLSS Approach and Copula Function

Abstract

Construction of an integrated drought index is a fundamental task to conducting drought disaster risk management and developing drought resistance planning strategies. Given the evident non-consistent features during the drought evolution process, firstly, the GAMLSS approach was utilized to construct multiple combination scenarios of time-variant parameters and corresponding probability distribution functions. Then, the time-variant comprehensive drought index integrating the variable characteristics of precipitation and soil moisture was established by means of the copula function. Finally, the reliability of the time-variant comprehensive drought index was verified through its application in frequency analysis and return period determination of drought hazard system in Huaibei Plain, China. The application results demonstrated that: (1) The variation of the time-variant integrated drought indicator presented strong consistency with both soil moisture and precipitation during historical years in Huaibei Plain. (2) The overall variation of the drought hazard system characterized by drought duration and severity presented a gradual mitigation trend from west to east and north to south in Huaibei Plain, which agrees with the geographic differences and water resources availability distribution features. (3) Drought recognition results, including the frequency of drought events and typical drought processes with extreme grades, are in agreement with the practical statistics and observed data series. On the whole, the proposed time-variant integrated drought indicator is capable of extracting complex variation characteristics within the drought hazard evolution process, and can be further applied in drought monitoring, recognition and assessment research fields.

1. Introduction

Drought disaster, which is usually characterized by long duration, wide influencing area and severe disaster loss, is regarded as one of the typical and frequently occurring natural disasters around the world [1,2,3]. Currently, due to the dual impacts of climate change and human activities, regional climate and underlying surface conditions have changed significantly, which will undoubtedly aggravate the uncertainties of hydrological cycle processes and spatiotemporal variation of extreme hydrological events [4,5]. It has been determined by the Sixth Assessment Report of the United Nations IPCC that nearly the entire global drought evolution has presented an obvious and aggravating trend since 2020 [6,7]. In addition, China has been listed as a typical country with high drought frequencies resulting in severe natural disaster loss [6].

Drought disaster risk evaluation is a fundamental task for the implementation of drought disaster precautions and controlling strategies [8,9]. The essence of drought disaster risk evaluation is to quantitatively determine the probability distribution of drought occurrences and potential for disaster losses. These are often determined by drought frequency and return period analysis [2]. Moreover, the selection of drought indicators used to monitor and describe drought process variation is of great significance to construct a reliable drought disaster risk evaluation model. Yuan et al. summarized primary drought indicators and their application conditions in 2014 and proposed the future research direction of drought indicator analysis [10]. On the whole, in terms of single drought indicator construction, Xiao et al. (2012) determined the drought risk and return period of 42 stations in the Pearl River basin of China through standardized precipitation index (SPI) and multivariate copula function [11]. In 2015, Gao et al. explored the spatiotemporal variation characteristics of drought events through standardized precipitation index (SPI), empirical orthogonal function (EOF) and rotated empirical orthogonal function (REOF) approaches over a historical period of 54 years in Liaoning Province, China [12]. Adnan et al. (2018) comparatively discussed the differences of 15 drought indicators including SPI, standardized precipitation evapotranspiration index (SPEI), reconnaissance drought index (RDI) and improved precipitation z-index, and proposed that SPI, SPEI and RDI are most suitable to monitor drought process variation in Pakistan [13]. Rashid et al. (2019) constructed the nonstationary standardized precipitation index (NSPI) based on a generalized hydrological modeling framework to capture the nonlinear characteristics of precipitation during drought evolution processes [14]. Additionally, in terms of the development of comprehensive drought indicators, Chang et al. proposed an integrated drought index in 2017 that combines hydrological elements of precipitation, runoff and soil moisture through principal component analysis (PCA), and further discussed its reliability compared to different single-type drought indicators [15]. Hu et al. (2019) constructed a comprehensive drought index combining different meteorological, hydrological, agricultural and remote sensing drought indicators, and verified its effectiveness in drought risk comprehensive evaluation analysis [16]. The same year, Shen et al. proposed a comprehensive drought index considering precipitation, potential evapotranspiration, temperature, soil moisture and vegetation factors, and eventually conducted application research in characteristics extraction of historical drought risk distribution in Inner Mongolia, China [17].

The abovementioned single-type and integrated drought indicators were all constructed based on the consistency hypothesis of drought index variation. Actually, because of the dual impacts of climate change and human activities, the evolution of observed drought indicators presents obvious spatiotemporal non-consistency features [18,19,20]. In other words, the design of the drought resistance standard of water conservancy project based on the consistency assumption of drought indicator variation, is not reliable. It will increase the risk of drought disaster losses. In addition, comprehensive drought indicators are usually obtained through weight combination, machine learning algorithm and multivariate combined distribution estimation approaches [21]. The weight combination method primarily relies on the estimated weight and linear correlation hypothesis of different single drought indicators. Machine learning algorithms are incapable of revealing the influencing mechanics of integrated drought indicators to actual drought evolution. Copula functions are frequently applied to develop a multi-variable combined distribution function, which has been widely applied in drought risk frequency analysis. Therefore, in this study, we applied the generalized additive models for location, scale and shape (GAMLSS) method and copula function to construct a meteorological (precipitation) and agricultural (soil moisture) comprehensive drought index (CTVDI) considering the time-variant characteristics of each, then the CTVDI index was utilized to conduct drought frequency analysis and return period calculation in Huaibei Plain, China to test its reliability and effectiveness. The main content of the manuscript is organized as follows: In Section 1, the introduction of the study area and framework of the manuscript were explicitly presented. In Section 2, based on the brief introduction of methodologies utilized in this study including GAMLSS and copula functions, the construction procedures of the integrated drought indicator CTVDI are elaborately presented. In Section 3, the proposed integrated drought indicator CTVDI is further applied in drought disaster risk analysis to verify its reliability, and the drought event occurring probability corresponding to different grades of drought duration and severity as well as return period of typical historical drought processes in Huaibei Plain, China are determined. The primary research findings will be beneficial and favorable for implementation of regional drought resistance schemes and strategies.

2. Study Area and Research Framework

2.1. Study Area and Data Sources

Huaibei Plain, China was utilized as the case study area in this manuscript. Huaibei Plain, located in the northern Anhui Province and central Huaihe river basin, is bordered by Shandong Province in the north, Jiangsu Province in the east and Henan Province in the west. It has jurisdiction over six districts and 27 counties including Bozhou, Fuyang, Suzhou, Huaibei, Huainan and Bengbu cities. The total area of Huaibei Plain is 37,400 km², which includes cultivated land area of 21,400 km² accounting for more than 50% to the total area of Anhui Province [22]. For the impacts of subhumid monsoon climate, the annual average precipitation of Huaibei Plain is approximately 869.6 mm but differs significantly within and between different years [2,22]. The annual precipitation of a wet year is about three to four times compared to dry years, and the precipitation during flood season (from June to September) accounts for more than 75% of the total [2]. In addition, the annual average runoff of Huaibei Plain is about 7772 million m³, while the annual variation difference is more obvious compared to precipitation [15,22]. For a long time, because of inappropriate land resources development, a sharp decrease in vegetation and excessive deforestation, etc., water and soil resources erosion became serious. The groundwater levels dropped prominently as well in Huaibei Plain, which resulted in frequent flooding and drought disasters. Drought disasters have occurred in 48 of the past 60 years, and drought duration, severity and losses all presented an aggravating trend especially after the 1990s. For instance, the drought disaster situation of Anhui Province in 2019 was the most serious during the past 40 years, and drought variation in nearly 51 counties reached severe grade. Therefore, a great variety of statistics and existing research findings indicate that it is urgent to conduct integrated drought index construction and practical application research in Huaibei Plain, Anhui Province.

The primary data series utilized in this study includes historical monthly precipitation and soil moisture data series as well as drought disaster loss statistics of different cities in Huaibei Plain from 1960 to 2014. We accessed the historical drought disaster loss statistics through Statistical Yearbook of Anhui Province, EPS statistic data platform of China (https://www.epsnet.com.cn/index.html#/Index) and China’s economic and social big data research platform (https://data.cnki.net/) on 10 November 2022. The historical monthly hydrological data series can be obtained from China Integrated Meteorological Information Sharing System (CIMISS, http://data.cma.cn, accessed on 10 November 2022) and land hydrological surface modeling model (Noah-MP). In addition, the original precipitation and soil moisture data series is grid data with spatial resolution of 0.25°, which was processed to obtain historical monthly hydrological data series of different cities in Huaibei Plain through spatial kriging interpolation calculation [22].

2.2. Research Framework

The explicit framework and relationship structure of different calculation modules used in this manuscript are clearly illustrated in Figure 1.

It can be demonstrated from Figure 1 that: (1) The primary task of Methodologies (top half of the figure) is to introduce the fundamental principles of the approaches applied in this study, including GAMLSS model and copula function, and then propose the derivation procedures of meteorological (precipitation) and agricultural (soil moisture) comprehensive time-variant drought indicator (CTVDI). Specifically, the key steps to establish the CTVDI comprise the derivation of time-variant probability density function (PDF) of meteorological and agricultural drought indicators through the GAMLSS model, the estimation of combined distribution function of the CTVDI through copula function, as well as the determination of CTVDI data series through inverse operation of combined distribution function. (2) In terms of Application (bottom half of the figure), the primary intention is to test the reliability of the CTVDI through its application in Huaibei Plain, China. Specifically, the main tasks of this section includes estimation error analysis of time-variant PDFs of precipitation and soil moisture using the GAMLSS model, drought process recognition analysis based on CVTDI through run theory and drought characteristic variables (duration and severity) extraction, multivariable combined frequency and return period determination through copula function, as well as the reliability verification of the proposed CVTDI through comparative analysis with the actual observed statistics. In addition, the application requirements, suggestions, limitations of the proposed index and future work plans are discussed as well.

3. Methodologies

3.1. Generalized Additive Models for Location, Scale and Shape (GAMLSS) Method

3.1.1. Introduction to GAMLSS

Generalized additive models for location, scale and shape (GAMLSS), proposed by Prof. Rigby and Stasinopoulos, is a semiparametric regression model. GAMLSS is frequently applied to derive the functional relationship between the parameters of probability distribution function of research object (namely response variable) and corresponding explanatory variables under nonstationary circumstances [23,24,25]. Assuming the PDF of observed data series y_t (t = 1, 2, …, n) of stochastic variable Y is denoted as f(y_t|θ_t), and its distribution parameter series varying with time parameter t is denoted as θ_t = (θ_t1, θ_t2, θ_t3, θ_t4) = (μ_t, σ_t, ν_t, τ_t), in which, variables μ and σ are defined as position parameter and scale parameter, corresponding to average vector and mean square deviation vector of stochastic variable Y, respectively. Variables ν and τ are uniformly defined as shape parameter, corresponding to the skewness vector and kurtosis vector of stochastic variable Y separately [25,26]. If denoting g_m(∙) as the monotonic function between parameter θ_m and its corresponding explanatory variables X_m as well as random effect items, then it becomes

$$\begin{array}{c}{g}_m\left({\theta }_m\right)=X_m·{\beta }_m+{\displaystyle \sum _{j=1}^{J_m}}{D}_{jm}·{\gamma }_{jm}\end{array}$$

(1)

where, parameter m = 1, 2, 3 and 4, (X_m)_n×Jm is the explanatory variables matrix, β_m is the regression coefficient vector with the length equaling J_m, J_m denotes the number of random effect variables of the m_th parameter, (D_jm)_n×qjm is the random effect variable matrix, γ_jm is normally distributed random variable vector with the length equaling q_jm, and q_jm is the number of random factors corresponding to j_th random effect variable [26]. Equation (1) is primarily utilized to describe the linear or nonlinear relationship of explanatory variables as well as linear relationship of random effect items for different distribution parameters, if neglecting the influences of random effect items, the original monotonic function of GAMLSS model can be simplified as follows

$$\begin{array}{c}{g}_m\left({\theta }_m\right)=X_m·{\beta }_m\end{array}$$

(2)

If assuming the explanatory variable of GAMLSS model satisfies a two-parametric probability distribution pattern, then the estimated value of two parameters will be obtained through RS algorithm, i.e., taking the maximum likelihood function of regression coefficient as the objective function [26,27]. Therefore, we primarily discuss the linear function relationship of distribution parameters vector μ_t and σ_t varying with time parameter t in this study, and the distribution parameter estimation method of GAMLSS model can be further simplified as a linear function mode varying with time parameter t [28,29], as follows

$$\begin{array}{c}\{\begin{array}{c}\mu \left(t\right)=a_0+a_{1}t\\ \sigma \left(t\right)=b_0+b_{1}t\end{array}\end{array}$$

(3)

where, parameters a₀ and a₁ are combined coefficients of position parameter μ_t, parameters b₀ and b₁ are combined coefficients of scale parameter σ_t, and parameter t is time concomitant variable.

3.1.2. Parameter Optimization of GAMLSS

Based on the analysis above, it can be concluded that the primary task of GAMLSS model is to derive the monotonic function g_m(∙) between explanatory variable X and its distribution parameters including time-variant position variable μ_t and scale variable σ_t [28]. Actually, the entire function derivation and simulation analysis of GAMLSS method is frequently accomplished by the program package of R language software, which can provide a variety of distribution functions for users. Combining previous application study findings of GAMLSS model in related meteorological and hydrological fields, we selected five types of two-parametric distribution functions including standard normal distribution (NO), log normal distribution (LOGNO), Weibull distribution (WEI), gamma distribution (GA) and Gumbel distribution (GU), to derive the time-variant PDF of drought indicators precipitation and soil moisture of different cities in Huaibei Plain [26,28], which is shown in

Table 1. Five types of two-parametric distribution functions of the GAMLSS model.

Name	PDF	Determination of Distribution Parameters
NO	$f_X\left(x|\mu ,\sigma \right)=\frac{1}{\sqrt{2\pi \sigma }}\cdot \exp \left[-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}\right]$	$E\left(X\right)=\mu$ $Var\left(X\right)={\sigma }^2$
LOGNO	$f_X\left(x|\mu ,\sigma \right)=\frac{1}{\sqrt{2\pi \sigma }}\cdot \frac{1}{x}\cdot \exp \left[-\frac{\log \left(x\right)-\mu )}{2{\sigma }^2}\right]$	$E\left(X\right)=\exp {\left({\sigma }^2\right)}^{1/2}\cdot e^{\mu }$ $Var\left(X\right)=\exp \left({\sigma }^2\right)\cdot \left[\exp \left({\sigma }^2\right)-1\right]$
WEI	$f_X\left(x|\mu ,\sigma \right)=\frac{\sigma x^{\sigma -1}}{{\mu }^{\sigma }}\cdot \exp \left[-{\left(\frac{x}{\mu }\right)}^2\right]$	$E\left(X\right)=\mu ·G\left(1/\sigma \right)$ $Var\left(X\right)={\mu }^2·\left\{G\left(\frac{2}{\sigma +1}\right)-{\left[G\left(\frac{1}{\sigma +1}\right)\right]}^2\right\}$
GA	$f_X\left(x|\mu ,\sigma \right)=\frac{1}{{\left({\sigma }^2\mu \right)}^{1/{\sigma }^2}}\cdot \frac{x^{1/\left({\sigma }^2-1\right)}\cdot e^{-x/\left({\sigma }^2\mu \right)}}{G\left(1/{\sigma }^2\right)}$	$E\left(X\right)=\mu$ $Var\left(X\right)={\mu }^2·{\sigma }^2$
GU	$f_X\left(x|\mu ,\sigma \right)=\frac{1}{\sigma }\cdot \exp \left[\left(\frac{x-\mu }{\sigma }\right)-\exp \left(\frac{x-\mu }{\sigma }\right)\right]$	$E\left(X\right)=\mu -\mu \sigma \approx \mu -0.57722\sigma$ $Var\left(X\right)={\pi }^2·{\sigma }^2/6$

Note: G(·) denotes the gamma function.

.

To fully reveal the nonlinear features of distribution parameters varying with the time variable, the distribution functions provided by the time-variant GAMLSS program package of R language software were divided into three categories in this study: Scenario 1, position variable μ_t and scale variable σ_t were all constant, revealing that the variation of distribution parameters μ_t and σ_t all satisfied stationary features. Scenario 2, position variable μ_t varied with time variable t while scale variable σ_t was constant. Scenario 3, position variable μ_t and scale variable σ_t were all varying with time variable t [26,29]. Furthermore, the fitting performance of different distribution functions was evaluated by a generalized Akaike information criterion (AIC) [26,27], and the distribution function with minimum AIC value was recognized as the optimal fitting function of distribution parameters. In other words, if the AIC value of distribution functions of Scenario 1 was minimum, the variation of distribution parameters μ_t and σ_t did not present nonstationary characteristics.

3.2. Copula Function

Copula function, firstly proposed by Sklar in 1959, has been widely applied in stochastic simulation and statistical analysis fields [30,31]. The primary purpose of copula function is to derive the multivariate joint distribution function through their univariate probability distribution functions. The widely applied univariate distribution pattern in hydrological frequency analysis comprises exponential distribution, normal distribution, gamma distribution, etc. In terms of the selection of copula joint function to derive the multivariate combined distribution function, the binary Archimedean copula function, for the advantages of rigorous logic structure, simple calculation formula and less estimated parameter, has been widely applied [31,32,33]. Specifically, binary Archimedean copula function primarily includes three types of joint functions, i.e., GH copula, Clayton copula and Frank copula, and the ultimate purpose of the application of copula function is to derive the joint probability distribution function of CTVDI based on the determination of PDFs of drought indicators including precipitation and soil moisture through GAMLSS model. Supposing U = u(x) and V = v(y) represent the derived PDFs of univariate drought indicators precipitation and soil moisture separately, and F(x, y) = c(u, v) denotes the joint probability distribution function of CTVDI, which can be denoted as follows [30,31]

$$\begin{array}{c}c\left(u,v\right)=exp\left\{-{\left[{\left(-\ln \left(u\right)\right)}^{\theta }+{\left(-\ln \left(v\right)\right)}^{\theta }\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\theta $}\right.}\right\},\tau =1-\frac{1}{\theta },\theta \in \left[1,\infty \right)\end{array}$$

(4)

$$\begin{array}{c}c\left(u,v\right)={\left[{\left(u\right)}^{-\theta }+{\left(v\right)}^{-\theta }-1\right]}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\theta $}\right.},\tau =\frac{2}{2+\theta },\theta \in \left[0,\infty \right)\end{array}$$

(5)

$$\begin{array}{c}c\left(u,v\right)=\frac{1}{\theta }ln\left[1+\frac{\left(e^{-\theta u-1}\right)·\left(e^{-\theta v-1}\right)}{e^{-\theta }-1}\right],\tau =1+\frac{4}{\theta }\left(\frac{1}{\theta }{{\displaystyle \int }}_0^{\infty }\frac{t}{e^t-1}dt-1\right),\theta \in R\end{array}$$

(6)

where, θ is an unknown parameter to describe the correlationship of precipitation and soil moisture and can be derived through its relationship with Kendall correlation coefficient τ [30,32,33]. Afterwards, the final CTVDI could be obtained through the inverse operation of joint distribution function F(x, y) [34].

3.3. Calculation Procedures of Integrated Index CTVDI

On the basis of the abovementioned discussion, it was obvious that the development of the CTVDI and its reliability verification analysis could be accomplished through the following steps:

Step 1: Determination of time-variant PDF of univariate drought indicators precipitation and soil moisture through GAMLSS model. We constructed three scenarios to explore the stationary or nonstationary relationship features of univariate drought indicators according to their different combination of position variable μ_t and scale variable σ_t varying with time variable, and then the optimal time-variant and monthly scale PDF. The u(x) and v(y) of precipitation and soil moisture were eventually derived under AIC principle for different cities in Huaibei Plain in this study.

Step 2: Determination of combined distribution function of integrated drought indicator CTVDI through copula function. As indicated in Section 3.2, three types of Archimedean copula functions, including GH, Clayton and Frank were utilized to derive the time-variant combined probability distribution function c(u, v) of precipitation and soil moisture in this manuscript. In addition, the function fitting performance of different cities was also verified through different error parameters.

Step 3: Derivation of CTVDI series. The CTVDI was determined through inverse operation of its joint probability distribution function c(u, v), as follows

$$\begin{array}{c}CTVDI={\phi }^{-1}\left(c\right)\end{array}$$

(7)

where, φ presents function inverse operation of standard normal distribution, which was realized through R language software in this study [34]. Furthermore, the derived CTVDI series was applied to conduct drought processes recognition during historical years in different cities of Huaibei Plain through run theory [35]. Run theory is a fundamental method to recognize drought event processes through single or coupled drought indicators by means of initial drought identification and drought merging analysis [31], in which, the determination of threshold values of R₁, R₂ and R₃ are crucial to accurately derive the starting and end time of drought process [31,35]. In this study, according to Standard for Hydrological Information and Hydrological Forecasting (GB/T22482-2008), the threshold values of parameters R₁, R₂ and R₃ are 0, −0.5 and −1, respectively, and then the historical drought event samples could be obtained based on the proposed CTVDI. Finally, drought characteristic variables, including drought duration and severity corresponding to different drought events could be obtained to describe drought variation conditions [35].

Step 4: Drought frequency analysis based on the CTVDI series and its reliability verification. The reliability verification of the proposed CTVDI was primarily accomplished by means of its application in future drought frequency analysis and return period calculation of typical historical drought events. Generally, drought frequency is defined as the occurrence probability of drought events with drought characteristic variables of duration and severity exceeding certain levels in a certain future period [33,35]. Suppose the derived probability distribution function of drought duration D and drought severity S were denoted as F_D(d) and F_S(s), respectively, then the drought frequency, i.e., the occurrence probability of future drought events with drought characteristic variables satisfying D > d and S > s simultaneously, could be represented as follows [19,35]

$$\begin{array}{c}P\left(D>d,S>s\right)=1-F_D\left(d\right)-F_S\left(s\right)+F_{D,S}\left(d,s\right)\end{array}$$

(8)

where, F_{D, S}(d, s) is the combined probability distribution function of drought duration and severity derived by copula function as introduced in Section 3.2. In addition, the return period of extreme hydrological events is another important concept in the field of hydraulic engineering designing processes. The return period of drought events with drought characteristic variables satisfying D > d and S > s, simultaneously, could be denoted as follows [19,36]

$$\begin{array}{c}P\left(D>d,S>s\right)=\frac{E\left(L\right)}{P\left(D>d,S>s\right)}=\frac{E\left(L\right)}{1-F_D\left(d\right)-F_S\left(s\right)+F_{D,S}\left(d,s\right)}\end{array}$$

(9)

where, E(L) denote the expectation of drought intervals, which could be represented as the sum of average of drought and non-drought duration.

4. Results and Discussion

To further validate the feasibility for the application of GAMLSS and copula function approaches in the field of time-variant evolution characteristics exploration of drought systems, we utilized the proposed CTVDI in Huaibei Plain, Anhui Province, China to conduct future drought frequency analysis and historical drought return period determination in this manuscript. Meanwhile, given the length constraint of the manuscript, taking August as an example, the establishment of the monthly CTVDI is emphatically discussed in this study, which is the same for other months.

4.1. Performance Analysis of GAMLSS Model

The derivation of time-variant PDF of univariate drought indicator precipitation (P) and soil moisture (SM) by means of GAMLSS model is the fundamental task for the development of integrated CTVDI. According to the AIC principle of time-variant PDF fitting through GAMLSS model, the variation of global fitting goodness of P and SM for different distribution patterns under different combination scenarios of position and scale parameters for different cities in Huaibei Plain, is indicated in

Table 2. Variation of fitting goodness under different distribution and scenarios in Huaibei Plain.

Name	Type	Huaibei		Suzhou		Bozhou		Bengbu		Fuyang		Huainan
		P	SM	P	SM	P	SM	P	SM	P	SM	P	SM
Scenario 1: constant μ and σ	NO	601	266	593	276	602	252	600	259	602	261	592	275
	LOGNO	602	268	598	279	596	255	603	264	589	266	604	280
	WEI	598	262	591	273	596	247	596	247	595	251	589	259
	GA	597	267	592	278	594	254	597	262	590	265	592	278
	GU	614	263	606	274	625	248	613	246	625	251	604	257
Scenario 2: μ = f(t) and constant σ	NO	602	257	595	263	603	246	601	259	603	256	593	271
	LOGNO	603	259	599	266	596	249	604	264	589	261	603	276
	WEI	598	253	593	259	598	241	598	243	596	246	591	259
	GA	599	258	594	265	595	248	599	262	591	259	593	275
	GU	615	254	607	258	627	240	615	241	626	245	606	256
Scenario 3: μ = f(t) and σ = s(t)	NO	604	259	597	264	604	248	601	257	605	257	592	271
	LOGNO	602	261	595	268	597	251	597	262	590	262	601	277
	WEI	599	254	594	259	598	242	596	245	597	247	586	256
	GA	599	260	593	267	595	250	595	261	592	260	585	275
	GU	617	254	609	260	624	242	613	243	627	246	603	254

Note: the bold and red data represent the minimum goodness of fit.

, then the corresponding optimal distribution pattern with minimum fitting goodness for precipitation and soil moisture could be recognized.

It can be seen from Table 2 that: (1) The optimal fitting distribution pattern of historical precipitation and soil moisture data series in August differs significantly for different cities in Huaibei Plain. The variation of historical precipitation series in Huaibei, Suzhou, Bozhou and Fuyang, belonging to Scenario 1 with the minimum PDF goodness of fit compared to other combination scenarios, presents obvious stationary characteristics. The variation of historical soil moisture series except for Huainan city, belonging to Scenario 2, presents distinct non-stationary characteristics for position parameter μ. The variation of both historical precipitation and soil moisture series in Huainan city, belonging to Scenario 3, present remarkable non-stationary features for both position parameter μ and scale parameter σ. (2) Compared with previous research findings satisfying the stationary assumptions for the variation of drought indicators, the derived time-variant PDF of precipitation and soil moisture in different cities of Huaibei Plain is expected to be more effective to describe its actual evolution processes.

To further verify whether univariate drought indicators precipitation and soil moisture satisfy the optimal distribution pattern, the residual distribution corresponding to different optimal distribution patterns are discussed as well by means of average (AV), variance (VA), skewness coefficient (SC), kurtosis coefficient (KC) and Filliben coefficient (FC) in this study, and if the variation of AV approaches 0, VA approaches 1, SC approaches 0, KC approaches 3 and FC approaches 1, it could be concluded that the residual distribution satisfies a standard normal distribution pattern and the corresponding optimal distribution pattern is also more reasonable [26,29]. The statistical results of AV, VA, SC, KC and FC corresponding to the optimal fitting distribution pattern and time-variant parameters of univariate drought indicators precipitation and soil moisture are indicated in

Table 3. Variation of residual distribution parameters of optimal distribution pattern in Huaibei Plain.

City	Name	Type	Distribution Parameter	Residual Distribution Parameter
				AV	VA	SC	KC	FC
Huaibei	P	GA	μ = exp(4.7862)	0.0018	1.0192	−0.2744	2.5630	0.9929
			σ = exp(−0.7287)
	SM	WEI	μ = exp(3.6082 − 0.0018·t)	−0.0019	1.0158	0.0465	2.4647	0.9922
			σ = exp(2.8581)
Suzhou	P	WEI	μ = exp(4.9376)	0.0003	1.0013	0.1000	2.5842	0.9938
			σ = exp(0.9520)
	SM	GU	μ = (36.8974 − 0.0835·t)	0.0005	0.9956	0.1559	2.4970	0.9908
			σ = exp(0.7262)
Bozhou	P	GA	μ = exp(4.7933)	0.0000	1.0185	−0.0009	2.8585	0.9966
			σ = exp(−0.7785)
	SM	GU	μ = (34.7393 − 0.0484·t)	0.0007	0.9987	0.1415	2.6827	0.9953
			σ = exp(−0.5639)
Bengbu	P	GA	μ = (4.6641 − 0.0053·t)	0.0054	1.0252	−0.0231	2.0978	0.9879
			σ = exp(−0.4236 − 0.0146·t)
	SM	GU	μ = (38.3270 − 0.0394·t)	−0.0063	1.0873	−0.3833	2.7571	0.9889
			σ = exp(0.5237)
Fuyang	P	LOGNO	μ = 4.6452	0.0000	1.0185	−0.1165	2.7462	0.9946
			σ = exp(− 0.7432)
	SM	GU	μ = 36.1412 − 0.0461·t	−0.0005	1.0271	−0.0491	2.8125	0.9968
			σ = exp(0.5921)
Huainan	P	GA	μ = exp(4.5209 + 0.0065·t)	0.0061	1.0253	−0.1981	2.3835	0.9901
			σ = exp(−0.3317 − 0.0162·t)
	SM	GU	μ = 37.1387 − 0.0348·t	−0.0109	1.1009	−0.4277	2.5492	0.9872
			σ = exp(0.2777 + 0.0122·t )

Note: AV, VA, SC, KC and FC represent average, variance, skewness coefficient, kurtosis coefficient and Filliben coefficient, respectively.

.

Meanwhile, the normal QQ and normal worm figures were also applied to verify the reliability of time-variant PDF derivation of precipitation and soil moisture in Huaibei Plain. The normal QQ figure actually represents quantile scatterplot with the estimated value of standard normal distribution and practical value of residual distribution as X and Y axis, respectively, and the normal worm figure is utilized to describe the deviation of practical residual distribution from critical line of standard normal distribution with confidence coefficient equaling 95% [25,26,29]. The practical residual distribution of optimal fitting distribution could be regarded as satisfying normal distribution pattern if the sample scatterplot of residual distribution is approximately linear or not exceeding the critical line, and the residual distribution parameters of AV, VA, SC, KC and FC could also be roughly estimated from normal QQ and worm figures. The normal QQ and worm figures representing the residual distribution of optimal fitting distribution patterns of precipitation and soil moisture in different cities of Huaibei Plain are indicated in Figure 2 and Figure 3, respectively.

It could be concluded from Table 3 and Figure 2 and Figure 3 that: (1) As for the residual distribution of optimal fitting distribution patterns of precipitation and soil moisture in August for different cities in Huaibei Plain, the variation of average (AV), variance (VA), skewness coefficient (SC) and kurtosis coefficient (KC) was approaching 0, 1, 0 and 3, respectively, gradually, which revealed that the residual distribution satisfies standard normal distribution assumption. In addition, the Filliben coefficient (FC) was also approaching 1, which further verified the independence assumption of residual distribution. (2) The residual samples of both precipitation and soil moisture were all uniformly distributed nearing the 45° line according to the normal QQ plot of residual sample, and all between the U-shaped and inverted U-shaped curves as well (i.e., threshold line of confidence coefficient equaling 95%) according to the normal worm plot of residual sample, indicating that the corresponding optimal fitting distribution patterns of precipitation and soil moisture all had good fitting performance with historical observed samples. (3) In the same way, the optimal time-variant fitting patterns of PDF of precipitation and soil moisture in other months for different cities were also derived using the GAMLSS model, which could provide better model foundation for the development of time-variant integrated drought indicator.

4.2. Derivation of CTVDI Series and Its Application in Drought Process Recognition

As presented in Section 3.3, based on the derivation of optimal time-variant PDF of monthly precipitation and soil moisture series of different cities in Huaibei Plain through GAMLSS model, the time-variant joint probability distribution function of precipitation and soil moisture was obtained utilizing three typical copula functions, and the corresponding time-variant monthly comprehensive drought indicator CTVDI could also be derived through function inverse operation according to Equation (7). The monthly variation of CTVDI from 1960 to 2014 of different cities in Huaibei Plain are shown in Figure 4. Moreover, the linear and Kendall correlation coefficients (denoted as LCC and KCC) were utilized as well to reveal the correlation features of proposed CTVDI with precipitation and soil moisture, as indicated in

Table 4. Correlation coefficient between different drought indicators in Huaibei Plain.

Name	Type	Huaibei	Suzhou	Bozhou	Bengbu	Fuyang	Huainan	Average
P and CTVDI	LCC	0.68	0.72	0.73	0.77	0.75	0.74	0.73
	KCC	0.68	0.67	0.71	0.79	0.68	0.73	0.71
SM and CTVDI	LCC	0.83	0.82	0.92	0.91	0.88	0.85	0.87
	KCC	0.73	0.71	0.87	0.82	0.76	0.77	0.78

Note: LCC and KCC denote linear and Kendall correlation coefficient, respectively.

.

It could be concluded from Figure 4 and Table 4 that: (1) The CTVDI had strong correlation with the variation of soil moisture series in comparison with precipitation, the average of LCC and KCC between historical CTVDI and soil moisture series was 0.87 and 0.78, respectively, for different cities in Huaibei Plain. (2) The average of LCC and KCC between historical CTVDI and precipitation series was 0.73 and 0.71, separately, in Huaibei Plain, belonging to strong correlation grade as well. Therefore, it is suggested that the variation of derived time-variant integrated drought indicator CTVDI is reasonable from the perspective of correlation analysis and is suitable for its application in drought process recognition and drought frequency analysis.

Historical drought process recognition is the foremost task of drought frequency analysis by means of the CTVDI, and the historical drought events recognition based on CTVDI index from 1960 to 2014 in Huaibei Plain was accomplished through run theory in this study, and the primary drought characteristic variables were extracted as well, in which, drought duration is defined as the accumulated length of time from the beginning to the end of certain drought event, and drought severity is defined as the accumulated value of drought indicator (i.e., CTVDI in this study) throughout the entire drought process. The statistical result of historical drought characteristic variable, 1960–2014 in Huaibei Plain, is given in

Table 5. Statistic result of historical drought characteristic variable, 1960–2014 in Huaibei Plain.

Name	Huaibei	Suzhou	Bozhou	Bengbu	Fuyang	Huainan
Drought event amount	91	82	80	80	74	76
Average of drought duration (month)	3.1	3.04	3.19	2.53	3.23	2.88
Average of non-drought duration (month)	4.13	4.37	4.93	5.86	5.04	6.13
Maximum of drought duration (month)	16	12	14	13	14	9
Average of drought severity	4.22	4.05	4.26	3.57	4.4	3.98
Maximum of drought severity	20.89	15.57	18.61	18.95	18.43	15.94
Kendall correlation coefficient of drought duration and severity	0.78	0.72	0.77	0.7	0.73	0.74

.

It can be concluded from Table 5 that: (1) The overall variation of drought duration and severity of historical drought events basically presented a consistent trend in Huaibei Plain. Further, drought hazard situation mitigated gradually from west to east and north to south during historical years, which agrees with the geographic differences and water resources availability distribution. (2) Comparatively, the variation of drought hazard in Huaibei city was the most severe during historical years with the maximum drought duration, severity and frequency as well. In addition, the variation of drought characteristic variables of typical drought events in Huaibei city was consistent with the actual statistic result [37]. Therefore, the reliability of the CTVDI was verified again from the perspective of comparative analysis between the variation of drought characteristic variables of typical drought events and actual historical statistic data, and afterwards, the future occurrence probability of drought events with different grades of drought duration and severity was estimated in this study.

4.3. Drought Frequency and Return Period Determination Analysis

Drought duration and severity are two significant characteristic variables to describe the variation of drought hazard evolution processes, and a great deal of previous research has suggested fitting the variation of drought duration and severity through Pearson-III, exponential and gamma distribution patterns [19,35]. Therefore, to further discuss the joint probability of future drought events exceeding certain grades of duration and severity, based on the obtain of samples of duration and drought of historical drought events according to empirical frequency formulas, the fitting PDFs of duration and severity in different cities of Huaibei Plain were determined on the assumptions of Pearson-III, exponential and gamma distribution patterns, as indicated in Figure 5A,B. Ultimately, based on the matching situation of different frequency curves with the actual sample scatters of drought duration and severity, the joint probability distribution of different combination scenarios these variables for different cities were determined by means of GH copula function, which has better fitting performance for describing the correlation of upper sample scatters [35], as indicated in Figure 5C, and the occurrence probability of future drought processes exceeding certain grades of duration and severity could be also determined according to Equation (8).

It could be determined from Figure 5 that the exponential and Pearson-III distribution were preferable to fit the distribution features of drought duration and severity of historical drought events in Huaibei Plain, with the minimum root mean square error (RMSE) of 0.8934 and 0.8232 compared to Pearson-III and gamma distribution patterns, respectively. Accordingly, the fitting PDF of drought duration and severity was more reliable to derive the combined probability distribution function through the GH copula function. In addition, based on the division of drought hazard grades by means of drought process return period (T) including light drought (T ≤ 2a), moderate drought (2a < T ≤ 6a), severe drought (6a < T ≤ 20a) and extreme drought (T > 20a) [38], the variation for the occurrence of historical drought events with different severity grades in Huaibei Plain was calculated, as shown in

Table 6. Frequency of historical drought events, 1960–2014 in Huaibei Plain.

City	T ≤ 2a (Light Drought)	2a < T ≤ 6a (Moderate Drought)	6a < T ≤ 20a (Severe Drought)	T > 20a (Extreme Drought)	Total Amount
Huaibei	60	18	8	5	91
Suzhou	52	18	7	5	82
Bozhou	51	16	9	4	80
Bengbu	50	17	9	4	80
Fuyang	46	16	11	6	79
Huainan	45	15	8	3	71

. In particular, the occurrence of extreme drought events was usually closely emphasized because of its serious drought disaster loss and socioeconomic impacts by the government to establish effective drought resistance schemes. Therefore, based on the estimation of joint frequency of drought duration and severity, the joint frequency, return period as well as primary drought characteristic variable of extreme drought events during historical years in different cities of Huaibei Plain were also extracted according to Equations (8) and (9), which is indicated in

Table 7. Statistic results of extreme drought processes, 1960 to 2016 in Huaibei Plain.

City	No.	Time of Drought Event	Drought Duration /Month	Drought Severity	Joint Frequency	Return Period /Year
Huaibei	1	August 1966–June 1967	11	17.56	0.0094	64.2
	2	November 1967–July 1968	9	14.19	0.0239	25.2
	3	October 1991–January 1993	16	14.08	0.0038	160.1
	4	September 1998–September 1999	12	13.23	0.0152	39.6
	5	March 2001–November 2001	9	18.43	0.0074	81.2
Suzhou	1	May 1966–January 1967	9	13.78	0.0183	34.2
	2	September 1967–July 1968	11	15.57	0.0088	70.3
	3	January 1978–December 1978	12	10.46	0.0074	87.9
	4	March 1981–April 1982	12	13.54	0.0215	28.7
	5	March 2001–November 2001	9	14.95	0.0133	47.5
Bozhou	1	May 1966–June 1967	14	18.62	0.0056	121.6
	2	March 1978–January 1979	11	12.17	0.0202	33.5
	3	October 1991–August 1992	11	11.51	0.0203	33.3
	4	March 2001–November 2001	9	18.26	0.0084	80.3
Bengbu	1	May 1966–June 1967	13	18.95	0.0044	157.9
	2	November 1967–June 1968	8	13.62	0.0206	34.3
	3	March 1978–March 1979	13	15.55	0.0064	109.9
	4	October 2010–June 2011	9	11.25	0.0253	27.6
Fuyang	1	December 1960–July 1961	8	11.05	0.0275	27.3
	2	October 1967–May 1968	8	15.94	0.0072	104.3
	3	August 1978–March 1979	8	7.62	0.0304	24.7
	4	November 1983–May 1984	7	11.11	0.0363	20.7
	5	October 1991–June 1992	9	11.04	0.0168	44.7
	6	March 2001–November 2001	9	12.98	0.0153	50.4
Huainan	1	May 1966–June 1967	14	20.89	0.0051	136.1
	2	February 1978–March 1979	13	15.6	0.0087	79.1
	3	March 2001–November 2001	9	20.49	0.0077	89.1

.

It was suggested from Figure 5 and Table 6 and Table 7 that: (1) From the perspective of the frequency of historical drought events, the drought hazard situation mitigated gradually especially for light and moderate drought grades from north to south, and the total frequency of drought events decreased from the maximum of 91 in Huaibei city to the minimum of 71 in Huainan city, which was consistent with the variation trend of drought duration and severity as shown in Table 5. (2) The difference for the frequency of historical drought events with extreme grade in Huaibei Plain was not evident. Moreover, the typical drought process recognition result of extreme grade through time-variant drought indicator CTVDI, i.e., 1966, 1978 and 2001, was basically consistent with the historical statistics according to Wen et al., 2007 [37]. For instance, according to historical statistics, the drought hazard variation in summer, autumn and winter in 1966 was tremendously serious in Huaibei Plain. The monthly precipitation in August, September and October decreased approximately 70% compared to the annual average level, and drought influencing area of Huaibei Plain in 1966 accounted for nearly 60% of the entire Anhui Province. Similarly, the variation in drought duration, drought-affected area and drought disaster loss of Huaibei Plain in 1978 was the most serious since the founding of China, and grain production reduced 714.24 million kg by drought in Huaibei Plain, accounting for 52% of the total grain reduction of Anhui Province.

In conclusion, the drought hazard variation characteristics of Huaibei Plain during the historical years derived through the time-variant integrated drought indicator are in good agreement with the practical statistics and observed data series. This revealed that the integrated drought indicator can better extract the features of drought evolution processes than single drought indicators and is suitable for application in the drought disaster risk management field. In addition, the variation process for the occurrence, development and mitigation of drought events is extremely complicated. Incorporation of hydraulic engineering regulation and groundwater resources utilization for the establishment of comprehensive drought indicator and drought disaster risk management in Huaibei Plain is indispensable. This represents a limitation of the study as well as a research direction for future work.

5. Conclusions

In this study, based on the non-consistent characteristics of drought indicators and deficiencies of drought process identification through single drought indicators, the time-variant integrated drought indicator, coupling precipitation and soil moisture characteristics was proposed by means of GAMLSS and copula functions, and the reliability and effectiveness of the proposed comprehensive drought indicator were verified through its application in drought hazard frequency and drought return period determination in Huaibei Plain, Anhui Province. Together, the primary findings of this study can be summarized as follows:

(1) Five types of probability density function and four patterns of combination scenarios of time-variant parameters were utilized to derive the time-variant PDF of precipitation and soil moisture by GAMLSS method. The optimal PDF of different single-drought indicators was determined based on the comparative analysis of fitting error, normal QQ and worm figures of fitting residual items. The variation of proposed time-variant comprehensive drought indicator presented strong correlation features with single precipitation and soil moisture indicator.

(2) According to the total frequency of drought events, drought frequency analysis and drought return period determination results through time-variant comprehensive drought indicator, the variation of drought hazard system presented the evident mitigation trend especially for light and moderate drought grades from north to south in Huaibei Plain. The typical drought processes of extreme grade occurred in 1966, 1978 and 2001, respectively, which is consistent with the statistics and observed data series of Anhui Province.

(3) The time-variant integrated drought indicator can better reveal drought evolution characteristic processes compared to single drought indicators. Determination of time-variant parameters as well as combination scenarios is crucial for derivation of time-variant PDFs of different single drought indicators. These are also difficult using the GAMLSS approach.

All in all, the time-variant integrated drought index can better reveal the uncertainty characteristics of the drought evolution process in comparison with different single drought indexes. Determination of time-variant parameters as well as combination scenarios is crucial for the derivation of time-variant PDFs of different single drought indicators. This is also a challenge for the establishment of comprehensive drought index. In addition, much work has been conducted concentrating on the reliability and applicability discussion of a new integrated drought index from the perspective of statistics, and further exploration focusing on the sensitivity analysis of time-variant drought index compared with other single drought indicators should be conducted in future work.