Joint Modeling of Severe Dust Storm Events in Arid and Hyper Arid Regions Based on Copula Theory: A Case Study in the Yazd Province, Iran

Natural disasters such as dust storms are random phenomena created by complicated mechanisms involving many parameters. In this study, we used copula theory for bivariate modeling of dust storms. Copula theory is a suitable method for multivariate modeling of natural disasters. We identified 40 severe dust storms, as defined by the World Meteorological Organization, during 1982–2017 in Yazd province, central Iran. We used parameters at two spatial vertical levels (near-surface and upper atmosphere) that included surface maximum wind speed, and geopotential height and vertical velocity at 500, 850, and 1000 hPa. We compared two bivariate models based on the pairs of maximum wind speed–geopotential height and maximum wind speed–vertical velocity. We determined the bivariate return period using Student t and Gaussian copulas, which were considered as the most suitable functions for these variables. The results obtained for maximum wind speed–geopotential height indicated that the maximum return period was consistent with the observed frequency of severe dust storms. The bivariate modeling of dust storms based on maximum wind speed and geopotential height better described the conditions of severe dust storms than modeling based on maximum wind speed and vertical velocity. The finding of this study can be useful to improve risk management and mitigate the impacts of severe dust storms.


Introduction
Evidence suggests that the rate of natural disasters has increased worldwide from 1993 to 2002. Of 2654 natural disasters that occurred during this time period, floods and dust storms were the most common types, accounting for 70% of the total, with drought, landslides, fires, etc., making up the remainder [1].
Dust storms, the result of particle transfer through the air, are an environmental problem and serious natural disaster in arid and semiarid regions. In general, dust storms are caused by wind turbulence, severe winds, and rapid ascents of dust to the upper atmosphere, which requires an increase in vertical and horizontal velocity over the area. The climate and environmental conditions of arid and semiarid regions contribute to dust emissions. Indeed, dust is an important indicator for assessing the degree of desertification. upper atmosphere (the free troposphere) and near surface level. Interactions between the upper atmosphere and near surface levels are the main reason for the dust storm. Therefore, the parameters of different atmospheric levels can play an important role in dust storms. For this purpose, the near surface maximum wind speed, geopotential height, and vertical velocity at different atmospheric levels (500, 850, and 1000 hPa) were selected as possible predictors of the return period of dust storms. Wind is a dynamic factor and a necessary condition for the creation and movement of dust from the surface of the earth. Many studies have considered wind speed as one of the most important parameters in the creation of dust storms [26,27,32,33]. Geopotential height indicates the atmospheric circulation pattern at upper atmosphere. In order to create dust storms, a low-pressure air mass is required. But if this low-pressure system is coordinated with the low geopotential height at intermediate atmospheric levels, and the thermodynamic relationship is established, the most suitable conditions for creating dust storms will arise. Additionally, the instability caused by the vertical velocity associated with a low pressure system can enhance conditions of the dust storm.
Most studies of dust storms have been performed using univariate methods. The hypothesis pursued here is that since dust storm occurrence and severity depends on several variables, such as wind speed, horizontal visibility, dust storm duration, and vertical velocity parameters at different levels, a multivariate approach would be useful. Most studies in the field of multivariate analysis of natural disasters are about hydrological and meteorological events such as floods and droughts [25,[34][35][36][37][38]. Few multivariate studies have been conducted on dust storms [7,29,30], although more researchers have recently discovered the importance of multivariate analysis of dust storms. [31] investigated the factors affecting storm dust by using correlation analysis. [29] calculated the return period of 79 severe dust storms in China over 1990-2008 using bivariate Archimedean copula functions. In that study, wind speed and dust storm duration were considered as the effective parameters of dust storm. The results of the study showed that the bivariate approach modeled the dust storm return period better than a univariate one. In [7], the parameters used for dust storm multivariate analysis include the 500-hPa atmospheric longitudinal circulation index, maximum wind speed in the 10-m high near ground, and surface soil moisture. The results showed that the bivariate Frank copula was suitable for return period of less than 10 years and trivariate Frank copula was suitable to estimate return period of more than 10 years. [30] investigated the dual effect of plant growth season and severe winds on spring dust storm in China using copulas. In this study, two pairs of variables of maximum wind speed-length of dust storm and growing season start date-number of intense wind events were used for bivariate analysis. The results of this study showed that the return period on the basis of the second pair of variables was closer to reality. In the current study area of Iran, a study of bivariate return period of dust storm has been carried out by [39]. There, the dust storm return period was calculated based on wind speed and geopotential height using copula theory for Yazd province in the period of 1982-2014.
Therefore, treating dust storm occurrence as a random variable, the bivariate return period of dust storm events will be estimated on the basis of the pairs of maximum wind speed-geopotential height and maximum wind speed-vertical velocity for the period of 1982-2017 in Yazd province, central Iran. For this purpose, the copula theory will be used for joint bivariate modeling of dust storm occurrence as a function of the above predictors. This theory will model the structure of the correlation between variables to provide bivariate predictions of dust storm return period. Based on the return period of severe dust storms, an early warning system can be created.

Materials and Methods
Located in the middle of the Loot desert and the Central desert of Iran, Yazd province extends from 29 degrees 52 min to 33 degrees and 27 min north latitude and 52 degrees 55 min to 56 degrees and 37 min east longitude and covers an approximate area of 131,575 square kilometers with the population of 990,818 people ( Figure 1). Due to low precipitation, high temperature, and geographical location of the province, about half of its area is covered by desert lands, which are always exposed to wind erosion and dust storms, and it is considered one of the critical hotspots of wind erosion in Iran [40]. Dust storms occasionally cause darkening air in the area due to their high intensity. In this region, severe dust storms cause irreparable damage each year, especially in the agricultural sector. Indeed, the farmers suffer from dust events more than others. Additionally, road accidents due to decreasing vertical visibility; air pollution; and the destruction of the buildings, roads, and channels are other effects of the dust storm that occur in this area every year [41]. In addition to sand dunes, there are numerous deserts or salty plains in the province, including the Siah Kooh desert in the north, the Daranjir Desert in the east, Abarkuh Desert in the southwest, and Saghand Desert in the northeast. Meteorological data and statistics indicate that the frequency of dust phenomena, including dust storms, is high in the study area [39]. Hence, recognizing the phenomenon of dust and the return period of severe storms in this region is important for combating its harmful effects. To this end, in order to estimate the return periods of dust storms, the stormy days were extracted based on the definition of the Meteorological Organization during the statistical period from 1982 to 2017 [42]. Forty dust storm events were identified during this period. Maximum wind speed variables at 10 m above ground level; geopotential height; and vertical velocity at three levels of 500, 850, and 1000 hPa were extracted, respectively, from synoptic station data and NCEP (National Centers for Environmental Prediction) data for analysis. temperature, and geographical location of the province, about half of its area is covered by desert lands, which are always exposed to wind erosion and dust storms, and it is considered one of the critical hotspots of wind erosion in Iran [40]. Dust storms occasionally cause darkening air in the area due to their high intensity. In this region, severe dust storms cause irreparable damage each year, especially in the agricultural sector. Indeed, the farmers suffer from dust events more than others. Additionally, road accidents due to decreasing vertical visibility; air pollution; and the destruction of the buildings, roads, and channels are other effects of the dust storm that occur in this area every year [41]. In addition to sand dunes, there are numerous deserts or salty plains in the province, including the Siah Kooh desert in the north, the Daranjir Desert in the east, Abarkuh Desert in the southwest, and Saghand Desert in the northeast. Meteorological data and statistics indicate that the frequency of dust phenomena, including dust storms, is high in the study area [39]. Hence, recognizing the phenomenon of dust and the return period of severe storms in this region is important for combating its harmful effects. To this end, in order to estimate the return periods of dust storms, the stormy days were extracted based on the definition of the Meteorological Organization during the statistical period from 1982 to 2017 [42]. Forty dust storm events were identified during this period. Maximum wind speed variables at 10 meters above ground level; geopotential height; and vertical velocity at three levels of 500, 850, and 1000 hPa were extracted, respectively, from synoptic station data and NCEP 1 data for analysis.

Copula Theory
To model multivariate probabilities, [43] proposed the copula function. The copula function makes it possible to combine univariate distributions with different families to create a bivariate or multivariate distribution, taking into account the dependence between variables. In other words, the copula function (C (FX1(x1), FX2(x2),…,FXN(xn))) is a connection function for the association of random variables X1, X2, ..
1 National Centers for Environmental Prediction

Copula Theory
To model multivariate probabilities, [43] proposed the copula function. The copula function makes it possible to combine univariate distributions with different families to create a bivariate or multivariate distribution, taking into account the dependence between variables. In other words, the copula function (C(F X1 (x1), F X2 (x2), . . . , F XN (x N ))) is a connection function for the association of random variables X 1 , X 2 , . . . , X N with marginal functions F X 1 (x 1 ), F X 2 (x 2 ), . . . ., F X N (x N ) that and the parameter of θ is defined in Equation (1) [44].
The most important advantage of the copula functions is the lack of limitations on the form of the marginal distributions, so that these functions can describe nonlinear and asymmetric dependence on these variables [30], whereas other types of probabilistic distributions are modeled by the assumption that the marginal distribution structure of all the variables is identical, which can cause errors. To create a realistic copula function for multiple variables, the correlation between the variables should be known. In this research, Kendall's correlation coefficient was used to estimate the correlation within the pairs of maximum wind speed-geopotential height and maximum wind speed-vertical velocity. The value of the correlation coefficient partly determines the type of copula function. Accordingly, if the correlation between variables is positive, different types of Archimedean (Frank, Clayton, and Gumbel), elliptical (t, Gaussian, and Normal), and other types of Copula function families can be used. For negative correlation values, a smaller number of copulas can be used for modeling (Table 1). Table 1. Copula functions used to bivariate analysis of dust storm.

Archimedean Family
Frank Joint CDF Generator Function

. Estimation of the Parameters of the Copula Functions
In order to estimate the parameters of the copula functions, nonparametric and parametric methods may be used. One nonparametric estimation approach is to use the relationship between the generator function of each copula and the Kendall correlation coefficient (Equation (2)) [45]. In the parametric method, the copula parameter is estimated by maximizing a logarithmic likelihood function (Equation (3)), where c θ is copula density function, F is marginal distribution function, and x 1k , x 2k , . . . , x pk (k = 1, . . . , n) are random variables.
where τ is Kendall correlation coefficient and L is a likelihood function that can be maximized.

Selecting the Copula Function
In order to select the most suitable copula function to fit the dust storm variables, the joint empirical probability values (JEPV) of the pair of variables were calculated by the empirical copula (Equation (4)) and then were compared with the values obtained from the theoretical copulas (Archimedes and elliptical families). Ordinary least squares method (OLS) was used to compare empirical copulas with each of theoretical copula functions. The OLS method, based on the square of the difference between the empirical and theoretical values, can be used to choose the best function form and parameter values (Equation (5)). Two information criteria-(AIC) and (BIC)-were also used (Equations (6) and (7)) [46,47]. In Equations (4)-(7), u and v are the empirical probability values of the maximum wind speed and the geopotential height/vertical velocity, respectively; P ei is the empirical copula value; P i is the theoretical copula value; k is the number of model parameters; n is the number of observations; and L is the maximum log-likelihood function.

Analysis of Bivariate Dust Storm Return Period
The bivariate return period is based on the mean interval time between dust storm events (E(L)) and the joint cumulative distribution function (C(U, V)). Based on the definition of the return period, if N is the length of the study period, n is the number of events, L is the interval time between events, and E(L) is the mean interval time of an event (E(L) = N/n). In this way, the bivariate return period of a dust storm is obtained by the primary ((T AND UV )) and secondary (T OR UV ) modes according to Equations (8) and (9). Accordingly, the return period T AND UV when two variables of maximum wind speed and geopotential height/vertical velocity are considered simultaneously, the bivariate return period of dust storm is based on the maximum wind speed that exceeds a certain value and the geopotential height/vertical velocity that exceeds a certain value (U >= u and V >= v) is determined, and the return period T OR UV when the maximum wind speed exceeds a certain value or the geopotential height/vertical velocity exceeds a certain value (U >= u or V >= v) is determined.
Additionally, the univariate return period is calculated for the comparison with bivariate values based on Equation (10). In this regard, F X (x) is the marginal function of the one of the variables of maximum wind speed, geopotential height, and vertical velocity.

Determine the Marginal Functions Of' Dust Storm Variables
It is necessary to use the copula functions for correlated predictor variables. For this purpose, the correlation between the pairs of maximum wind speed-geopotential height and maximum wind speed-vertical velocity at the three levels of 500, 850, and 1000 hPa was determined by the Kendall coefficient on the day of the dust storm and one day before the dust storm. The Kendall coefficient showed that on dust storm days, maximum wind speed as the main factor has a negative and significant relation with the geopotential height at 500 hPa and has a positive and significant relation with the vertical velocity at 1000 hPa at 5% significant level (p value < 0.05). There is no significant correlation between variables on the day before the dust storm. Table 2 shows the Kendall correlation coefficients within pair of dust storm predictor variables. In order to determine the marginal functions of the pair of variables, more than 20 univariate distribution functions were fitted to the variables, and the most appropriate functions based on Kolmogorov-Smirnov, Anderson Darling, and chi square tests were selected. Based on the results, the Wakeby distribution function was identified as the most suitable fit for maximum wind speed data and generalized extreme value function (GEV) as the most suitable fit to the geopotential height of 500 hPa and the vertical velocity at 1000 hPa (Table 3).

Variables CDF Parameters
Maximum wind speed Wakeby

Choosing the Best Copula Function for Bivariate Modeling of Dust Storms
In this study, the negative correlation between the pair of maximum wind speed-geopotential height resulted in candidate copulas for bivariate modeling of the dust storm that included Frank, Gaussian, Rotated Clayton, Rotated Gumbel, Student t, and Joe. Similarly, based on the positive correlation between the pair of variables maximum wind speed-vertical velocity, the copula types considered were the Frank, Gumbel, Clayton, Gaussian, and Student t ( Table 1). The parameters of the copula functions were estimated by both parametric and nonparametric methods. The Student t copula parameters, of which there are two parameters, are estimated only by the parametric method. The nonparametric method, which is based on the relationship between the Kendall coefficient (t) and the generator function of each copula, applies only to functions with one parameter. In Tables 4  and 5, the values of the copula parameters are presented based on both parametric and nonparametric methods for the pair of dust storm variable. After estimating the parameters, all selected copula functions were fitted to a pair of dust storm variables and the best fit was obtained by comparing the values of the empirical copula function against each of the above copulas as well as the OLS, AIC, and BIC criteria for the bivariate modeling of dust storms. Based on the selection criteria, Student t and Gaussian copulas were selected for bivariate modeling because of having the lowest values of AIC, BIC, and OLS (Tables 4 and 5). Figure 2 shows the q-q plot of the Gaussian copula and Student t copula for the pair of dust storm variables. The high correlation between the empirical and theoretical copulas in each pair of variables expresses the appropriateness of the Student t and Gaussian functions for the bivariate modeling of dust storm. Figure 3 shows the probability density values of the selected copula functions (Student t and Gaussian) for each of variables. The joint probability values (Figure 3) for the dust storm variables, maximum wind speed, and geopotential height have a symmetric correlation in

Joint and Conditional Probability of Dust Storm
Using the fitted copula function, the probability of dust storm can be determined based on its effective parameters. Here, the joint probability is the probability that the parameters affecting the dust storm will exceed a certain value (P (U ≥ u, V ≥ v)). Knowing of this possibility could be helpful in assessing the risk of severe dust storm and creating a dust storm warning system. The joint probability of dust storm is defined by the following relation based on the copula theory [34]).
For example, for a dust storm event based on the pair of maximum wind speed-geopotential height, when wind speed and geopotential heights exceed 13 m/s and 5681 m, the probability of dust storm occurrence is 0.22. Similarly, based on the pair of variable, when the values exceed 13 m/s and -0.12, respectively, the wind speed-vertical velocity occurrence probability is 0.021 (Figure 4). Additionally, the conditional probability of the dust storms by the copula function was computed based on the threshold velocities of 6.5 m/s in the study area [40] (Equation (12)). Accordingly, the probability of dust storm occurrence by dividing the wind speed, according to the threshold velocities set at 3 levels (3, 5.6 and 13 m/s), was obtained conditionally ( Figure 5). As the results showed, the probability values of the dust storm are different for the threshold levels of wind speed. Based on the wind speed-geopotential height, with increasing wind speed, the probability of dust storm

