Copula-Based Hazard Risk Assessment of Winter Extreme Cold Events in Beijing

Extreme cold events (ECEs) have occurred more frequently over the last few winters in China, associated with large losses of human life and increasing costs. Here, copulas are used to establish a bivariate copula distribution model for ECE variables of duration and intensity, based on observed daily surface air temperatures in winter from 1978 to 2015 at 20 meteorological stations in Beijing. We demonstrate that durations of ECEs follow Weibull distributions, while their intensities fit a generalized extreme value distribution at most stations. The Gumbel–Hougaard copula best described the relationship between duration and intensity of ECEs at most stations. The joint and conditional return periods based on the bivariate copula described both ECE frequency and the corresponding hazard risk. A high risk was calculated for northern and western areas of Beijing, while a lower risk was calculated for urban and southeastern areas. Although the risk of a low temperature event of greater than 3 days with intensity in the range from −12 ◦C to −15 ◦C decreased, the risk of extreme low temperature events with durations greater than 2 days and intensity lower than −15 ◦C increased over the last 18 years. These probabilistic properties provide useful information for both climate change and hazard risk assessments.


Introduction
Arctic amplification (the greater warming of the Arctic compared with lower latitudes) is increasing the likelihood of the type of weather patterns that lead to mid-latitude cold extremes [1][2][3][4].Cold extremes, which are becoming more likely under global warming, are associated with large losses of human life and exponentially increasing costs [5][6][7][8].Thus, occurrences, characteristics, and risk evaluation of cold extremes have received increasing attention in recent decades [1,5,9,10].
Many previous studies have investigated the typical features and formation mechanisms of extreme cold events (ECEs) [7,11,12].For example, Choi,et al. [11] examined the spatial and temporal patterns of changes in extreme events of temperature over the Asia-Pacific Network countries.Cholaw, et al. [12] discussed large-scale circulation features of ECEs.Zuo, et al. [5] and Sun, et al. [10] investigated the characteristics of ECEs in China.Honda, et al. [13] and Screen, et al. [1] discussed the determinants of the future risk of cold extremes.Although previous studies have paid much attention to cold extremes, there are few studies dealing with their probabilistic behavior or risk assessments.
Identifying the probabilistic behavior and risk of an ECE is crucial to an in-depth understanding of climatic mechanisms responding to climate changes, as well as providing evidence for forecasting of future climate-induced risks.Risk assessment of extreme events is typically carried out using the peaks over threshold (POT) extreme value analysis [14], a generalized extreme value (GEV) [15] distribution, or a percentile threshold method [16].All these methods are probabilistic and based on univariate information, but can only be used for univariate statistical analysis of extreme temperature events.They cannot be used to analyze the multivariate characteristics of such events.An ECE is a multivariate phenomenon characterized by a temperature threshold, duration, and interval.It is difficult to make a comprehensive evaluation of the risk and probability of such cold extremes using only univariate information.The copula technique makes it possible to model dependence among multivariate data, independent of their marginal probability distributions.Thus, it allows multivariate random events to be described using different types of marginal distributions [17].The copula function was initially introduced to explore abnormally high risks associated with insurance and finance, but has been widely used in studies of extreme events, such as drought [18], rainfall [19], and floods [20].Previous studies have shown that the copula function can provide better multivariate joint distributions and an analysis of the frequency of related incidents.The corresponding return periods is a common criterion employed in the risk assessment of extreme events.The return period provides a very simple, yet efficient means for conducting risk analysis because it is able to concentrate a large amount of information into a single number [17].Few studies have applied copula functions and the return periods to study ECEs.
Beijing, the capital of China, is an international metropolis; it will host the 2022 winter Olympics.Concurrent with the rapid expansion of the city and the rapid increase in population, the intensity and frequency of ECEs have increased.Anomalous cold winters are important ECEs that have adverse effects on industrial and agricultural production, people's livelihoods, and social security.For instance, in 2004, a low temperature event occurred in Beijing in winter, which caused a severe heating energy shortage.The persistent extreme low temperature events in Beijing in November 2009 led to severe impacts on human life and the economy [21].Clearly, there is a strong need to identify characteristics of cold extremes and a risk assessment through comprehensive analysis.Here, we attempt to study ECEs in Beijing using copulas.The principal aim of this study was to evaluate the risk of ECEs in Beijing using copula functions and provide useful information for both climate change and risk assessments.The remainder of this paper is organized as follows.Section 2 describes the study area and the data used.The method used to establish a bivariate copula distribution model for ECE variables is presented in Section 2. The results of the copula and risk analysis for ECEs are given in Section 3. Finally, our conclusions are stated in Section 4.

