Exploration of Copula Models Use in Risk Assessment for Freezing and Snow Events: A Case Study in Southern China

Abstract

Due to cold waves, low and extremely low temperatures occur every winter. Sudden cooling can cause freezing and snow disasters, which seriously affect transportation, power, safety, and other activities, resulting in serious economic losses. Based on precipitation and average temperature data from 258 national meteorological stations over the past 70 years, this study established a historical freezing and snow event data set, extracting the accumulated precipitation intensity (API) and accumulated temperature intensity (ATI). We selected the optimal distribution function and joint distribution function for each station and calculated the univariate and bivariate joint return periods. The return period accuracy plays an important role in risk assessment results. By comparing the calculations with the real return period for historical extreme events, we found that the bivariate joint return period based on a copula model was more accurate than the univariate return period. This is important for the prediction and risk assessment of freezing and snow disasters. Additionally, a risk map based on the joint return period showed that Jiangsu and Anhui, as well as some individual stations in the central provinces, were high-risk areas; however, the risk level was lower in Chongqing and the southern provinces.

1. Introduction

Global warming increases instability in the climate system, thereby increasing the frequency and intensity of extreme weather events [1,2,3]. Although the probability of extreme events is small, their impacts are serious, which is an important issue in current research. Currently, extreme high temperature events receive more attention than cold events, but recent weather extremes have shown an increase in cold air outbreaks and/or severe snowfalls across the Northern Hemisphere from 1990 to the recent past [4,5]. Moreover, global warming may exacerbate the cold events or lead to extremely cold winters in China [6]. However, cold wave disasters are frequently dismissed as having minimal probability, resulting in insufficient disaster prevention preparation, such as social underestimation of risks and lack of emergency equipment, ultimately increasing the risks [7]. Therefore, the study of extremely low temperature events in such areas should not be ignored.

However, extremely low temperature events have frequently occurred worldwide in recent years. The most recent example was the anomalously cold weather of January and February 2021 in Asia, Europe, and North America [2,6]. China has also suffered from extremely low temperature events, resulting in heavy losses [8]. Cold winter temperatures in southern China have become more frequent, while the preparation is insufficient, thereby resulting in more loss and increased risks. For example, in 2008, Southern China experienced the coldest winter in nearly 30 years. This disaster seriously impacted transportation, energy supply, power transmission, communication facilities, agricultural production, and daily life, causing RMB 151.65 billion of direct economic losses. Scholars estimated that the return period for this event was 50 to 100 years [9]. However, in January 2016, a more serious freezing and snow disaster occurred in southern China. According to the China Meteorological Administration, an area of 7.86 million km² experienced a temperature drop exceeding 6 °C from January 20 to 25 in 2016, changing the record of the southernmost boundary for snow [10]. From the perspective of the disaster formation mechanism in southern China, the anomaly of atmospheric circulation caused by the La Niña event led to a wide range of sustained low temperatures. At the same time, the cold air from the polar region meets the warm air from the tropical ocean in the middle and lower reaches of the Yangtze River; this is considered the main reason for freezing and snow events [11,12]. On the other side, the cold events are related to the impact of sea surface temperature anomalies in the tropical western Pacific and Northwest Atlantic on mid-high latitude and subtropical circulation anomalies [13]. Therefore, it is important to assess the risk of low-temperature disasters in these areas.

Analyzing the frequency and the probability of natural disasters facilitates a better understanding of their mechanisms. The return period has become an important tool for analyzing the possibility of occurrence and assessment of natural disasters. It can provide policy support for disaster prevention, which aids in building disaster reduction facilities [14], designing project standards [15,16], and making decisions on insurance and investment [17]. Additionally, it can aid decision-making departments and the public in understanding disaster risks and reducing disaster losses and impacts. The return period is widely used in meteorological, hydrological, and other natural disaster fields [18,19]. As the return period accuracy directly affects the risk assessment results, the method for calculating the return period is important. Various calculation methods of the return period showed different effects with disasters; therefore, it is important to select a more suitable method for disaster assessment. Traditional disaster studies used multivariate analysis based on linear correlation or common distribution functions [20]. However, the copula function model has been widely used in natural disaster analyses recently. It can flexibly describe the relationships among multivariate factors and has no marginal distribution form restrictions. Studies have shown that copula models can improve evaluation accuracy, and, thus, are widely used to assess the risks of drought, floods, earthquakes, sandstorms, rainstorms, and typhoons [10,21,22].

Owing to anthropogenic warming, freezing and snow disasters are easily ignored; therefore, many studies have focused on northern or northwestern China, where temperatures are lower and more freezing and snow disasters occur frequently [23]. However, after the disaster in Southern China in 2008, several studies focused on the mechanisms, characteristics, impacts, and risk assessment of such disasters in this area [11,24,25]. Although these studies assessed risks using various methods, none have evaluated the return periods based on the copula model. Therefore, we chose southern China as the representative to study freezing and snow events. We calculated the univariate and bivariate return periods based on historical meteorological data. We evaluated their spatial characteristics by comparing the calculated return periods with the average interval of historical events at each station. The discussion on the return period accuracy of freezing and snow disasters can provide a risk assessment method for countries and regions facing low temperature, freezing and snow disasters. Moreover, our results will provide references to take effective action for risk prevention and disaster reduction, help relevant departments accurately respond to disaster events, and support the construction of facilities and the formulation of emergency management standards.