Joint and Conditional Probability of Dust Storm
Using the fitted copula function, the probability of dust storm can be determined based on its effective parameters. Here, the joint probability is the probability that the parameters affecting the dust storm will exceed a certain value (P(U ≥ u, V ≥ v)). Knowing of this possibility could be helpful in assessing the risk of severe dust storm and creating a dust storm warning system. The joint probability of dust storm is defined by the following relation based on the copula theory [34]).
For example, for a dust storm event based on the pair of maximum wind speed-geopotential height, when wind speed and geopotential heights exceed 13 m/s and 5681 m, the probability of dust storm occurrence is 0.22. Similarly, based on the pair of variable, when the values exceed 13 m/s and −0.12, respectively, the wind speed-vertical velocity occurrence probability is 0.021 ( Figure 4). Additionally, the conditional probability of the dust storms by the copula function was computed based on the threshold velocities of 6.5 m/s in the study area [40] (Equation (12)). Accordingly, the probability of dust storm occurrence by dividing the wind speed, according to the threshold velocities set at 3 levels (3, 5.6 and 13 m/s), was obtained conditionally ( Figure 5). As the results showed, the probability values of the dust storm are different for the threshold levels of wind speed. Based on the wind speed-geopotential height, with increasing wind speed, the probability of dust storm occurrence will increase, while based on the wind speed-vertical velocity, the probability will decrease at higher wind speed. For example, when the wind speed is higher than 3 m/s (less than the threshold level) and the geopotential height and vertical velocity are less than a certain value (probability of occurrence equal to 0.2, F(v) = 0.2), the probability of dust storm occurrence are, respectively, 0.20 and 0.19. If the wind speed is higher than 6.5 m/s (threshold velocity) and the geopotential height and vertical velocity are less than a certain amount, the probability is 0.21 and 0.17, respectively. Additionally, if the wind speed is higher than the threshold speed in the study area (above 13 m/s) and the geopotential height and vertical velocity are less than the specified value, the probability is 0.37 and 0.90, respectively ( Figure 5).
Climate 2020, 8, x FOR PEER REVIEW 10 of 16 (above 13 m/s) and the geopotential height and vertical velocity are less than the specified value, the probability is 0.37 and 0.90, respectively ( Figure 5).