Descriptions of Study Area and Meteorological Data
Beijing was selected as the study area.The city borders the North China Plain and the Inner Mongolia Plateau; it is surrounded by the Taihang Mountains on the west and the Yanshan Mountains to the north (Figure 1).Beijing is influenced by the Siberian anticyclone, which brings colder, windier, and drier winters.The minimum absolute temperature can be lower than −20 • C in bitterly cold winters.
The observational daily mean 2-m temperature (T2m) data from 20 meteorological stations in Beijing were used in this study (Table 1, Figure 1).The meteorological data are provided by the Chinese Meteorological Administration, which carries out data quality checks.The data for December, January, and February for the period 1978-2015 were used.There were a few missing values in the daily temperature data sets, although missing values comprised less than 0.01% of the total values.These missing values were filled using the average values of neighboring days.This gap filling method has no influence on long-term temporal trends.Locations and elevations of meteorological stations are given in Table 1.These stations are sparsely distributed in mountainous areas, but are densely distributed on the plain, especially within urban areas.The observational daily mean 2-m temperature (T2m) data from 20 meteorological stations in Beijing were used in this study (Table 1, Figure 1).The meteorological data are provided by the Chinese Meteorological Administration, which carries out data quality checks.The data for December, January, and February for the period 1978-2015 were used.There were a few missing values in the daily temperature data sets, although missing values comprised less than 0.01% of the total values.These missing values were filled using the average values of neighboring days.This gap filling method has no influence on long-term temporal trends.Locations and elevations of meteorological stations are given in Table 1.These stations are sparsely distributed in mountainous areas, but are densely distributed on the plain, especially within urban areas.

Selection of Extreme Cold Events
Exposure to −15 • C can result in frostbite, while long outdoor exposure to −10 • C can result in freezing of superficial tissues.In this study, ECEs were divided into low temperature events and extreme low temperature events.A low temperature event occurs when the daily mean temperature falls below −10 • C but is above −15 • C, while an extreme low temperature event occurs when the daily mean temperature falls below −15 • C. In previous studies, cold extremes were described by few parameters: frequency, duration or length, severity, and/or intensity [22][23][24].The principal aim of this study was the risk assessment of ECEs that can provide useful information to help the 2022 Winter Olympics.Therefore, the ECEs with long lasting times and extreme low temperatures should be paid more attention.In this study, for each ECE, we assigned a duration and an intensity.Duration refers to the number of days the event lasted, while intensity is computed as the mean temperature over the entire event.
Subsequently, over the period from 1978 to 2015, we computed the frequency of occurrence, maximum duration, and average intensity of ECEs at each station.The maximum duration refers to the maximum time these events lasted, while all averaged quantities refer to the average value computed for events occurring over an entire year.

Estimation of Marginal Probability Distributions
Marginal distributions of ECE duration and intensity were analyzed using a conventional univariate approach for each station.Six distributions were employed, including normal, gamma, generalized extreme value (GEV), extreme value (EV), three-parameter lognormal, and Weibull distributions.Parameters for these distribution functions were estimated using maximum likelihood estimation.The goodness-of-fit of the probability function was evaluated using the Kolmogorov-Smirnov (K-S) test [25] at a >90% confidence level.The probability distribution function having the best goodness-of-fit with observational data was selected for each station.

Copula Functions
Sklar [26] showed that if two random variables x and y follow arbitrary marginal distribution functions F X (x) and F Y (y), then there exists a copula C to combine these two marginal distribution functions, as follows: where the function F X,Y (x, y) is a two-dimensional distribution function with marginal distributions F X (x) and F Y (y).If the marginal distribution F X (x) and F Y (y) are continuous, then the copula function C is unique, and the joint probability density function can be written as follows: where f X (x) and f Y (y) denote the density functions corresponding to F X (x) and F Y (y), respectively.In this case, c is the density function of C, expressed as follows: where u and v denote the cumulative distribution functions, that is, F X (x) and F Y (y), in Equations ( 1) and ( 2), the values of which range from zero to one.
Copulas provide a convenient way to separately fit each variable to a marginal distribution and then join them together.Multivariate normality is only one option in a wide range of copula-based models that can capture the principal features of variable data, such as non-symmetry, nonlinear dependence, or heavy-tail behavior.There are numerous different copulas to choose from, with various correlation properties such as symmetry, tail dependence, and range of dependence [27].