2. Materials and Methods

2.1. Study Area

We chose 12 provinces (Anhui, Fujian, Guangdong, Guangxi, Guizhou, Hubei, Hunan, Jiangsu, Jiangxi, Shanghai, Zhejiang, and Chongqing) in southern China as the study area (Figure 1). These areas are economically developed and densely populated. In 2020, the population was 693 million, and the GDP was RMB 57.53 trillion, accounting for 49.08% and 56.62%, respectively, of the national total. Although the overall winter temperatures are mild, freezing temperatures, snow, and extremely low temperatures also occur. For example, a rare freezing and snow disaster occurred in 2008, and heavy snow and blizzards occurred in January 2016 and 2018 in the central and northern portions of the study area. Some provinces have experienced record snowfall in the last 10 years.

2.2. Data Sources and Factor Selection

This study selected 258 national meteorological stations in the study area (Figure 1). The data included daily precipitation and daily average temperature during the winter months (December, January, and February) from 1951 to 2020 and were downloaded from the China Meteorological Data Service Centre (http://data.cma.cn, accessed on 15 June 2021). The data were screened according to the data quality control code. Digital elevation data (“SRTMDEMUTM 90M”) was provided by the Geospatial Data Cloud site from the Computer Network Information Center of the Chinese Academy of Sciences (http://www.gscloud.cn, accessed on 22 June 2021).

We defined freezing and snow events according to the temperature and precipitation thresholds in previous studies. The daily average temperature of ≤0 °C and the daily precipitation of ≥0.1 mm were considered as freezing and snow days, and a freezing and snow event was further determined by adjacent freezing and snow days [24,26,27]. Based on this, a set of historical events was established for each station. Later, the accumulated precipitation (AP), accumulated temperature (AT), and duration (D) were extracted through each event, after which the accumulated precipitation intensity (API) and accumulated temperature intensity (ATI) for each event were calculated according to Equations (1) and (2).

(1)$\mathrm{API}=\frac{AP}{D}$

(2)$\mathrm{ATI}=\left|\frac{AT}{D}\right|$

where AP is the accumulated precipitation, AT is the accumulated temperature, and D is the duration of the event. Using the annual maximum sampling method, the annual extreme values of the above two indicators were extracted.

2.3. Methods

2.3.1. Fitting Distribution Functions

The annual extreme values of API and ATI were continuous random variables. Considering the entire area and different spatial properties in each region, we selected the best fitting distribution function for each station from among several common functions instead of adopting one distribution function for all stations. Thus, different functions were used to improve the accuracy of probability calculation. Nine distribution functions commonly used in the fields of hydrology and meteorology were selected, specifically gamma, exponential, lognormal, Weibull, normal, generalized extreme value, extreme value (Gumbel), Poisson, and generalized Pareto distribution functions. The optimal distributions for the factors were determined using the Kolmogorov–Smirnov test at a significance level of 0.05. The maximum likelihood method was used for the parameter estimations.

2.3.2. Construction of the Copula Model

The copula method, which was based on Sklar’s theorem [28], was first proposed in 1959 to explain nonlinear correlations among multiple factors. It can obtain a joint distribution by using multiple marginal distributions, and, thus, is closer to the actual events. Therefore, it has been widely used in probabilistic calculations and disaster risk assessments.

The copula function plays an important role in connecting the joint distribution functions of random variables (x, y) with their corresponding marginal distribution functions $F_X\left(x\right) \mathrm{and} F_Y\left(y\right)$, namely, copula function $C_{\theta }\left(U,V\right)$, which can be illustrated as follows:

(3)$F\left(x,y\right)=C_{\theta }\left(F_X\left(x\right), F_Y\left(y\right)\right),\forall x,y$

where $F_X\left(x\right) \mathrm{and} F_Y\left(y\right)$ are the corresponding marginal distribution functions of random variables (x, y), and C is the copula function with the parameter θ.

The Gaussian copula, t-copula, and Archimedean copula are widely used because of their simple calculations. Clayton, Frank, and Gumbel are the most common Archimedean copula. Thus, we chose these five types of copula functions for this study, and their formulas are presented in Table 1. The maximum likelihood method was used for parameter estimation, and the goodness-of-fit was calculated using the root mean square error (RMSE) and Akaike information criterion (AIC) [29]. The computation was completed in MATLAB. The minimum goodness-of-fit value indicates the best fitness. Later, the optimal copula function for each station was selected.

2.3.3. Calculations and Comparisons of Return Periods

We defined N as the time sequence length of the freezing and snow events and n as the number of events. Thus, E(L) is the average interval of events.

(4)$E\left(L\right)=\frac{N}{n}$

The univariate return period, T_API or T_ATI, is calculated as follows:

(5)$T_x=\frac{E\left(L\right)}{1-F\left(x\right)}, F_X\left(x\right)=P_X\left(X\le x\right)$

where T represents the return period, E(L) is the average interval of events, $F\left(x\right)$ is the corresponding marginal distribution, and X is the random variable (API and ATI).

The bivariate joint return period T_API&ATI, based on the copula function, can be expressed as

(6)$T\left(x,y\right)=\frac{E\left(L\right)}{P\left(X\ge x{\displaystyle \cup }Y\ge y\right)}=\frac{E\left(L\right)}{1-P\left(x,y\right)}=\frac{E\left(L\right)}{1-C\left(F_X\left(x\right),F_Y\left(y\right)\right)}=\frac{E\left(L\right)}{1-C\left(u,v\right)}$

where T represents the return period, E(L) is the average interval of events,$F_X\left(x\right)$ and $F_Y\left(y\right)$ are the corresponding marginal distribution, x and y represent the random variables API and ATI, and C is the copula function.

Using our historical freezing and snow event data set, we defined the 75% quantile of the API and ATI index values for each station as the threshold for extreme events and, thereby, extracted the extreme event database [30].

The average interval of the past extreme events was calculated as the return period (T_real) to check the relevant results [10,21]. The univariate return period for API (T_API) and ATI (T_ATI) and the joint return period for the two indicators (T_API&ATI) were calculated based on Equations (4)–(6). The error ratio between the theoretical return periods and the historical return period was calculated as follows to determine the accuracy of each method.

(7)$\mathrm{Error} \mathrm{Rate}=\frac{T-T_{real}}{T_{real}}\times 100\%$

where T represents the univariate/joint return period, and T_real is the average interval for the historical events.

3. Results

3.1. Spatial Distribution of Frequency and Duration

The historical spatial distribution of freezing and snow events was analyzed by calculating the historical average duration (Figure 2a), historical maximum duration (Figure 2b), and frequency (Figure 2c). Considering the average duration, regional differences among stations were small, and the average duration varied from 0 to 3.7 d. The duration was longer at higher elevations in the west and was shorter in Northeast Jiangsu, northwest Chongqing, southern Guangxi, and Guangdong. The average duration was approximately 1.5 to 2.0 d in the northern and central eastern regions of the study area. The duration was over 2.5 d for sixteen stations, including 12 in Guizhou, two in Jiangxi, and one each in Hunan and Hubei. Due to the relatively high elevation in Guizhou, this area was prone to lower temperatures for longer periods. The other four stations (Wufeng, Nanyue, Jinggangshan, and Lushan) were also high-elevation areas, with an average elevation of approximately 973.3 m. A comparison of the maximum duration revealed large differences in the values ranging from 0 to 33 d. The duration was over 20 d for 17 stations, including 8 in Guizhou, 4 in Hunan, and other stations in the central parts. Most stations that experienced an extremely long duration were in Guizhou and Hunan; this spatial distribution is also similar to the average duration.

There were significant regional differences in average annual frequency, ranging from 0 to 14.3. The spatial distribution was regular, except for a few high-value points, and it gradually decreased from north to south. The high-value stations were in the north and scattered throughout the center of the study area. There were 11 stations with an average of more than seven freezing and snow events per year. However, 69 stations had one or fewer yearly events; these were concentrated in Chongqing and the southern areas, mainly in Guangxi, Guangdong, and Fujian.

3.2. Construction of Distributions

3.2.1. Univariate Distributions

Among the selected distributions for the API data, the optimal distribution was the generalized extreme value distribution for 96 stations, gamma distribution for 37 stations, generalized Pareto distribution for 34 stations, and lognormal, Weibull, and other distributions for the remaining stations. For the ATI data, the optimal distribution was the generalized extreme value distribution for 93 stations, Weibull distribution for 51 stations, gamma distribution for 23 stations, and lognormal and normal distributions for the remaining stations. We randomly selected two stations as examples (Weining station and Puan station). The optimal distributions and parameter estimations for these two stations are listed in Table 2 and the fitness results are shown using a quantile-quantile plot in Figure 3.

3.2.2. Bivariate Joint Distribution

To construct the bivariate joint distribution, the parameters of the copula function were calculated, and a goodness-of-fit test was carried out using the RMSE and AIC. The relevant parameter estimation and goodness-of-fit results are listed in Table 3. Subsequently, the joint distribution was established, and the results are shown in Figure 4a. The optimal copula model, according to the goodness-of-fit test, was the Clayton copula for 108 stations, Frank copula for 66 stations, Gumbel copula for 30 stations, T-copula for 19 stations, and Gaussian copula for the remaining stations.

3.3. Calculation of the Return Period

The univariate return period was calculated using Equation (5), while the joint return period for API and ATI was calculated according to Equation (6); the result is shown in Figure 4b. Four different return periods were used: 5, 10, 30, and 50 a. Table 4 shows that the joint return period is always less than the univariate return period. As the univariate return periods increased, the bivariate joint return period also increased.

3.4. Comparison of Return Period Accuracies

The accuracy of the return period directly affects risk evaluations; therefore, improving the accuracy is crucial. Three types of return periods (real, bivariate, and univariate) were estimated to select the most suitable for assessment. Taking the Badong station as an example (Table 5), the 75% quantile of API and ATI were 4.13 mm and 1.58 °C, respectively, and extreme events occurred in 1955, 1956, 1957, 1964, 1984, and 1991. The time intervals between each adjacent event were 1, 1, 7, 20, and 7 years. Thus, the historical average return period was 7.20 a, which was considered the real return period. Additionally, the bivariate joint return period (T_API&ATI) and univariate return periods (T_API and T_ATI) were estimated as 6.75, 35.35, and 12.86, respectively. The error rates of univariate and bivariate return periods were 6%, 391%, and 79%, respectively.

Approximately three-quarters of the calculated joint return periods were less than the real return period, and the error rate was less than that of the univariate return periods. In terms of return period accuracy, 108 stations had a joint return period error rate of less than 50%, compared with 61 stations for API and 74 stations for ATI (Figure 5). There were 204 joint return periods with error rates of less than 100%, compared with 105 and 132 for API and ATI, respectively. Thus, the joint return period had a much lower error rate than that for the univariate types and can improve the return period accuracy of more stations. Comparing the univariate types, the return period error rate for ATI was less than that for API; thus, its accuracy was slightly higher. Therefore, the joint return period calculation including both API and ATI can greatly improve the accuracy of risk assessments for freezing and snow events.

3.5. Risk Analysis of Severe Freezing and Snow Events

The average API and ATI values for freezing and snow events in the study area were obtained, and the joint return periods for each station were calculated to evaluate the regional risk. The larger the return period was, the smaller the probability of extreme events. In contrast, the smaller the joint return period was, the higher the probability of the event. Based on this, the joint return period for each station was used to create a risk map (Figure 6). The map shows that the risk level in Jiangsu and Anhui, in the northern portion of the study area, is high. In addition, the risk levels of individual stations in central Hunan, Hubei, Guizhou, Jiangxi, and Fujian are high. Chongqing and the southern provinces in the study area have large return periods and thus, low probabilities. Except for Chongqing and a few central stations, the overall risk gradually decreases from north to south, which is consistent with existing research [31,32,33].

4. Discussion

Based on the temperature and precipitation data of 258 meteorological stations, this study defined a data set of freezing and snow events over the past 70 years. By determining the optimal univariate distribution and bivariate joint distribution model for each station, the return period was evaluated to assess risks. From the spatial distribution of historical freezing and snow events, the average elevation of the stations in Guizhou is approximately 958 m, which is significantly higher than that of other provinces; thus, the events lasted longer. Low frequencies and long durations characterized the freezing and snow events in Hunan; those in Jiangsu, northern Hubei, and Anhui had high frequencies and short durations. In the south, Guangxi, Guangdong, and Fujian provinces in the south had higher winter temperatures; thus, there were fewer freezing and snow events.

For the risk assessment results, Anhui, Jiangsu, and Hubei in northern areas and Guizhou, Hunan, and Jiangxi in central and western areas are at higher risk. This is similar to the results of previous studies [31,34]. Owing to latitude and climate, the temperatures in northern areas are much lower than those in southern areas, as shown through analysis of the winter temperatures. The average winter temperature from 1951 to 2020 in Jiangsu and Anhui were −0.6 °C and −0.18 °C, respectively, which are the only two provinces in the study area that have temperatures lower than 0 °C; therefore, the frequency of low temperatures below 0 °C in winter in these two areas is relatively high. Furthermore, the historical precipitation is relatively high [26] resulting in a high risk level. However, the average winter temperature (1951–2020) was 5.59, 6.48, 9.56, and 10.86 °C in Chongqing, Fujian, Guangxi, and Guangdong, respectively, where the risk is low. Another direct factor affecting the freezing and snow disaster is elevation, as Guizhou, Hunan, and Anhui have high elevations with more extreme events happening. The distribution of mountains and hills will also strengthen the impact of freezing and snow events [34]. La Niña is also one of the important reasons that may expand the freezing area. The anomaly of atmospheric circulation caused by La Niña led to a sustained low temperature, while the meeting of cold air from the polar region and the warm air from the tropical ocean in the study area provided plenty of moist air; the cold events or extreme disasters occurred in this manner [11,12].

When evaluating the goodness-of-fit for the univariate and joint distributions, this study selected the optimal distribution for each station, instead of selecting one distribution for all stations. This method is better for expressing the differences between regional stations and increases the accuracy of the probability calculation. The generalized extreme value distribution matched most stations, followed by the Weibull and gamma distributions. In the bivariate joint distribution, five types of copula models were considered, of which the Clayton copula best fit most stations, followed by the Frank copula and Gumbel copula. This provides a reference for the distribution selection of cold events, especially the freezing and snow events.

The existing risk assessment in China, especially in southern China, mainly relies on meteorological indicators such as temperature and precipitation. Some studies established comprehensive freezing and snow indexes through the analytic hierarchy process, fuzzy comprehensive evaluation, and compound weight method to analyze the risk distribution [26,34,35,36]. Referring to other disaster assessment methods, we assessed the freezing and snow events with the help of the copula method. By screening historical extreme freezing and snow disasters events and then comparing the univariate and joint return periods with the historical return period, it was found that the joint return period greatly improved the accuracy of the return period estimation for each station. The API and ATI are good indicators for freezing and snow risk assessments. Therefore, multivariate evaluation by way of the copula function should be considered when evaluating freezing and snow disaster events in the future.

Risk assessment for the study area was performed based on the joint return periods. Jiangsu and Anhui, in the northern study area, had small joint return periods and thus, a high probability of event occurrence and high risk level. At the same time, high-risk individual stations were distributed throughout Hunan, Hubei, Guizhou, Jiangxi, and Fujian. Chongqing, in the northwestern study area, and the southern provinces had large return periods and low probabilities of occurrence. Furthermore, the risk map based on joint return periods is useful for developing appropriate strategies for transportation, power supply, construction of facilities, and risk evaluation.

5. Conclusions

Owing to global warming, the frequency of extremely high temperature meteorological events is increasing, which has gained considerable attention. However, large-scale and high-intensity cold waves have frequently occurred worldwide recently; nevertheless, the prevention of such low-temperature disasters is generally neglected. Thus, when such disasters occur, they can cause serious losses, affect social production and economic activities, and even cause casualties. Using a joint return period with the copula model could aid regions or countries in accurately assessing the risks, providing a reference for disaster prevention and the corresponding construction of facilities, such as reinforcing traffic facilities and protecting power transmission lines in high-risk areas. It is also important to raise awareness of disaster prevention because enough public awareness and general knowledge of risk can help reduce the damage from the disaster. Further, it can provide scientific support to management authorities to improve disaster prevention and reduction and formulate risk response strategies. However, the distribution function selection, parameter estimation, return period calculation, and test methods need to be further improved. In future studies, we will continue to explore risk assessment methods.

Author Contributions

Q.L.: methodology, writing–original draft, and formal analysis; L.C.: data curation and investigation; Z.Y.: editing and validation; Y.X.: conceptualization, methodology, and supervision. All authors have read and agreed to the published version of the manuscript.

Funding

This study was funded by the National Natural Science Foundation of China (41977406) and the National Key Research and Development Program of China (2017YFC1502505).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

In this study, the data were mainly obtained from the China Meteorological Data Service Centre (http://data.cma.cn, accessed on 15 June 2021) and the Geospatial Data Cloud site from the Computer Network Information Center of the Chinese Academy of Sciences (http://www.gscloud.cn, accessed on 22 June 2021).

Acknowledgments

The authors thank the members of the disaster research group for their assistance.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Location of the study area and the meteorological stations. The red points represent the 258 stations.

Figure 2. Historical distributions of freezing and snow events in southern China. Distributions for (a) average duration, (b) maximum duration and (c) average frequency at each station.

Figure 3. Comparisons of the theoretical distribution function with the empirical distribution function for (a) API and (b) ATI. The data were from Anshun station in Guizhou Province.

Figure 4. Joint cumulative probability (a) and joint return period(b) for the API and ATI at Puan station.

Figure 5. Number of stations with different error rates for the three return period types.

Figure 6. Risk map for freezing and snow events. The historical average API and ATI data were used to calculate the joint return period to characterize the risk distribution.

Table 1

Copula functions used in this study.

Copula	Functions
Gaussian Copula	$C\left(u,v;\rho \right)={{\displaystyle \int }}_{-\infty }^{{\varnothing }^{-1}\left(u\right)}{{\displaystyle \int }}_{-\infty }^{{\varnothing }^{-1}\left(v\right)}\frac{1}{2\pi \sqrt{\left(1-{\rho }^2\right)}}\exp \left[-\frac{x^2-2\rho xy+y^2}{2\left(1-{\rho }^2\right)}\right]dxdy$
t-Copula	$C\left(u,v;\rho ,\delta \right)={{\displaystyle \int }}_{-\infty }^{t_{\delta }^{-1}\left(u\right)}{{\displaystyle \int }}_{-\infty }^{t_{\delta }^{-1}\left(v\right)}\frac{1}{2\pi \sqrt{\left(1-{\rho }^2\right)}}{\left[1+\frac{x^2-2\rho xy+y^2}{\delta \left(1-{\rho }^2\right)}\right]}^{-\left(\delta +2\right)/2}dxdy$
Clayton Copula	$C\left(u,v;\theta \right)={\left(u^{-\theta }+v^{-\theta }-1\right)}^{-\frac{1}{\theta }}$
Frank Copula	$C\left(u,v;\theta \right)-\frac{1}{\theta }\ln (1+\frac{(e^{-\theta u}-1)(e^{-\theta v}-1)}{e^{-\theta }-1}$)
Gumbel Copula	$C\left(u,v;\theta \right)\exp (-{[{(-\ln u^{\theta })}^{\theta }+{(-\ln v)}^{\theta }]}^{\frac{1}{\theta }})$

Where C represents the copula function with different parameters of ρ, δ and θ, x and y are the random variables, and u and v are the corresponding marginal distribution.

Table 2

Distributions functions of API and ATI for Weining and Puan Stations.

Variable	Station	Type	Distribution Function	Parameters
API	Weining	Generalized Extreme Value	$f=\left(\frac{1}{\sigma }\right)\exp \left(-{\left(1+k\frac{\left(x-u\right)}{\sigma }\right)}^{-\frac{1}{k}}\right){\left(1+k\frac{\left(x-u\right)}{\sigma }\right)}^{-1-\frac{1}{k}}$	k = 0.64μ = 0.91σ = 1.24
ATI				$k=$ −0.20$\mu =$ 1.29$\sigma =$ 4.91
API	Puan	Lognormal	$f=\frac{1}{x\sigma \sqrt{2\pi }}\exp \left\{\frac{-{\left(\log x-\mu \right)}^2}{2{\sigma }^2}\right\}$	μ = 1.08$\sigma =$ 0.79
ATI		Generalized Extreme Value	$f=\left(\frac{1}{\sigma }\right)\exp \left(-{\left(1+k\frac{\left(x-u\right)}{\sigma }\right)}^{-\frac{1}{k}}\right){\left(1+k\frac{\left(x-u\right)}{\sigma }\right)}^{-1-\frac{1}{k}}$	$k=$ −0.12$\mu =$ 0.84$\sigma =$ 2.06

Table 3

Parameters and goodness-of-fit tests for the copula function.

Station	Copula	Parameter Estimation	RMSE	AIC	Selected Copula
Puan	Gaussian Copula	0.2442	0.0211	−458.822	t-Copula
	t-Copula	0.2516	0.0210	−459.668
	Clayton Copula	0.2716	0.0221	−455.337
	Frank Copula	1.4752	0.0214	−459.145
	Gumbel Copula	1.1657	0.0233	−449.345

Table 4

Return periods defined separately by API and ATI, and the corresponding joint return periods.

Univariate RP	Weining Station		Puan Station		Panxian Station		Joint RP
	API/mm	ATI/°C	API/mm	ATI/°C	API/mm	ATI/°C	TARI&ATI	TAPI	TATI
5	3.52	6.58	5.70	3.22	4.89	2.85	2.84	3.01	2.81
10	5.81	7.24	8.07	3.73	8.45	3.61	5.39	5.75	5.30
30	12.20	8.08	12.50	4.42	16.91	4.86	15.66	16.78	15.29
50	17.08	8.40	14.87	4.70	22.50	5.48	25.93	27.82	25.29

Table 5

Error rate and comparisons among different return period accuracies.

Station	Extreme Event Threshold		Return Period				Error Rate/%
	API/mm	ATI/°C	Real	TARI&ATI	TAPI	TATI	TARI&ATI	TAPI	TATI
Badong	4.13	1.58	7.20	6.75	35.35	12.86	6%	391%	79%
Jingzhou	9.05	3.53	4.50	4.06	6.49	20.31	10%	44%	351%
Wuhan	11.20	4.45	4.50	3.90	6.79	10.79	13%	51%	140%
Laifeng	5.24	2.86	4.75	3.39	8.28	5.43	29%	74%	14%
Renhuai	3.60	1.52	5.57	5.357	12.89	10.00	4%	131%	80%

References

1. Ma, S.; Zhu, C.; Liu, J. Combined Impacts of Warm Central Equatorial Pacific Sea Surface Temperatures and Anthropogenic Warming on the 2019 Severe Drought in East China. Adv. Atmos. Sci.; 2020; 37, pp. 1149-1163. [DOI: https://dx.doi.org/10.1007/s00376-020-0077-8]

2. Cohen, J.; Agel, L.; Barlow, M.; Garfinkel, C.I.; White, I. Linking Arctic Variability and Change with Extreme Winter Weather in the United States. Science; 2021; 373, pp. 1116-1121. [DOI: https://dx.doi.org/10.1126/science.abi9167]

3. Cohen, J.; Screen, J.A.; Furtado, J.C.; Barlow, M.; Whittleston, D.; Coumou, D.; Francis, J.; Dethloff, K.; Entekhabi, D.; Overland, J. et al. Recent Arctic Amplification and Extreme Mid-Latitude Weather. Nat. Geosci.; 2014; 7, pp. 627-637. [DOI: https://dx.doi.org/10.1038/ngeo2234]

4. Overland, J.; Francis, J.A.; Hall, R.; Hanna, E.; Kim, S.J.; Vihma, T. The Melting Arctic and Midlatitude Weather Patterns: Are They Connected?. J. Clim.; 2015; 28, pp. 7917-7932. [DOI: https://dx.doi.org/10.1175/JCLI-D-14-00822.1]

5. Cohen, J.; Zhang, X.; Francis, J.; Jung, T.; Kwok, R.; Overland, J.; Ballinger, T.J.; Bhatt, U.S.; Chen, H.W.; Coumou, D. et al. Divergent Consensuses on Arctic Amplification Influence on Midlatitude Severe Winter Weather. Nat. Clim. Change; 2020; 10, pp. 20-29. [DOI: https://dx.doi.org/10.1038/s41558-019-0662-y]

6. Zheng, F.; Yuan, Y.; Ding, Y.; Li, K.; Fang, X.; Zhao, Y.; Sun, Y.; Zhu, J.; Ke, Z.; Wang, J. et al. The 2020/21 Extremely Cold Winter in China Influenced by the Synergistic Effect of La Niña and Warm Arctic. Adv. Atmos. Sci.; 2022; 39, pp. 546-552. [DOI: https://dx.doi.org/10.1007/s00376-021-1033-y]

7. Barnett, A.G.; Hajat, S.; Gasparrini, A.; Rocklöv, J. Cold and Heat Waves in the United States. Environ. Res.; 2012; 112, pp. 218-224. [DOI: https://dx.doi.org/10.1016/j.envres.2011.12.010]

8. Wang, C.; Yao, Y.; Wang, H.; Sun, X.; Zheng, J. The 2020 Summer Floods and 2020/21 Winter Extreme Cold Surges in China and the 2020 Typhoon Season in the Western North Pacific. Adv. Atmos. Sci.; 2021; 38, pp. 896-904. [DOI: https://dx.doi.org/10.1007/s00376-021-1094-y]

9. Ye, Q. Building Resilient Power Grids from Integrated Risk Governance Perspective: A Lesson Learned from China’s 2008 Ice-Snow Storm Disaster. Eur. Phys. J. Spec. Top.; 2014; 223, pp. 2439-2449. [DOI: https://dx.doi.org/10.1140/epjst/e2014-02218-7]

10. Feng, J.; Li, N.; Zhang, Z.; Chen, X. The Dual Effect of Vegetation Green-up Date and Strong Wind on the Return Period of Spring Dust Storms. Sci. Total Environ.; 2017; 592, pp. 729-737. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2017.02.028] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28336085]

11. Zhang, Z.; Gong, D.; Hu, M.; Guo, D.; He, X.; Lei, Y. Anomalous Winter Temperature and Precipitation Events in Southern China. J. Geogr. Sci.; 2009; 19, pp. 471-488. [DOI: https://dx.doi.org/10.1007/s11442-009-0471-8]

12. Liao, Z.; Zhai, P.; Chen, Y.; Lu, H. Differing Mechanisms for the 2008 and 2016 Wintertime Cold Events in Southern China. Int. J. Climatol.; 2020; 40, pp. 4944-4955. [DOI: https://dx.doi.org/10.1002/joc.6498]

13. Zhang, H.; Qin, J.; Li, Y. Climatic Background of Cold and Wet Winter in Southern China: Part I Observational Analysis. Clim. Dyn.; 2011; 37, pp. 2335-2354. [DOI: https://dx.doi.org/10.1007/s00382-011-1022-4]

14. de Silva, K.; Tanaka, N.; Yagisawa, J. Return Period of Flood Disturbance That Increases Diversity on Gravel Bars in Middle Stream of Rivers. Landsc. Ecol. Eng.; 2015; 11, pp. 1-18. [DOI: https://dx.doi.org/10.1007/s11355-013-0229-4]

15. de Biagi, V.; Barbero, M.; Borri-Brunetto, M. A Reliability-Based Method for Taking into Account Snowfall Return Period in the Design of Buildings in Avalanche-Prone Areas. Nat. Hazards; 2016; 81, pp. 1901-1912. [DOI: https://dx.doi.org/10.1007/s11069-016-2161-z]

16. Shiau, J.T. Return Period of Bivariate Distributed Extreme Hydrological Events. Stoch. Environ. Res. Risk Assess.; 2003; 17, pp. 42-57. [DOI: https://dx.doi.org/10.1007/s00477-003-0125-9]

17. Mochizuki, J.; Vitoontus, S.; Wickramarachchi, B.; Hochrainer-Stigler, S.; Williges, K.; Mechler, R.; Sovann, R. Operationalizing Iterative Risk Management under Limited Information: Fiscal and Economic Risks Due to Natural Disasters in Cambodia. Int. J. Disaster Risk Sci.; 2015; 6, pp. 321-334. [DOI: https://dx.doi.org/10.1007/s13753-015-0069-y]

18. Kao, S.C.; Govindaraju, R.S. A Copula-Based Joint Deficit Index for Droughts. J. Hydrol.; 2010; 380, pp. 121-134. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2009.10.029]

19. Alam, M.A.; Emura, K.; Farnham, C.; Yuan, J. Best-Fit Probability Distributions and Return Periods for Maximum Monthly Rainfall in Bangladesh. Climate; 2018; 6, 9. [DOI: https://dx.doi.org/10.3390/cli6010009]

20. Yue, S.; Ouarda, T.B.M.J.; Bobée, B. A Review of Bivariate Gamma Distributions for Hydrological Application. J. Hydrol.; 2001; 246, pp. 1-18. [DOI: https://dx.doi.org/10.1016/S0022-1694(01)00374-2]

21. Li, Y.; Gu, W.; Cui, W.; Chang, Z.; Xu, Y. Exploration of Copula Function Use in Crop Meteorological Drought Risk Analysis: A Case Study of Winter Wheat in Beijing, China. Nat. Hazards; 2015; 77, pp. 1289-1303. [DOI: https://dx.doi.org/10.1007/s11069-015-1649-2]

22. Zhang, L.; Singh, V.P. Bivariate Rainfall Frequency Distributions Using Archimedean Copulas. J. Hydrol.; 2007; 332, pp. 93-109. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2006.06.033]

23. Huang, X.; Deng, J.; Ma, X.; Wang, Y.; Feng, Q.; Hao, X.; Liang, T. Spatiotemporal Dynamics of Snow Cover Based on Multi-Source Remote Sensing Data in China. Cryosphere; 2016; 10, pp. 2453-2463. [DOI: https://dx.doi.org/10.5194/tc-10-2453-2016]

24. Wang, X.; Yang, F.; Gao, X.; Wang, W.; Zha, X. Evaluation of Forest Damaged Area and Severity Caused by Ice-Snow Frozen Disasters over Southern China with Remote Sensing. Chin. Geogr. Sci.; 2019; 29, pp. 405-416. [DOI: https://dx.doi.org/10.1007/s11769-019-1041-3]

25. Hu, A.; Xie, W.; Li, N.; Xu, X.; Ji, Z.; Wu, J. Analyzing Regional Economic Impact and Resilience: A Case Study on Electricity Outages Caused by the 2008 Snowstorms in Southern China. Nat. Hazards; 2014; 70, pp. 1019-1030. [DOI: https://dx.doi.org/10.1007/s11069-013-0858-9]

26. Wang, Y.; Wang, X.; Jiang, Z.; Ceng, X.; Ma, Q.; Cheng, C.; Liu, Q. Assessment and Zoning of Low Temperature, Rain/Snow and Freezing Disasters in China. Meteorol. Mon.; 2013; 39, pp. 585-591. (In Chinese)

27. Wang, L.; Gao, G.; Zhang, Q.; Sun, J.; Wang, Z.; Zhang, Y.; Zhao, S.; Chen, X.; Chen, Y.; Wang, Y. Analysis of the Severe Cold Surge, Ice Snow and Frozen Disasters in South China During January 2008: I. Climatic Features and Its Impact. Meteorol. Mon.; 2008; 34, pp. 95-100. (In Chinese)

28. Sklar, A. Fonctions de Repartition a n Dimensions et Leurs Marges. Publ. L’institut Stat. L’université Paris; 1959; 8, pp. 229-231.

29. Akaike, H. A New Look at the Statistical Model Identification. IEEE Trans. Autom. Control; 1974; 19, pp. 716-723. [DOI: https://dx.doi.org/10.1109/TAC.1974.1100705]

30. Aihaiti, A.; Jiang, Z.; Zhu, L.; Li, W.; You, Q. Risk Changes of Compound Temperature and Precipitation Extremes in China under 1.5 °C and 2 °C Global Warming. Atmos. Res.; 2021; 264, 105838. [DOI: https://dx.doi.org/10.1016/j.atmosres.2021.105838]

31. Gao, J. Analysis and Assessment of the Risk of Snow and Freezing Disaster in China. Int. J. Disaster Risk Reduct.; 2016; 19, pp. 334-340. [DOI: https://dx.doi.org/10.1016/j.ijdrr.2016.09.007]

32. He, F.; Shao, B. Risk Regionalization of Meteorological Disasters Caused by Cryogenic Freezing Rain and Snow in Shanghai. Sci. Meteorol. Sin.; 2011; 31, pp. 33-39. (In Chinese)

33. Hu, A.; Li, N.; Zhu, Y.; Wu, J.; Li, C. Integrated Risk Governance Paradigm for Meteorological Disasters: Thoughts on the Low-Temperature Freezing Rain and Snow Disaster in Southern China of 2008. Prog. Geogr.; 2010; 29, pp. 159-165. (In Chinese)

34. Mao, S.J.; Li, D.L. Comprehensive Assessment of Low Temperature, Snow and Freezing Weather in Southern China Based on Meteorological Elements. J. Glaciol. Geocryol.; 2015; 37, pp. 14-26. (In Chinese)

35. Wei, Y.; Wang, S.; Fang, Y.; Nawaz, Z. Integrated Assessment on the Vulnerability of Animal Husbandry to Snow Disasters under Climate Change in the Qinghai-Tibetan Plateau. Glob. Planet. Change; 2017; 157, pp. 139-152. [DOI: https://dx.doi.org/10.1016/j.gloplacha.2017.08.017]

36. Han, S.; Guo, W.; Ke, C.; Zhao, J.; Ma, D. Comprehensive Assessment of Snow Disaster Risk in Northeast China. Gaojishu Tongxin/Chin. High Technol. Lett.; 2020; 30, pp. 91-100. (In Chinese) [DOI: https://dx.doi.org/10.3772/j.issn.1002-0470.2020.01.011]