Bivariate Return Period of Dust Storm
According to equations (8) and (9), the joint return period of dust storm was calculated by fitting the Student t and Gaussian copula functions to the pair of variables. The dust storm return period, as defined, is a function of the interval time between dust storms in a study period, which is determined by the mean of these intervals. The results of calculating the return period and for the pair of variables are shown in Figure 6. The return period of most of the dust storm events based on the maximum wind speed-geopotential height in the study period is less than 1 year (Figure 7(a)) and based on the maximum wind speed-vertical velocity is less than 0.5 year (Figure 7(b)) Additionally, the return period of most events based on the pair of dust storm variables (maximum wind speed-geopotential height and maximum wind speed-vertical velocity) is less than 1 year (Figure 7, c and d). The maximum return period was based on the maximum wind speedgeopotential height for a severe dust storm event with a horizontal visibility of 500 m and a wind speed of 29 m/s on May 20, 2016, while the maximum return period on the basis of the maximum wind speed-vertical velocity was on February 20, 2003, which was a less severe dust storm,

Bivariate Return Period of Dust Storm
According to Equations (8) and (9), the joint return period of dust storm was calculated by fitting the Student t and Gaussian copula functions to the pair of variables. The dust storm return period, as defined, is a function of the interval time between dust storms in a study period, which is determined by the mean of these intervals. The results of calculating the return period T AND UV and T AND UV for the pair of variables are shown in Figure 6. The T AND UV return period of most of the dust storm events based on the maximum wind speed-geopotential height in the study period is less than 1 year (Figure 7a) and based on the maximum wind speed-vertical velocity is less than 0.5 year (Figure 7b) Additionally, the T OR UV return period of most events based on the pair of dust storm variables (maximum wind speed-geopotential height and maximum wind speed-vertical velocity) is less than 1 year (Figure 7c,d). The T AND UV maximum return period was based on the maximum wind speed-geopotential height for a severe dust storm event with a horizontal visibility of 500 m and a wind speed of 29 m/s on 20 May 2016, while the T AND UV maximum return period on the basis of the maximum wind speed-vertical velocity was on 20 February 2003, which was a less severe dust storm, with horizontal visibility of 1000 m and a wind speed of 22 m/s. The values of the univariate return period based on the maximum wind speed showed that the dust storm event has a maximum value on the same date as the T AND UV return period (20 May 2016). However, based on the geopotential height variable, the maximum return period for the dust storm event occurred on 19 June 2010 with a lower intensity and wind speed of 7 m/s, and for the vertical velocity, the maximum return period for the dust storm event occurred on 20 February 2003. Figure 6 shows univariate return period of the dust storm events based on each variable (maximum wind speed, geopotential height, and vertical velocity). The maximum joint return period of T OR UV was based on the maximum wind speed-geopotential height and maximum wind speed-vertical velocity on the same dates as the return period of T AND UV . Therefore, considering that the maximum return period obtained from the maximum wind speed-geopotential height matches the intensity of the dust storm events, it can be concluded that the bivariate modeling of dust storms based on maximum wind speed and geopotential heights is more suitable than that based on the maximum wind speed and vertical velocity. In other words, the return period resulting from the variables of maximum wind speed and geopotential height at the occurrence time of a severe dust storm conforms to the maximum wind speed recorded during the study period, while the return period resulting from the variables of maximum wind speed and vertical velocity at the time of the dust storm conforms with less intensity and lower wind speed. It can be concluded that dust storm modeling based on maximum wind speed and geopotential height can better describe the situation of severe dust storms in the study region. the return period resulting from the variables of maximum wind speed and geopotential height at the occurrence time of a severe dust storm conforms to the maximum wind speed recorded during the study period, while the return period resulting from the variables of maximum wind speed and vertical velocity at the time of the dust storm conforms with less intensity and lower wind speed. It can be concluded that dust storm modeling based on maximum wind speed and geopotential height can better describe the situation of severe dust storms in the study region.