Bivariate Archimedean Copulas
In this study, the Archimedean copulas (Table 2), including Clayton, Gumbel-Hougaard (GH), and Frank copulas, were selected to analyze the joint probability of ECEs because of their simplicity and wide representation [28,29].The copula parameter θ is used to measure the degree of association between u and v.The inference functions for the marginal distribution method (IFM) proposed by Joe (1997) [30] employed for fitting copulas involves estimating the parameters of the copula.The IFM is a fully parametric method, so the misspecification of the marginal distributions may affect the performance of this estimator.Thus, the joint distribution is estimated expeditiously by the marginal distributions and the parameters in the copula function are estimated using these maximum likelihood values individually.As noted by Joe (2005) [31], use of the technique can prevent computational difficulty for high dimensional models in the maximum likelihood estimation where the marginal distribution parameters and the copula parameters are estimated simultaneously.

. Goodness-of-Fit Tests for Copula Functions
Copula functions were evaluated using an empirical copula function [28].This empirical copula was defined as follows: where x i , y i (i=1, 2, . . ., n) are samples from a two-dimensional distribution; n is the sample size; and (X, Y), F n (x), and H n (y) are the empirical distribution functions of X and Y, respectively.Here, The root-mean-square error (RMSE) was used to estimate the degree-of-fit between the empirical copula and the theoretical copula [32].The RMSE can be written as follows: where n is the sample size, C p denotes the computed values of the theoretical copula, and C e denotes the observed values of the probability obtained from the empirical copula.The copula function with the lowest RMSE value was selected at each station.

Return Periods
The return period and related probability distributions were derived for risk assessment.In this study, the return periods for a bivariate distribution can be defined in two ways.The first method defines the joint return periods using one random variable equaling or exceeding a certain magnitude and/or another random variable equaling or exceeding another given magnitude.The second method defines the conditional return periods for one random variable given that another random variable equals or exceeds a specific value [33].
The two variables describing both types of ECEs are duration and intensity.Low temperature/extreme low temperature events with long durations and low temperatures may have a considerable influence on human society.In this case, the joint events considered in this study are the following: Let D denote the duration and T the intensity, then joint return periods related to these two events can be computed using the following: where ∧ denotes "and", and R {D>d∧T<t} denotes the joint return period when D exceeds the specific threshold and T is less than the specific threshold; P(D > d ∧ T < t) denotes the probability of this event; and G(t) and F(d) are the marginal distribution functions of intensity and duration, respectively.In this case, µ denotes the mean time interval between two successive events, with the unit of µ defined as a year in this study [17].
The conditional probability of ECEs describes the possibility of certain extreme cold episodes.The conditional return period of ECEs, when the given duration exceeds a certain threshold, is given by the following: and the conditional return period of ECEs with a given intensity threshold is expressed as follows: Return periods were interpolated to cover the whole area of Beijing using inverse distance weighted (IDW) interpolation [34].The risks of low temperature and extreme low temperature events, respectively, were evaluated in this study.To analyze the changes in risk characteristics for these two types of ECEs, data from 1978 to 2015 of each event were divided into two periods, covering from 1978 to 1998 and from 1999 to 2015.

