A 3D Copula Method for the Impact and Risk Assessment of Drought Disaster and an Example Application

Droughts have more impact on crops than any other natural disaster. Therefore, drought risk assessments, especially quantitative drought risk assessments, are significant in order to understand and reduce the negative impacts associated with droughts, and a quantitative risk assessment includes estimating the probability and consequences of hazards. In order to achieve this goal, we built a model based on the three-dimensional (3D) Copula function for the assessment of the proportion of affected farmland areas (PAFA) based on the idea of internally combining the drought duration, drought intensity, and drought impact. This model achieves the “internal combination” of drought characteristics and drought impacts rather than an “external combination.” The results of this model are not only able to provide the impacts at different levels that a drought event (drought duration and drought intensity) may cause, but are also able to show the occurrence probability of impact at each particular level. We took Huize County and Mengzi County in Yunnan Province as application examples based on the meteorological drought index (SPI), and the results showed that the PAFAs obtained by the method proposed in this paper were basically consistent with the actual PAFAs in the two counties. Moreover, due to the meteorological drought always occurring before an agricultural drought, we can get SPI predictions for the next month or months and can further obtain more abundant information on a drought warning and its impact. Therefore, the method proposed in this paper has values both on theory and practice.


INTRODUCTION
Droughts are a complex and recurrent natural disaster that result in widespread effects on humans and natural systems, including agriculture, ecosystems, energy, economics, and so forth [1][2][3][4][5][6][7][8]. A drought can be commonly classified into four types: meteorological, agricultural, hydrological, and socioeconomic, and the meteorological drought has received more attention as it usually occurs before the other three types of drought [1,9].
Over the past decades, droughts around the world have caused tremendous economic losses on an annual basis [10,11]; furthermore, droughts are a major cause of unexpected crop failure (in corn, rice, and wheat). In China, problematically frequent occurrences of drought have resulted in great impacts on social and economic development and human life [12,13], and historical disaster data show that the average annual area of crop damage due to drought accounted for 50% of all areas affected by meteorological disasters over the period of 2004-2013 [14]. Due to rising temperatures and increasing climatic variability [15][16][17][18], it is widely recognized that droughts may become more frequent and more serious globally in the coming decades on regional and temporal levels [17,19,20]. In this context, the challenges faced by agricultural production will become more severe. In addition, the possibility of soil desertification will increase under prolonged drought, so sand prevention and control also face major challenges [21,22]. Therefore, drought risk assessment, especially quantitative drought risk assessment, is required in order to reduce the negative impacts associated with droughts [23][24][25][26].
Quantitative drought risk assessment includes estimating the probability and consequences of a drought disaster [27]. To achieve this, we need on three levels to (i) characterize the drought statistically in terms of its relative severity; e.g., exceeding probability [28], (ii) quantify the consequences of a given drought event [29], and (iii) estimate the consequences of droughts with varying severity. The key point here is to separate "droughts" from "non-droughts" [30]. Although studies have proven that risk assessment results on the impacts or losses provided by statistical methods and physical models are reliable and valid [31][32][33][34][35][36], the accurate assessment of drought impacts on agricultural production is not easy to calculate, and there are still various data-related and methodological problems that need to be solved [37,38]. In particular, these include the requirement of (i) the time series of drought index with sufficient spatial and temporal details in order to obtain enough information on the local meteorological and hydrological conditions, and (ii) the accurate, high-resolution, long-term yield or economic loss databases [36]. Therefore, a complete description of a drought disaster requires multiple related variables, and the appropriate option is a multivariate analysis using the Copula function [39]. The Copula function is a statistical tool, which can used to construct a multivariate joint distribution function for analyzing the statistical characteristics of dependent variables. However, the applications of the Copula technique in assessing the impacts of droughts on agricultural production are very limited in published literature, and almost all research has been performed to obtain the qualitative relationships between drought characteristics and the degree of a potential disaster [32,[40][41][42][43][44][45][46]. Moreover, most previous studies use two-dimensional Copula functions, and crop damage data are not often used as a variable to be substituted into the Copula function, but are instead used as an independent variable in combination with the comprehensive features of a drought obtained by the Copula function.
The goal pursued by many international researchers is to obtain objective results and reduce subjective judgments through the use of statistical methods and physical models to provide more accurate drought information and warnings for the prevention of disasters and to reduce losses. It is of more practical significance to assess the specific impacts of droughts with different severities and durations, and to see how crop damage data are substituted directly into the Copula function as one of the variables in assessing the impact of droughts. However, there is no published work that uses the above ideas to study the impact of droughts, so it is, therefore, worth studying in further detail. Therefore, the objective of this paper is to make a breakthrough in this respect.
Combining the above reasons, in this paper we propose a three-dimensional (3D) Copula model for the assessment of drought risks in terms of agriculture. We built a model based on the three-dimensional (3D) Copula function for the assessment of the proportion of affected farmland areas for the purpose of quantitative drought risk assessment. We calculated the proportion of affected farmland areas (PAFA) based on the farmland areas affected by droughts and the annual agricultural planting area, and then obtained the PAFA time series caused by each drought event. The reason behind using PAFA as the assessment object is that the affected area is a result of drought; thus, the impacts of "droughts" can be clearly separated from those of "non-droughts." Next, we extracted drought events from a monthly standardized precipitation index (SPI) time series, which is calculated based on precipitation data and satisfies a standard normal distribution. After, we calculated the parameters of the duration and the severity of a drought, as well as the PAFA, and thus, the marginal distribution functions of duration, severity, and PAFA were obtained, and their joint distribution function was then computed using a three-dimensional Copula function. The joint probability and the classification of these three variables were both achieved using this 3D statistical model. Different types of drought events were identified based on the combination of these variables, and the probabilities of occurrence of the drought events were calculated. Thus, the estimation of the PAFA of a drought event could be obtained by integrating this method with the monitoring of current droughts and the prediction of future droughts.

Study Area
Yunnan Province, located in the southwest of China, is influenced primarily by the South Asian monsoon and also impacted by the East Asian monsoon, the Tibetan Plateau monsoon, and westerlies [47]. The soil in Yunnan Province is covered by a large area of karst landforms and characterized by poor water conservation. Due to the complicated land surface terrain, as well as the control and influence of different general circulation patterns, the seasonal and geographical distributions of precipitation are extremely uneven in Yunnan, and droughts have become one of the most extensive, frequent, and severe natural disasters in Yunnan Province [48][49][50]. Yunnan's crops mainly include rice, wheat, and corn, among others, and the crops are normally fed and cultivated with rain.
Because in this article we are trying to build a Copulabased three-dimensional (3D) risk analysis model, including the duration and intensity of the drought, as well as the PAFA, only the areas where the three types of data (SPI, affected farmland areas caused by droughts, and annual agricultural planting area) that meet the following three requirements can therefore be used as the research area in this article: (i) the daily precipitation series as long as possible without missing values, (ii) the number of disaster loss records (affected farmland areas caused by droughts) from 1984 to 2012 should be as much as possible, and (iii) the annual agricultural planting area from 1984 to 2012 requires as much continuous data as possible. Figure 1 shows the spatial distribution of the number of disaster loss records caused by droughts from 1984 to 2012 submitted by each county in Yunnan Province, and the stations that comply with the requirement in Yunnan Province are also shown in Figure 1. After a survey of all the stations (counties) in Yunnan Province, it was found that two stations (counties) met the above three requirements: Huize county and Mengzi county.
In this paper, we take Huize County and Mengzi County in Yunnan Province as examples to demonstrate the threedimensional (3D) Copula model for the assessment of drought risk on agriculture.

Data
The data used in this study on damages caused by droughts come from China's Meteorological Disaster Loss Databases at County Level (1984-2012, including county name, geographic location, category of meteorological disasters, starting and ending time of meteorological disaster, number of affected populations, affected farmland areas, and direct economic losses) compiled by the China Meteorological Administration (CMA). Currently, the database is stored in the National Climate Center, and applied to the scientific research [51] and operational works (http://10.28. 107.46:8084/MDMIS_oneMap/) of meteorological disaster risk management. The data on the affected farmland areas caused by drought was obtained between 1986 and 2012; of these, there are 20 records in Huize County and 21 records in Mengzi County. The two counties did not experience drought in 1984 and 1985.
It is known that agricultural yields vary continuously with the hydro-climatic conditions; it is also influenced by many other factors such as field management and specific varieties of crops. In addition, the direct economic loss of agriculture caused by droughts is also related to the current price of crops, and the factors affecting prices are so complex that they far exceed the drought itself. Therefore, in this study, we used the data of affected farmland areas caused by droughts because they have an intuitive cut-off between "drought" and ''non-drought" conditions.
The exposure of crops to drought will also change with changes in the agricultural planting areas. Even if a drought occurs with the same duration and intensity, the larger the agricultural planting area, the larger affected the farmland area is. This makes it difficult to compare and analyze the impact of droughts when the data on the affected farmland areas in different years are used directly. Thus, the PAFA in each year was used for analysis in this study. In order to cooperate with the data on the affected farmland areas, we also need the continuous data of annual agricultural planting areas to calculate the PAFA. The data on the annual agricultural planting areas were taken between 1986 and 2012 (Figure 2) in Huize and from 1995 to 2012 in Mengzi. These data were obtained from the yearbook of Huize and Mengzi. Because of the establishment of the marginal distribution functions of duration and intensity of drought, it is necessary to use a time series of drought index. The monthly SPI [52] was calculated after processing daily precipitation data into monthly data. It should be noted that the Gamma distribution is used as the distribution function of precipitation in the calculation of SPI. The Kolmogorov-Smirnov distribution test results show that all stations in Figure 1 have passed the distribution test. In order to obtain more accurate and stable marginal distribution functions of drought duration and drought intensity, the daily precipitation data without missing values from 1961 to 2012 from Huize County and Mengzi County were acquired from the National Meteorological Information Center of the China Meteorological Administration (http://www.nmic.cn).

Description of Drought Characteristics
We described drought events based on the two major characteristics: drought duration and severity, which were both extracted from monthly SPIs. The most applied index is the SPI that concerns meteorological drought. As a powerful index, the SPI has been widely used to analyze droughts in different parts of the world, and the advantage of the SPI over other indices is that the SPI depends only on precipitation, the flexibility of timescales at which this index can be calculated, as well as being comparable in time and space. The SPI obtained based on the precipitation of 1 month is called a monthly scale index, the SPI obtained based on the precipitation of 3 months is called a seasonal scale index, and the SPI obtained based on the precipitation of 6 months is called an semi-annual scale index. SPI for 1 month is suitable for describing meteorological drought, SPI for 3 months is suitable for describing agricultural drought, and SPI for more than 6 months is suitable for describing drought in watersheds or groundwater. Since the SPI can characterize short-and long-term droughts, various research on drought risk assessment has been carried out using the SPI, and the results show that the SPI is ideal for performing a risk assessment in comparison to other drought indices [53][54][55][56]. Recently, a piece of research applied six drought indices to estimate the drought onset, and the results showed that meteorological drought indices predict the onset of a drought earlier than hydrological and agricultural drought indices [57]. Therefore, in this paper, we use the SPI to allow the method to be more practical in terms of risk warnings of droughts.
For an SPI time series, the definitions of drought duration and drought intensity that are adopted in the run theory [58] are shown in Figure 3. Run theory is a method of time series analysis, which is widely used in the identification of drought events [44,[59][60][61]. Drought duration (d) represents the number of months during which the SPI index is continuously below the threshold value (s 0 = 0). Because the SPI of the monthly scale is used in this paper, the drought duration lasts for at least 1 month, and the range of the d value is within [1, + ∞] month. Drought severity refers to the accumulation of the value of the SPI index within the ranges of the drought duration, and can be calculated by using the formula s = − d i=1 SPI i , and the range of the s value is within (0, +∞). Thus, the duration (d 1, d 2, d 3 · · · ) and severity (s 1, s 2, s 3 · · · ) of drought events can be extracted from an SPI time series [62,63]. Based on previous studies, the classification of drought duration and severity was defined, and is shown in Tables 1, 2.

The Distribution Characteristics of Three Variables
The proportion of affected farmland area (PAFA) in each year was used for analysis. The mathematical expression is: where I refers to PAFA, A d is the affected farmland area, and A P is the annual agricultural planting area. An exponential function is applied to describe the distribution of the index I, and the mathematical expression is defined as: Some studies [62,64] have reported that drought duration and severity follow exponential and Gamma distributions, respectively. The mathematical expressions of the two distributions are: where D and S refer to the sample set of drought duration and severity, respectively, d and s denote the element in the sample set of drought duration and severity, respectively, F D d refers to the probability that drought duration d is equal to or less than D, F S (s) denotes the probability that drought severity s is equal to or less than S, and λ, α, and β are the distribution parameters. According to Equations (2-4), the Copula function was used to establish the joint distribution function of the three variables (the duration and intensity of the drought, as well as the PAFA caused by drought) in this study.

Copula Function
The Copula function is an effective method that uses the marginal distribution functions of different random variables to build a joint distribution function. As each marginal distribution of the variables is known, the joint distribution can be constructed with the Copula function. According to the theory of Sklar [39], the joint distribution function F (x 1 , x 2 , · · · , x n ) can be decomposed into marginal distribution functions F (x 1 ), F (x 2 ), · · · , F (x n ) and  a Copula function C (·). The mathematical expression is: The Copula function can be classified into three classes: the elliptic type, Archimedean family, and quadric form. The Archimedean Copula function has been widely used [62,64], and therefore, it has also been used in our study. The common 3D Archimedean Copulas functions [40,65] are as follows. Clayton Frank Gumbel-Hougaard Where, µ, v, ω are the marginal distribution functions, and θ is the parameter of the Copula function.
In this study, the parameter θ is estimated through the maximum log-likelihood estimation. The function can be written as follows: where µ, v, ω are marginal distributions, and c(µ, v, ω, θ ) is defined as c (µ, v, ω, θ ) = ∂C(µ,v,ω,θ ) ∂µ∂v∂ω . The log-likelihood function ln L can be maximized in order to obtain the estimation of Copula parameterθ . Furthermore, the root-mean-square error method [66,67] was used to select the optimal Copula function: where Pe i is the empirical probability value obtained directly from the sample size, and P i is the theoretical probability value under different Copula functions. We compared and analyzed three Archimedes Copula functions, suggesting that the difference between the three Copula functions is not significant. However, the Gumbel-Hougaard function has an upper tail dependency, which is more suitable for the analysis of the dependence among three variables (the duration and intensity of a drought, as well as the PAFA caused by a drought). Thus, our subsequent analyses were all based on the Gumbel-Hougaard type Copula function.

Joint Distribution
In recent years, some researchers have begun to use threedimensional Copula functions in their research [68,69]. The joint distribution function can be obtained from the marginal distribution functions F D d , F S (s), and F I (i), as well as the Copula function of drought duration, severity, and PAFA.
The joint distribution describes the probability that the drought duration, severity, and PAFA are all equal to or less than a given value of the drought event.

Analysis Process
The main analysis process in this paper is as follows: first, the Run theory was used to extract the drought events by using the SPI sequences in Huize County and Mengzi County from 1961 to 2011. Then, the marginal distribution functions and their parameters were obtained using the data on drought durations and drought intensities, respectively. After that, we calculated the PAFA, and the marginal distribution functions and their

RESULTS AND DISCUSSION
The Establishment and Analysis of the 3D Joint Probability Distribution The SPI sequences of the Huize and Mengzi counties were applied to extract drought events during the period between 1961 and 2011. Statistical characteristics were obtained and shown in Table 3, with s 0 = 0 set as the threshold. The number of drought events extracted from the SPI sequence were 161 (Huize County) and 153 (Mengzi County) from 1961 to 2011. The statistical indicators of drought duration and severity in Mengzi County were slightly larger than in Huize County. The average drought duration of both counties was close to 2 months, and the average drought severity of both counties reached the level of severe drought. The longest drought duration in Mengzi County was 12 months, and its cumulative drought severity reached the maximum recorded value (12.12) in history. The Pearson correlation coefficient of drought duration and drought severity indicated that there was a strong correlation between the two variables in the two counties, showing that the duration and severity of the droughts in the two counties were both strongly synchronized (Table 3). Moreover, each pairwise relationship among drought duration, drought severity, and PAFA is measured using Kendall's tau coefficient (Table 4), the results indicating that there exists a positive interrelated relationship between any two of these three variables in the two counties because the values of Kendall's tau are all > 0.   Figure 4 shows the scatter diagrams and statistical histograms of the drought duration and severity in the two counties, respectively. The distribution characteristics of drought duration and severity in the two counties were basically consistent (Figure 4). In most cases, the drought durations were all under 2 months and the drought severities were lower than the level of a severe drought in the two counties.
In order to confirm that the distributions of the drought durations of the two counties satisfy the exponential distribution, drought severities of the two counties satisfy the gamma distribution, as well as having the distributions of the PAFA of the two counties satisfy the exponential distribution, the Kolmogorov-Smirnov (KS) test was used to analyze the difference between the empirical distribution function and the theoretical distribution function at the 0.05 significance level. The results according to the KS test are shown in Table 5. Table 5 shows that the results had passed the statistical test and were reliable; namely, the distributions of the drought durations of the two counties were subject to exponential distribution, the distributions of the drought severities of the two counties were subject to gamma distribution, and PAFA of the two counties were subject to exponential distribution.
In this study, Copulas are employed to construct the joint distribution function. A three-variable joint distribution function was obtained by using Equations (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11), and the probability distributions of PAFA were then calculated based on different drought durations and severities. Figure 5 shows the threevariable joint probability distribution of different PAFAs in Huize County, and the distribution characteristics were similar in Mengzi County.
The Z coordinate in Figure 5 represents the probability when the random variables are less than given values. Taking Figure 5A as an example, the meaning of the curved surface in the figure is the probability P D ≤ d, S ≤ s, PAFA ≤ 10% . The meaning in Figure 5B is the same as Figure 5A, except that the PAFA is less than a different value. As shown in Figure 5A, with regard to all drought events, 40% of the drought events caused PAFAs of no more than 10%; in other words, the farmland areas affected by 40% of the drought events did not exceed 10% of the planted areas. In addition, 80% of the drought events caused PAFAs of no more than 30%, meaning that the drought events, which caused PAFA to exceed 30% occurred rarely in history ( Figure 5B). Then, the 3D joint probability distribution and the statistical relationships between the PAFA, the drought duration, and the drought severity could be obtained.
In this paper, more attention was actually paid to the practical application of the 3D joint probability distribution rather than the statistical relationships between the PAFA, the drought duration, and the drought severity, so we propose a classification scheme. By combining the classification method of drought duration and severity (Tables 1, 2) with the classification standard of PAFA ( Table 6) (National Standard of the People's Republic of China: Grade of drought disaster GB/T 34306-2017), the probability distribution of various combinations of different drought durations, severities, and PAFAs could be determined, as shown in Figure 6. Figure 6 shows the joint probability distributions of Huize County that drought duration and drought severity are within a given interval. Figure 6A shows the occurrence probability (P d 1 < D ≤ d 2 , s 1 < S ≤ s 2 , 10% < PAFA ≤ 30% ) of drought events with different durations and severities when the PAFA is at the level of 2. The durations of drought events with a high probability of occurrence are at levels 1 (0 < D ≤ 1 within 1 month) and 2 (1 < D ≤ 3 within one season); that is, the durations do not exceed one season. If the drought duration is in level 1 (0 < D ≤ 1 within 1 month), the severity of the most likely drought events is in level 1 (0.5 < S ≤ 1 slight drought), but the severity of the most likely drought events is in level 4 (2.0 < S extreme drought) while the drought duration is in level 2 (1 < D ≤ 3 within one season). The probability of drought events with other durations and severities is very small. Similarly, Figures 6B-D gives the occurrence probability of drought events with different durations and severities when the PAFA is at three other levels. When the PAFA is at level 3 (30% < I ≤ 50%) (Figure 6B), the durations of drought events with a high probability of occurrence are at levels 2 (1 < D ≤ 3 within one season) and 3 (3 < D ≤ 6 cross-quarter), as well as severity level 4 (2.0 < S extreme drought); when the PAFA is at level 4 (50% < I ≤ 80%) (Figure 6C), the durations of drought events with a high probability of occurrence are at levels 3 (3 < D ≤ 6 cross-quarter) and 4 (6 < D over 6 months), as well as severity level 4 (2.0 < S extreme drought). The conclusion of Figure 6D is similar to Figure 6C but for the PAFA being at level 5 (80% < I). As the PAFA increases, the joint probability of drought duration and severity also increases with the increase in their levels.  Table 7 shows the disaster data during this period.

The Application of the 3D Joint Probability Distribution
There was one piece of recorded disaster drought data in Huize County in 2012 (during the period between January and March); at the same time, there were three records in Mengzi County in 2012 ( Table 7). By comparison, it was found that the start and end times of the drought disasters given by the recorded disaster data were basically consistent with the calculation results of the SPI index sequences in the two counties. The calculation result of the SPI index shows that the drought event in Huize County lasted for 2 months from January to February 2012 and the severity was 1.90, and the drought event in Mengzi County lasted for 5 months from February to June 2012 and the severity was 3.09, although the drought intensity in April-June was weak. The levels of drought duration and severity in Huize County were 2 (1 < D ≤ 3 within one season) and 3 (1.5 < S ≤ 2 Moderate drought), respectively, while the levels of drought duration and severity in Mengzi County were 3 (3 < D ≤ 6 cross-quarter) and 4 (2 < S Extreme drought), respectively. Equation (8) was used to estimate the occurrence probability of PAFAs that correspond to the two drought events in the two counties in 2012, respectively, and the results are shown in Table 8. Table 8 shows that in Huize County, the highest probability of an occurrence of a PAFA was 10% < I ≤ 30%; that is, the maximum risk corresponded to a PAFA of 10%; and in Mengzi County, the highest probability of an occurrence of PAFA was  Table 6. On the one hand, the underestimation of Mengzi County can be understood from a statistical point of view because we gave the occurrence probabilities of different levels of PAFA in a certain drought event, and while the probability is precisely an expression of risk, and on the other hand, the drought index used in this study (i.e., the SPI) is a meteorological drought index rather than a direct indicator of agricultural drought. Therefore, there is still a certain difference between the duration/severity of a drought obtained through the SPI index and the actual agricultural drought. However, a meteorological drought always occurs before an agricultural drought, so a meteorological drought index could be used to forecast and forewarn farmers of an agricultural drought and its effects in advance. In practical applications, compared with simple drought monitoring indicators, we can get SPI predictions for the following month or months based on monthly or seasonal precipitation forecasts, and can further obtain more abundant information about a drought warning and its impact through the method proposed in this paper so that relevant countermeasures can be taken according to the probability of occurrence of different levels of PAFA in order to achieve the purpose of disaster prevention, mitigation, and relief.

CONCLUSIONS AND DISCUSSIONS
Copula-based distribution has been used by researchers from different countries around the world in multivariate analyses of hydrological or meteorological events such as droughts due to its advantages in modeling the non-linear dependence structure of variables regardless of their marginal distributions. However, the applications of the Copula technique in assessing the impacts of a drought on agricultural production are very limited in published literature, and the damage data are often not used as a variable to be substituted into the Copula function, but are instead used as an independent variable in combination with the comprehensive features of a drought obtained by the Copula function. This assessment is actually based on an "external combination" of drought characteristics and drought impacts.
In order to achieve the "internal combination" of drought characteristics and drought impacts, as well as assessing the specific impacts of droughts with different severities and  durations, we used three variables of drought events and their impacts; namely, the drought duration, drought intensity, and the proportion of affected farmland areas (PAFA), and built a model based on the three-dimensional (3D) Copula function for the assessment of the PAFA for the purpose of quantified drought risk analysis. The result of the model can provide occurrence probabilities for different levels of PAFA in a certain drought event, and while probability is precisely an expression of risk, we took Huize County and Mengzi County in Yunnan Province as application examples, and the results showed that the PAFAs obtained by the method proposed in this paper were both consistent with the actual PAFAs in the two counties. Based on this model, the loss size and uncertainty (probability) under a given drought intensity can be well-expressed, which is the expected result of risk analysis. What's more, according to the division of drought grade, what kind of loss and its probability caused by a certain type of drought event can be clearly given. Moreover, in order to take risk control and mitigation measures, it is necessary to assess the impacts of a drought on agricultural production based on current drought conditions, and, more importantly, based on the forecast and predicted drought conditions in the future. A meteorological drought always occurs before an agricultural drought. In practical applications, we can obtain SPI predictions for the following month or months and can obtain further abundant information on drought warnings and their impact according to the probability of occurrence of different levels of PAFA. Therefore, the method proposed in this paper has values both on theory and practice.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
WH, DZ, and PY: methodology. WH and DZ: writing original draft preparation. GF and WH: writing review and editing. DZ and PY: visualization. All authors contributed to the article and approved the submitted version.