Conclusions
Dust storms are occurring frequently as a natural disaster, partly due to global warming and desertification. Due to the repetition of this phenomenon in recent years, information about its return period is a key reference for risk assessment, as well as the forecast and warning of the medium and

Conclusions
Dust storms are occurring frequently as a natural disaster, partly due to global warming and desertification. Due to the repetition of this phenomenon in recent years, information about its return period is a key reference for risk assessment, as well as the forecast and warning of the medium and long-term risks [48]. The estimation of the return period of natural disasters, such as dust storm, will determine the course of the occurrence of this phenomenon in the future and allow for the assessment of its severity purposefully and quantitatively. In fact, calculating the return period can predict the occurrence time of a phenomenon on different time scales for the future.
In this study, the return period of 40 dust storm events during the statistical period was estimated using copula theory. Because the generator factors of dust storm have different probability distribution function, there is a nonlinear relationship between them. Therefore, copula theory was used to solve this problem. Bivariate modeling was performed based on the effective factors of dust storm events. These parameters were selected at two spatial levels (near-surface and upper-atmosphere levels). At near-surface level, the maximum wind speed was applied to joint analysis of severe dust storm events, while vertical velocity and geopotential height were considered as upper atmospheric level parameters.
There are few research papers that can compared with our finding in this study. However, the results of correlation analysis indicate that there is a significant negative relationship between maximum wind speed and geopotential height, and a positive and significant relationship between the maximum wind speed and vertical velocity. Based on the results of our best fit, the Student-t and Gaussian copula functions were used for bivariate modeling of storm storms. These findings were not similar to the results of [7,29], and [30], due to the different dust storm variables. Those selected Clayton and Frank copulas for bivariate modeling of dust storm.
However, our finding indicated that the bivariate return period of dust storm is closer to reality than univariate return period. These results are in agreement with [29], who showed bivariate return period of dust storm was smaller than the univariate return period.
Estimation of the occurrence time of dust storms in the future is very important to prevent the probability of hazards and for risk management. In this study, considering the random nature of dust storm, bivariate modeling was conducted to evaluate this phenomenon by copula functions.
The results of the dust storm return period showed that the values obtained from the pair of maximum wind speed-geopotential height can describe the dust storm conditions better than the pair of maximum wind speed-vertical velocity, because the occurrence time of severe dust storm and maximum wind speed, recorded at the study period, is almost similar. Generally, the finding of this study indicated that the joint modeling of dust storm events using the pair variable of maximum wind speed-vertical velocity can help provide the practical strategies to reduce the potential hazards as a consequence of severe dust storm. Finally, we propose using other copula functions such as Gaussian mixture copula to joint modeling of dust storm events. Due to nonlinear relationship between the effective parameter of dust storm event, the Gaussian mixture copula can model highly nonlinearity in the data [49].

Conflicts of Interest:
The authors declare no conflict of interest.