Statistics of Low Temperature Events
The frequency, maximum duration, and mean intensity of low temperature events for each station are presented in Table 3.Clearly, there are regional difference in frequency.Zhaitang station, with 666 ECE occurrences, had the highest frequency.Generally, low temperature events occurred more often during the period 1978-1998 than during the period 1999-2015.For example, at Yanqing station, 393 events occurred before 1998, while 263 events occurred after 1998.
The maximum duration of ECEs at each station in Beijing averaged 10 d.The longest events were recorded at Miyun station, with a low temperature event lasting 14 days in 1978-1998, and at Pinggu station, with a low temperature event lasting 14 days in 1999-2015.There was little difference in the mean intensity of ECEs among stations.For example, the mean intensity of 148 events at Shunyi station from 1978 to 1998 was −11.53 • C, while the mean intensity of 170 events at Haidian station from 1978 to 1998 was −11.35 • C. For all stations, the mean intensity of events for the period 1999-2015 was lower than for the period 1978-1998.Both frequency and mean intensity were higher from 1978 to 1998 than from 1999 to 2015, which suggests that frequency and intensity of winter low temperature events are in decline.Yanqing station (ID 54506) was selected as a case study, with the aim to verify computations using the copula function in this study.Yanqing was selected for illustration because some new venues planned for the 2022 Winter Olympic will be built here and the ECEs are important for the Olympic Games.The goodness-of-fit of the marginal distributions for duration and intensity of low temperature events at Yanqing station were evaluated (Figure 2).Using the K-S test, the optimal probability distributions for duration and intensity for 20 stations in Beijing were determined.These distributions are listed in Table 4.
There were 393 and 263 low temperature events recorded at Yanqing station during the periods 1978-1998 and 1999-2015, respectively.At a 90% confidence level, the K-S test yielded D values of 0.0822 and 0.1005, respectively [25].To model the intensity of these events, all probability distributions performed well, although the GEV distribution best fitted the observed data.However, duration of these events was only modeled by the Weibull distribution, based on the K-S test (Figure 2).Similar results were found for the other 19 stations in Beijing (Table 4).Typically, for event duration, only one or two distributions passed the K-S test, with marginal distributions of duration best fitted by the Weibull distributions at the majority of stations.At a few stations, such as Miyun and Daxing, only the normal distribution passed the K-S test; while at Xiayunling station, only the gamma distribution passed the K-S test for the period 1978-1998.In the case of intensity, the GEV distribution performed well for most stations.At the Haidian station, an EV distribution performed better than other distributions in both study periods.In general, the marginal distributions of duration and intensity of low temperature events were best fitted by the Weibull and GEV distributions, respectively.
After selecting the best-fit marginal distribution functions, the parameters of each copula function were computed.The copula families were determined based on RMSE values.The copula that produced the least error (lowest RMSE) was selected (Table 5).According to Table 5, the GH copula provided the best-fit joint distribution for duration and intensity at most stations.The Frank copula was the best-fit for only three stations during the period 1978-1998, but four stations during the period 1999-2015.The Clayton copula distribution was selected for few stations, including Foyeding and Tanghekou stations, suggesting that the distributions of low temperature events may be regionally specific.Similar results were found for the other 19 stations in Beijing (Table 4).Typically, for event duration, only one or two distributions passed the K-S test, with marginal distributions of duration best fitted by the Weibull distributions at the majority of stations.At a few stations, such as Miyun and Daxing, only the normal distribution passed the K-S test; while at Xiayunling station, only the gamma distribution passed the K-S test for the period 1978-1998.In the case of intensity, the GEV distribution performed well for most stations.At the Haidian station, an EV distribution performed better than other distributions in both study periods.In general, the marginal distributions of duration and intensity of low temperature events were best fitted by the Weibull and GEV distributions, respectively.For the marginal distributions of the intensity and duration of low temperature events, the joint cumulative distribution function for the Yanqing station is shown in Figure 3.In the next section, these joint distributions are used to predict the return periods of winter low temperature events.
After selecting the best-fit marginal distribution functions, the parameters of each copula function were computed.The copula families were determined based on RMSE values.The copula that produced the least error (lowest RMSE) was selected (Table 5).According to Table 5, the GH copula provided the best-fit joint distribution for duration and intensity at most stations.The Frank copula was the best-fit for only three stations during the period 1978-1998, but four stations during the period 1999-2015.The Clayton copula distribution was selected for few stations, including Foyeding and Tanghekou stations, suggesting that the distributions of low temperature events may be regionally specific.For the marginal distributions of the intensity and duration of low temperature events, the joint cumulative distribution function for the Yanqing station is shown in Figure 3.In the next section, these joint distributions are used to predict the return periods of winter low temperature events.

Return Period and Risk Analysis
Here, we consider only one scenario in our analysis of joint return periods, that is, the duration of events greater than 3 days duration with an intensity below −12 • C but above −15 • C. The results of the joint return period (Figure 4) suggest there is spatiotemporal variation to the risk of low temperature events within the study area.

plain.
Comparison of Figure 4(a) and 4(b) shows a general decrease in the risk of low temperature events after 1999, particularly in the central area of Beijing.Low temperature events were more prevalent during the period 1978-1998, with frequencies that were twice as high as during the period 1999-2015.Although there was still a high risk of low temperature events in Huairou, Miyun, and Pinggu from 1999 to 2015, the area of high risk has shrunk.The spatial patterns for conditional return periods are shown in Figures 5 and 6.Two scenarios were considered in our analysis of conditional return periods, that is, (1) with an intensity below −12 °C but above −15 °C, and given duration exceeding 3 d (Figure 5); and (2) a duration greater than 3 days with a given intensity threshold below −12 °C but above −15 °C (Figure 6).

The main concentration of low values for the conditional return period { 12| 3}
T D R   was over Huairou, Yanqing, and Mentougou stations (Figure 5), which means that events with an intensity in the range from −12 °C to −15 °C with a duration exceeding 3 days occurred frequently in these regions.The northern parts of Beijing experienced a shorter conditional return period compared Larger joint return periods imply a smaller probability that low temperature events will occur, and vice versa.Figure 4a,b show similar spatial patterns to return periods for both study periods.The highest low temperature risk was calculated for areas north and west of Beijing.Return periods exhibited a northeast-southwest trend, consistent with the topography of the Beijing region (Figure 1).This indicates that under the influence of topography, low temperature events occur more frequently in the northern and western mountainous areas of Beijing and less frequently on the plain.
Comparison of Figure 4a,b shows a general decrease in the risk of low temperature events after 1999, particularly in the central area of Beijing.Low temperature events were more prevalent during the period 1978-1998, with frequencies that were twice as high as during the period 1999-2015.Although there was still a high risk of low temperature events in Huairou, Miyun, and Pinggu from 1999 to 2015, the area of high risk has shrunk.
The spatial patterns for conditional return periods are shown in Figures 5 and 6.Two scenarios were considered in our analysis of conditional return periods, that is, (1) with an intensity below −12 • C but above −15 • C, and given duration exceeding 3 days (Figure 5); and (2) a duration greater than 3 days with a given intensity threshold below −12 • C but above −15 • C (Figure 6).
The main concentration of low values for the conditional return period R {T<−12|D>3} was over Huairou, Yanqing, and Mentougou stations (Figure 5), which means that events with an intensity in the range from

Statistics of Extreme Low Temperature Events
Descriptive statistics for extreme low temperature events are presented in Table 6.There is a strong spatial variability in their frequency.Extreme low temperature events occurred 411 times at Tanghekou station during the period from 1978 to 2015, while only 12 extreme low temperature events occurred elsewhere in Beijing during this period.There were large differences in the maximum duration among stations.For instance, there was an extreme low temperature event lasting 10 days at Tanghekou station in the period from 1978 to 2015.However, elsewhere in Beijing, The spatial distributions for R {T<−12|D>3} (Figure 5) and R {D>3|T<−12} (Figure 6) are relatively consistent, implying there was a high probability of concurrence of low temperature events in winter with longer duration and lower temperatures.The results for R {D>3|T<−12} showed that the southern areas of Beijing had an extended conditional return period, compared with northern areas.Hence, the northern areas of Beijing experienced recurrent events more frequently, compared with southern areas.Comparing Figure 6a,b, a general increase in conditional return period of R {D>3|T<−12} after 1999 is observed, particularly at Tongzhou station.However, Miyun, Yanqing, and Mentougou stations experienced shorter return periods during the study period 1999-2015 than during the period 1978-1998, which contrasts with the conditional return period for R {T<−12|D>3} .This indicates that the risk of a low temperature event lasting a long time was higher at Miyun, Yanqing, and Mentougou stations after 1998.
In addition, by comparing the spatial distributions for R {T<−12|D>3} (Figure 5) and R {D>3|T<−12} (Figure 6), it can be found that the conditional return period of R {D>3|T<−12} is smaller than that of R {T<−12|D>3} .It indicates that the risk of the first scenario corresponding the conditional return period of R {T<−12|D>3} (Figure 5) happening is lower than the risk of the second scenario corresponding to the conditional return period happening.

Statistics of Extreme Low Temperature Events
Descriptive statistics for extreme low temperature events are presented in Table 6.There is a strong spatial variability in their frequency.Extreme low temperature events occurred 411 times at Tanghekou station during the period from 1978 to 2015, while only 12 extreme low temperature events occurred elsewhere in Beijing during this period.There were large differences in the maximum duration among stations.For instance, there was an extreme low temperature event lasting 10 days at Tanghekou station in the period from 1978 to 2015.However, elsewhere in Beijing, extreme low temperature events lasted for a maximum of 2 days.For all stations, the mean intensity of events for the period 1978-2015 was about −16 • C.Although the frequency of events during the period 1978-1998 was higher than for the period 1999-2015 for all stations.In contrast, the mean intensity of events was lower for the period 1978-1998 than for The best copula model was selected and listed for each station in Table 8.Generally, the GH copula was selected as the best copula to represent the relationship between extreme low temperature intensity and duration.However, the Frank copula yielded the smallest RMSE value for nine stations during the period 1978-1998 and for two stations during the period 1999-2015.The Clayton copula was the best-fit for only one station during the period 1978-1998 but three stations during the period 1999-2015.The best-fit copula-based joint distributions of Yanqing station for both study periods are plotted in Figure 8.Clearly, the copula distribution reflects the dependence of correlated extreme low temperature variables.Comparing Figures 3 and 8, it can be found that the best-fit copula-based joint distributions of Yanqing station for two types of ECEs are different.The Clayton copula and GH copula were the best-fit for winter low temperature events for two time periods, respectively.However, the GH copula was the best-fit for winter extreme low temperature events for both two time periods.
The best-fit copula-based joint distributions of Yanqing station for both study periods are plotted in Figure 8.Clearly, the copula distribution reflects the dependence of correlated extreme low temperature variables.Comparing Figure 3 and Figure 8, it can be found that the best-fit copula-based joint distributions of Yanqing station for two types of ECEs are different.The Clayton copula and GH copula were the best-fit for winter low temperature events for two time periods, respectively.However, the GH copula was the best-fit for winter extreme low temperature events for both two time periods.

Return Period and Risk Analysis
The return period calculations were supported by IDW interpolation mapping, as shown in Figures 9-11. Figure 9 shows the spatial patterns of joint return periods, corresponding to events of duration longer than 2 days and intensity below −15.5 °C.The spatial distributions of the return periods varied between study periods.Short return periods were mainly concentrated in northern and western areas of Beijing, especially in mountainous areas.Long return periods mainly occurred in urban centers.During the period 1978-1998, the return period varied from 0.2 yr to almost never.However, during the period 1999-2015, the return period varied from 0.15 to 45 yr, showing a general increase in the risk of extreme low temperature events compared with the earlier period.During the period 1978-1998, extreme low temperature events were unlikely in urban centers, while during the period 1999-2015, these centers were likely to experience such an event once every 20 yr on average.

Return Period and Risk Analysis
The return period calculations were supported by IDW interpolation mapping, as shown in Figures 9-11. Figure 9 shows the spatial patterns of joint return periods, corresponding to events of duration longer than 2 days and intensity below −15.5 • C. The spatial distributions of the return periods varied between study periods.Short return periods were mainly concentrated in northern and western areas of Beijing, especially in mountainous areas.Long return periods mainly occurred in urban centers.During the period 1978-1998, the return period varied from 0.2 year to almost never.However, during the period 1999-2015, the return period varied from 0.15 to 45 years, showing a general increase in the risk of extreme low temperature events compared with the earlier period.During the period 1978-1998, extreme low temperature events were unlikely in urban centers, while during the period 1999-2015, these centers were likely to experience such an event once every 20 years on average.
The spatial distribution of the conditional return periods corresponding to an intensity below −15.5 • C with a given duration exceeding 2 days is shown in Figure 10.From 1978 to 1998, the return period of R {T<−15.5|D>2}varied from 0.1 to 50 years.Short return periods were mainly distributed in northern areas of Beijing, with an average return period of less than 1 year.Meanwhile, southeast areas of Beijing had longer return periods of greater than 10 years.From 1999 to 2015, the return period of R {T<−15.5|D>2}varied from 0.1 to 12 years, with an average of 4 years.Short return periods were also mainly distributed in northern areas of Beijing, covering an expanded area compared with the earlier period.The return period in urban centers was less than 10 years, which is much shorter than the return period for the earlier study period.This indicates that events with an intensity below a threshold of −15.Extreme low temperature events were most likely to occur in suburbs, on average every 2 yr.The risk of extreme low temperature events was relatively low in all urban areas.By comparing Figure 10 and Figure 11, it can be found that the conditional return period of

Discussion and Conclusions
Quantifying natural disaster risks using multivariate indicators can lead to an assessment of major natural disasters that is more accurate, and allows preventive measures to be taken earlier to reduce losses of life and property [18][19][20].In this study, a methodology to establish probabilities and return periods for risk analysis and disaster mitigation using copulas was presented.A risk assessment of two types of winter ECEs in Beijing, namely low temperature and extreme low temperature events, was carried out.An event was described using two variables: duration and intensity.These variables were separately modeled by different marginal distributions and joined using copula functions.Based on the return period function, the joint and conditional return periods were calculated, and their corresponding spatial distributions were discussed.
Six marginal distributions were considered.Generally, the Weibull and GEV distributions provided the best goodness-of-fit for duration and intensity, respectively (Table 4 and Table 7).Three copula functions were applied to construct the bivariate joint probability distribution.The GH copula proved to be the best model for calculating the joint distribution of low temperature event variables (Table 5 and Table 8).Both the GH and Frank copulas were suited to construct the joint distribution of extreme low temperature event variables.
The spatial distributions of different return periods highlighted significant regional differences (Figures 4-6, 9-11).Generally, both low temperature and extreme low temperature events had short joint return periods and conditional return periods in northern and western mountainous areas of Beijing.Events with low temperature and long duration were less likely to occur in urban centers and southeastern areas of Beijing.The spatial distributions of return periods for low temperature and extreme low temperature events were approximately the same, but there were differences in the frequencies of these two kinds of events over the two study periods.These temporal changes suggest that the frequency, duration, and intensity of low temperature events have generally decreased since 1998.In contrast, extreme low temperature events have become longer and more intense in Beijing over the last 18 years, especially in southeastern areas of Beijing, which generally have had a low risk of extreme low temperature events.We found that Beijing experienced fewer events driven by low temperatures, but more events driven by extremely low temperatures over the most recent study By comparing Figures 10 and 11, it can be found that the conditional return period of R {T<−15.5|D>2} is smaller than that of R {D>2|T<−15.5} .This indicates that the risk of the first scenario corresponding to the conditional return period of R {T<−15.5|D>2}(Figure 10) happening is higher than the risk of the second scenario corresponding to the conditional return period of R {D>2|T<−15.5} happening (Figure 11), which is different from the comparison between the two conditional return periods for low temperature events.

Discussion and Conclusions
Quantifying natural disaster risks using multivariate indicators can lead to an assessment of major natural disasters that is more accurate, and allows preventive measures to be taken earlier to reduce losses of life and property [18][19][20].In this study, a methodology to establish probabilities and return periods for risk analysis and disaster mitigation using copulas was presented.A risk assessment of two types of winter ECEs in Beijing, namely low temperature and extreme low temperature events, was carried out.An event was described using two variables: duration and intensity.These variables were separately modeled by different marginal distributions and joined using copula functions.Based on the return period function, the joint and conditional return periods were calculated, and their corresponding spatial distributions were discussed.
Six marginal distributions were considered.Generally, the Weibull and GEV distributions provided the best goodness-of-fit for duration and intensity, respectively (Tables 4 and 7).Three copula functions were applied to construct the bivariate joint probability distribution.The GH copula proved to be the best model for calculating the joint distribution of low temperature event variables (Tables 5 and 8).Both the GH and Frank copulas were suited to construct the joint distribution of extreme low temperature event variables.
The spatial distributions of different return periods highlighted significant regional differences (Figures 4-6 and 9-11).Generally, both low temperature and extreme low temperature events had short joint return periods and conditional return periods in northern and western mountainous areas of Beijing.Events with low temperature and long duration were less likely to occur in urban centers and southeastern areas of Beijing.The spatial distributions of return periods for low temperature and extreme low temperature events were approximately the same, but there were differences in the frequencies of these two kinds of events over the two study periods.These temporal changes suggest that the frequency, duration, and intensity of low temperature events have generally decreased since 1998.In contrast, extreme low temperature events have become longer and more intense in Beijing over the last 18 years, especially in southeastern areas of Beijing, which generally have had a low risk of extreme low temperature events.We found that Beijing experienced fewer events driven by low temperatures, but more events driven by extremely low temperatures over the most recent study period, compared with data from prior to 1998.Such findings support a tendency toward more extreme events as global temperatures rise [35][36][37][38].Our results provide useful modeling for authorities given that increases in extreme low temperature events will likely become more frequent, despite global warming.
The purpose of this study was to use copula functions to establish a joint probability distribution of meteorological data.Using these distributions, we calculated the joint return period and conditional return period to establish risks related to ECEs.Our results have shown that calculations of the return period based on the copula function can act as an accurate and convenient framework for disaster assessment.Although the study period of 48 years considered herein is relatively short, it demonstrates the advantages of using a joint probability distribution.We plan to extend the time series of ECEs to continue to verify the accuracy of calculated return periods using a larger data set.

Figure 1 .
Figure 1.The Beijing study area showing the location of the 20 meteorological stations used herein.

Figure 1 .
Figure 1.The Beijing study area showing the location of the 20 meteorological stations used herein.

Figure 3 .
Figure 3. Joint cumulative distributions of intensity and duration for winter low temperature events at Yanqing station for the study periods: (a) 1978-1998, and (b) 1999-2015.
−12 • C to −15 • C with a duration exceeding 3 days occurred frequently in these regions.The northern parts of Beijing experienced a shorter conditional return period compared with southeastern parts.The conditional return periods of R {T<−12|D>3} during the period 1999-2015 were larger than those during the period 1978-1998.The largest differential between study periods occurred at Tongzhou station, with R {T<−12|D>3} values ranging from 0.35 to 1.77 years.Atmosphere 2018, 9, x FOR PEER REVIEW 11 of 21 with southeastern parts.The conditional return periods of { 12| 3} T D R   during the period 1999-2015 were larger than those during the period 1978-1998.The largest differential between study periods occurred at Tongzhou station, with { 12| 3} T D R   values ranging from 0.35 to 1.77 yr.

Figure 5 .R
Figure 5. Spatial distributions for the conditional return periods for low temperature events with a duration exceeding 3 days for study periods: (a) 1978-1998; (b) 1999-2015 (unit: year).

Figure 6 .
Figure 6.Spatial distributions of conditional return periods for low temperature events with a given intensity threshold of below −12 °C and above −15 °C for study periods: (a) 1978-1998; (b) 1999-2015 (unit: year).

Figure 6 .
Figure 6.Spatial distributions of conditional return periods for low temperature events with a given intensity threshold of below −12 • C and above −15 • C for study periods: (a) 1978-1998; (b) 1999-2015 (unit: year).

Figure 8 .
Figure 8. Joint cumulative distributions of intensity and duration of winter extreme low temperature events at Yanqing station for study periods: (a) 1978-1998, (b) 1999-2015.

Figure 8 .
Figure 8. Joint cumulative distributions of intensity and duration of winter extreme low temperature events at Yanqing station for study periods: (a) 1978-1998, (b) 1999-2015.
5 • C and a given duration exceeding 2 days occurred with increased frequently during the period 1999-2015.The spatial distributions of conditional return periods of R {D>2|T<−15.5} are shown in Figure 11 for both study periods.Clearly, the spatial distribution of the return periods of R {D>2|T<−15.5} changed significantly.The return periods from 1978 to 1998 were longer than the return periods for the study period 1999-2015.From 1978 to 1998, extreme low temperature events rarely occurred in most regions of Beijing, although they occurred once every 2 years at Yanqing, Huairou, and Miyun stations.In contrast, the return period varied from 0.15 to 40 years for the study period 1999-2015.Extreme low temperature events were most likely to occur in suburbs, on average every 2 years.The risk of extreme low temperature events was relatively low in all urban areas.

Figure 10 .
Figure 10.Spatial distributions of conditional return periods for extreme low temperature events of a given duration exceeding 2 days for study periods: (a) 1978-1998; (b) 1999-2015 (unit: year).

Figure 10 .
Figure 10.Spatial distributions of conditional return periods for extreme low temperature events of a given duration exceeding 2 days for study periods: (a) 1978-1998; (b) 1999-2015 (unit: year).

Figure 11 .
Figure 11.Spatial distributions of conditional return periods for extreme low temperature events with a given intensity threshold of below −15.5 °C for the study periods: (a) 1978-1998; (b) 1999-2015 (unit: year).

Figure 11 .
Figure 11.Spatial distributions of conditional return periods for extreme low temperature events with a given intensity threshold of below −15.5 • C for the study periods: (a) 1978-1998; (b) 1999-2015 (unit: year).

Table 1 .
Information on stations and their locations used in the study.

Table 3 .
Statistics of low temperature events for 20 meteorological stations in Beijing.
3.1.2.Joint Distribution of Winter Low Temperature Based on the Copula Function

Table 4 .
The goodness-of-fit based on the Kolmogorov-Smirnov (K-S) D statistic for marginal distributions of duration and intensity for low temperature events at 20 meteorological stations in Beijing.EV-extreme value; GEV-generalized EV.

Table 4 .
The goodness-of-fit based on the Kolmogorov-Smirnov (K-S) D statistic for marginal distributions of duration and intensity for low temperature events at 20 meteorological stations in Beijing.EV-extreme value; GEV-generalized EV.StationsThe

Table 5 .
Parameter estimation and root-mean-square error (RMSE) for selected copula distributions for 20 meteorological stations in Beijing.

Table 5 .
Parameter estimation and root-mean-square error (RMSE) for selected copula distributions for 20 meteorological stations in Beijing.

Table 6 .
Statistics of extreme low temperature events for 20 meteorological stations in Beijing.

Table 8 .
Parameter estimation and root-mean-square error (RMSE) for selected copula distributions for 20 meteorological stations in Beijing.

Table 8 .
Parameter estimation and root-mean-square error (RMSE) for selected copula distributions for 20 meteorological stations in Beijing.