1 INTRODUCTION

The South China Sea is one of the three major marginal seas in China, lies on the coast of South China, its long and winding coastline makes it the most devastated area by typhoon and storm surge disasters (Yin et al., 2019; Zhang et al., 2019). According to statistics (Dong et al., 2014; Zhang, 2019), the number of typhoons in the South China Sea and the Northwest Pacific accounts for 36% of the total number of typhoons in the world, taking the first place among the eight typhoon generating regions in the world. During 1950–2020, there had been 319 typhoons storm surges landed in the northern coastal areas of the South China Sea, with an annual average of 4.5 in counts. Among them, a typhoon storm surge water increment of 5.94 m was recorded at Nandu Station in Haikang County, Guangdong Province (July 22, 1980), which is the maximum water increment ever recorded in Chinese history. According to the "China Marine Disaster Bulletin" statistics, in 2014, the South China Sea district had encountered the most serious loss caused by a storm surge disaster in the past 21 years. In that year, the direct economic loss of the district amounted to 125.29 million yuan, accounting for 92.27% of the total direct economic loss of the whole country in that year, and the disaster had caused 6 people died (including missing), as shown in Fig. 1. "Rammasun" (the strongest typhoon appears in China since its landing in 1949) and the landing of "Hato" result in the disaster in South China Sea district in 2014, the maximum water increment during the disaster was 3.92 m and 4.95 m respectively. In 2008, the storm surge disaster in 3 provinces (Guangdong, Guangxi, and Hainan of China) caused the death of (including missing) 55 people, with 17 355 million yuan direct economic loss, accounting for about 90% of the national storm surge disaster losses in that year; among them, Guangdong province was the most severely affected, accounting for about 80% of the total loss. What can be learnt from that is the storm surge disaster should not be underestimated. At the same time, the northern coastal areas in the south of China are also facing a serious threat from storm surge disaster. For that reason, establishing an assessment system that can accurately evaluate storm surge disasters and determine the risk level to provide a direct reference basis for coastal city planning, emergency evacuation, hazard insurance, etc. is crucial and urgent (Zeng et al., 2017; Liu et al., 2019b, 2020; Chen et al., 2021).

Storm surge disaster assessment is the basis of its risk management (Wang, 2014; Xu et al., 2018). The assessment of its risk level is mainly carried out in terms of disaster-causing intensity and major disaster losses. In recent years, the major method used for disaster loss assessment is the indicator evaluation method (Guo and Li, 2018; Zhao et al., 2020). It is mainly based on economic loss, number of casualties, the number of people killed and injured, and the usage of various regression models (e.g., neural network (Feng and Liu, 2017), vector machines (Wang and Liu, 2018), etc.) to construct an evaluation system. The evaluation mechanism of indicator evaluation method is relatively simple. However, since it is based on disaster data, it is more suitable for post-disaster assessment. Starting from the causes and results of disasters, evaluation methods based on numerical simulation (ADCIRC model (Gao et al., 2012; Guo et al., 2020)) and GIS technology (Li et al., 2017; Wang et al., 2017b) are gradually proposed. Assessment based on probability statistics is a traditional assessment method widely used in engineering (Zhao et al., 2012, 2013), which can effectively avoid the problem of regional differences, which cannot be avoided in the classification of disaster factor intensity or disaster loss. The designing of water increment under the influence of storm surge was calculated by the extreme value model, and designing of water increment was used to map the storm surge disaster level (Hisamatsu et al., 2020). Wang et al. (2021) used the Gumbel distribution model and the Gumbel-logistic model to calculate the return period of water increment and the joint water increment-wave height return period, individually, and determined the storm surge intensity level based on the return period level. However, the potential influence of typhoon frequency on the return period is also significant, and it is far enough to analyze the intensity of the disaster factor alone. Liu et al. (2021) constructed Poisson-Gumbel distribution and Poisson-Weibull distribution successively for long-term extreme value prediction of wave height in typhoon-affected waters, and obtained better designing values, which were widely promoted and applied (Wang et al., 2013, 2017a; Liu et al., 2021; Zhang et al., 2021). Mínguez et al. (2013) proposed that the Pareto-Poisson distribution model gives a more accurate estimate of the regression level. Thus, it is necessary to consider the effect of typhoon frequency in the extreme value model when assessing the storm surge risk level.

The basis for using the extreme value model to derive the designing water increment is to conduct data sampling. According to Liu et al. (2016), the longer period of the input data, the more stable the designing results of storm surge water increment. However, due to the short years of hydrological observation in China, it is difficult to collect data. As a result, the sample sequence obtained by the traditional sampling method (the annual extreme value method) is often short. Moreover, the actual observation often occurs in a year when the second largest value is greater than the annual maximum value. In this case, if the annual extreme value method is used for sampling, the second largest value of the year will be eliminated, which leads to the observation data in waste. In contrast, the Peak Over Threshold model (POT model) can obtain more sample points by setting reasonable thresholds, and make greater use of the limited measured information (Liu et al., 2019c; Xu et al., 2019; Deng et al., 2021). In recent years, POT model theory is also widely used in risk analysis (Song and Wang, 2020; Zhang and Duan, 2020), vehicle and machinery (Yuan et al., 2016; He et al., 2019) insurance as well as, and other fields (Jaladri et al., 2021). Therefore, sampling based on POT theory is expected to obtain more effective data.

In addition to sampling, whether an appropriate extreme value distribution is selected to fit the sample series would largely affects the results of probability analysis. It is generally believed that the over-threshold sample approximation obeys the Generalized Pareto Distribution (GPD) (Chen et al., 2019). Petrov et al. (2019) showed that the GPD distribution is much more sensitive to thresholds than the maximum entropy method, and higher thresholds improve the rationality of model expansion (Mazas and Hamm, 2011). The GPD distribution improves the fit of the "thick tail" of the over-threshold sample, but does not optimize the fit of both the peak and the right tail (Wang et al., 2016; Liu et al., 2019a). In contrast, the classical extreme value distribution can describe the peaks of the samples better, but the tails of its probability density function exhibit a power function form, and the fit to the right tail of the data is obviously insufficient (Muraleedharan et al., 2016; Mackay and Johanning, 2018). If the classical extreme value distribution can be effectively combined with the generalized Pareto distribution, it will be the key to improve the fitting accuracy. Aljarrah et al. (2014) proposed a method to generate a generalized cumulative distribution function of continuous distribution functions—the theory of family of distribution functions. Based on this theory, many scholars have constructed various combinatorial forms of this type of extreme value expansion. For example, Ganji et al. (2016) constructed the Weibull-Rayleigh distribution, and Alzaatreh et al.' s(2012, 2013) Gamma-Pareto distribution and Weibull-Pareto distribution.

Based on above analyses, this paper proposes to construct a typhoon storm surge disaster rating assessment system based on the POT theory and the composite extreme value expansion model, in order to improve and enhance the scientificity of the assessment model in terms of both data sampling and extreme value distribution. Specifics: firstly, a new sample series is obtained by sampling in the POT method. Then, with the method of generating a family of distribution functions proposed by Aljarrah et al. (2014), the Weibull and Pareto distributions are effectively combined to form an extreme value expansion model to improve the fitting of the thick tail of the sample distribution. Then, the Poisson-Weibull-Pareto composite extreme value expansion model is constructed based on the composite extreme value theory by considering the influence of typhoon frequency on the level of return period. The design value of the multi-year water increment calculated by the new model is also used as the standard for typhoon storm surge level zoning. Finally, the zoning results are compared with those evaluated based on the criteria of the SOA (2012) for the 10 typhoon storm surges that made landfall in the Guangdong area during 2010–2020.

2 CONSTRUCTION OF RATING EVALUATION SYSTEM

The construction of typhoon storm surge disaster level assessment system is mainly divided into three stages: sampling stage, calculation of designing water increment stage, and level zoning stage. The sampling process mainly conduct based on the peak over threshold theory, which is extracted from the original data series after setting a reasonable threshold to obtain a new sample series. The second stage is taking advantage of a composite extreme value expansion model constructed to calculate the multi-year designing water increment. The last stage is the delineation partition of the risk level based on the values from different return periods.

2.1 Threshold sampling criteria determination

Peak over threshold (POT) is a sampling method that selects a fixed larger value as the "threshold" through a certain method or technical means, and selects the points that meet or exceed this "threshold" to form a new sample series. In order to ensure a scientific and reasonable analysis sample, this paper adopts a relatively objective and quantitative sampling method. That is: based on the mean excess function of Pareto distribution, draw the average remaining life of the water increment series, and then determine the threshold value by observing the starting and ending position of the straight line segment on the graph.

If the water increment X(x>0) obeys the Pareto distribution, its cumulative density function is

Its probability density function is

where α and σ are the shape and scale parameter, respectively.

At this point, the mean excess function e(u) of the Pareto distribution can be determined as

That is, there is a linear correlation between e(u) and u. where u is the threshold value of u∈D(α, σ), σ < 1. E(x–u|x>u) is the expectation of over the threshold excess. According to the above equation, a scatter plot of the mean value of the threshold u and the observed excess (x–u|x>u) (called the mean excess) can be made. When the shape parameter σ is stable, the graph is approximated as a straight line. That is, when the threshold u is the horizontal axis and the mean excess is the vertical axis, the slope and intercept of the line are and , from it here, the horizontal coordinate corresponding to the straight line segment in the average remaining life graph can be used as the optional range of the threshold value.

In particular, the average excess function e_n(u) for samples x₁, x₂, ···, x_n can then be expressed as

At this point, the average remaining life can be plotted according to {(u, e_n(u)), u < x_n}, the linear variation of e_n(u) after the threshold u₀ can be analyzed, and if the slope of e_n(u) can be approximately constant, then u₀ is considered as the optimal threshold. In addition, the scale parameter α and shape parameter σ of the Pareto distribution can be estimated using the maximum likelihood estimation method (MLE) method within a reasonable threshold range. When the parameters α and σ are stable, u₀ can be considered the optimal threshold.

2.2 Composite extreme value expansion model—threshold embodied extreme value expansion model construction

The construction of the extreme value expansion model is based on a family of distribution functions as the underlying theory. The method of construction is to combine two extreme value distribution functions so that they add additional parameters to the existing distribution. In this paper, the statistical properties of the random variable X are discussed in terms of the auxiliary variable T. The family of continuous function probability distributions is constructed as follows.

Theorem 1 Let the random variable T be a continuous random variable defined on the interval [a, b] with r(t) as the probability density function, -∞ ≤a < b≤+∞, x be a random variable with h(x) as the probability density function and the distribution function H(x), and there exists a distribution function G(x)

This family of probability function distributions is called the extreme expansion mode, and is denoted as the extreme extension distribution (EED) mode.

Specifically, if the marine environment is to be T defined on the interval [a, b] and obeys the generalized extreme value distribution (GED) R(t), -∞≤a < b≤+∞, with r(t) being the probability density function of the random variable T.

in which, ξ is the shape parameter; σ₁ is the scale parameter.

The random variable X obeys a Pareto distribution with a distribution function of H(x) and a probability density function h(x) = ασ^αx^-α–1, α > 0, x≥σ > 0, where α is the shape parameter and σ is the scale parameter and also the threshold value. Then the probability density function g(x) corresponding to the distribution function G(x) can be expressed as

Here, α represents the shape parameter of Pareto distribution, and ξ is the shape parameter in generalized extreme value distribution.

Let α/σ₁=β. The above equation can be simplified

Name the Eq.8 as the GED-Pareto distribution, with the notation as EED (β, σ, ξ).

In fact, the generalized extreme value distribution is a collective term for three typical extreme value distributions, which are the Gumbel distribution, the Weibull distribution and the Frechet distribution. Distributions of the maximum values of independent random variables basically converge to these three types of extreme value distributions.

The Weibull distribution is also known as the extreme type Ⅲ distribution, when

then the probability density function

Let α/γ=β, then the simplification gives

Name the Eq.11 as the Weibull-Pareto distribution with the notation EED-Ⅰ type (c, β, σ).

The probability density of EED-Ⅰ with different parameters is shown in Fig. 2.

In particular, it has the following two special cases.

(1) When σ=1, Eq.11 becomes:

At this point, it becomes a logarithmic Weibull distribution.

(2) When c=1, Eq.11 becomes:

At this point, the Weibull-Pareto distribution also degenerates to the Pareto distribution.

The design parameters of various marine environmental elements in the sea affected by typhoons and storm surges are usually calculated in the way of composite extreme value distribution. It is a distribution pattern formed by the discrete distribution of typhoon frequency and the continuous distribution of various marine environmental elements under the influence of typhoons.

If n denotes the frequency of typhoons and obeys the Poisson distribution, then

where G(x) in the above equation can be any kind of one-dimensional continuous distribution.

If the newly constructed extreme expansion model is calculated by substituting into Eq.14, the composite extreme expansion model Poisson-EED-Ⅰ type is obtained.

2.3 Storm surge risk classification

As the storm surge risk classification criteria, China's State Oceanic Administration (SOA) has given a classification scheme based on storm surge in the technical guidelines for its disaster assessment and zoning in 2012 (SOA, 2012) (Table 1). Since there are significant differences in storm surge disasters in different regions, the size of storm surge is affected by various factors, especially by microtopography, and storm surge may vary by tens or even hundreds of centimeters in different river sections in the same estuary. Therefore, it is not possible to avoid regional differences by directly using the intensity of the causative factor (storm surge) as an indicator for risk classification. The frequency of occurrence is a measure of the likelihood of disaster occurrence, and the frequency of disaster occurrence is an important indicator of its severity, which can objectively reflect the regional disaster risk and avoid the above problems. In this paper, a typhoon storm surge classification scheme based on frequency of occurrence is given (Table 1), and different levels of annual occurrence are used to portray the disaster level of different regions, the higher the frequency, the shorter the return period, the higher the level, and the lighter the disaster level.

3 ENGINEERING APPLICATION AND CALCULATION

Typhoon and storm surge disasters happen frequently in the South China Sea region of China. Guangdong Province, bordering the South China Sea, is one of the hardest hit areas in China suffering from the devastating threat of typhoon and storm surge disasters. Among them, typhoon and storm surge disasters occur most frequently in western Guangdong and Zhujiang port area, with an average of three times per year. Therefore, in this paper, the western part of Guangdong is selected as the study area, and the data of water increment and typhoon statistics of Naozhou Hydrological Observatory for a total of 25 years from 1990 to 2016 (missing 2004 and 2007) are selected as the original data for the study. The statistics of the frequency of typhoon disasters in the western part of Guangdong are shown in Table 2; the frequency of typhoons is consistent with the Poisson distribution by χ² test, and the parameter λ of Poisson distribution is 3.24 (the Naozhou Hydrological Station is located on the northeast coast of Naozhou Island, with the coordinates 110.550 0°E and 20.900 0°N).

3.1 Sample analysis and model comparison

In this paper, the average remaining life of over-threshold water increment samples was plotted based on the mean excess function of the Pareto distribution, and the values of the thresholds were analyzed and determined in combination with the shape and scale parameter stability plots (Fig. 3).

As shown above, under 95% confidence interval, the water increment threshold can be selected in range of 0.1–0.5 m according to the linear trend of the mean excess in the mean remaining life graph relative to the threshold. On this occasion, the determination of the specific value of the threshold still has a relatively large subjectivity. For this reason, the parameter stability assessment was introduced into the mean remaining life diagram method, and the stability was finally determined to be 0.24 m by combining the modified scale parameter (modified scale) and the maximum likelihood estimation of the shape parameter (shape), under the premise of ensuring the minimum sampling error of the POT method.

Now the extreme water increment data of Naozhou station for a total of 25 years are divided into two groups based on two different methods. Group A is the annual extreme water increment series from 1990 to 2016; Group B is the water increment over threshold series from 1990 to 2016 with a threshold value of 0.24 m, as shown in Figs. 4–5.

Comparing Fig. 4 with Fig. 5, it can be referred that the series length of the over-threshold sample is 140% longer than that the annual extreme sample series. The over-threshold sample contains more actual measured point data, allowing more valid observations to reflect useful information in the subsequent analysis.

The results of the test for the water increment series A and B using the K-S test are shown in Table 3.

From the test results in Table 3, it is clear that both sets of water increment series passed the K-S test.

For further analysis, the following test plots of the extreme value distribution of the A and B water increment series samples were made (Figs. 6–7). The solid lines in the plots indicate the model curves and the circles indicate the data points.

Analysis of the test plots of the extreme value distribution of the water-increase samples of the A and B series shows that most of the data points in the Quantile-Quantile plot fall on the diagonal line, and some of them are more evenly scattered on both sides of the diagonal line, which shows that the observed data points and the model curve are in good agreement. From the return level plot, it can be seen that the data points are completely distributed within the 95% confidence interval of the return level of the model.

In summary, the analysis shows that the selected water increment data sets are in good agreement with each extreme value distribution model and can be used as samples for the marginal analysis of the new distribution.

In general, the skewness coefficient S and kurtosis coefficient K can be used to describe the degree of symmetry and the tail shape of the random sample distribution, respectively. Specifically, they can be defined as the third and fourth central moment of the sample. If the sample has a positive S value, the sample distribution is right-skewed, while the opposite is left-skewed; if the K value of the sample is greater than 3.0, the probability distribution of the sample is "spiky and thick-tailed", which means that more extreme values of the random sample are located in the tails of the sample distribution. Table 4 shows the kurtosis coefficient K and skewness coefficient S for the water increment samples.

The kurtosis value greater than 3 in Table 4 is for the over-threshold sample, while the kurtosis value of the annual extreme sample is less than 3, indicating that the over-threshold sample has the characteristic of "sharp peak and thick tail". In terms of skewness: both samples are positively skewed (right skewed), and the skewness value of the over-threshold sample is much larger than that of the annual extreme value sample. This also indicates that the right skewness of the over-threshold sample is more pronounced.

To further analyze the samples and compare the distribution models, the sample histograms and the probability density curves of the fitted samples for each distribution are plotted in the same figure (Fig. 8). The left and right plots in Fig. 8 show the distribution test plots for the A (annual extremes) water increment series and B (over-threshold) water increment series, respectively. From the test results: compared to the left panel, the peak form of the probability density distribution curve in the right panel is still present and the peaks are larger overall. The right tail of the probability density curves of Gumbel distribution and Weibull distribution in the right figure shows the form of power function, which quickly converges to a small probability; while the right tail of the probability density function curve of EED-Ⅰ type has a longer and slower convergence fit. This demonstrates that the Weibull-Pareto distribution is more conducive to fitting the "thick tail" due to the characteristics of "peaks and thick tails" become more pronounced after sampling with the threshold method.

3.2 Designing water increment calculation and storm surge risk level classification

The POT method can extract as much effective information as possible from the observed data, but only when the sampling method is chosen properly to the corresponding statistical model that one can have a good fitting effect. Since there are few distribution models applicable to over-threshold samples and the computational accuracy of the existing models needs further engineering verification, the following is to verify the reasonableness and applicability of the over-threshold Weibull-Pareto model, which reference to the computed results of Gumbel, Weibull, and Weibull-Pareto distributions based on annual extreme data. The results of the goodness-of-fit tests for each distribution type under each sampling series are shown in Table 5. The results of the designing water increment calculation are shown in Table 6.

The goodness-of-fit test results in Table 5 show that each extreme value distribution passed the A-D test. According to the criterion of "the optimal distribution corresponds to the smaller A-D value", the Weibull-Pareto distribution is in good agreement with the distributions of the annual extreme water increment sample and the over-threshold water increment sample, respectively.

The Poisson-Weibull-Pareto model exhibits a significant "thick tail" in fitting the annual extreme water increment, i.e., the result of new model is significantly larger than the other two extreme models when the return period is greater than 50 years. Comparing the results of A and B series, we found that the difference between the Poisson-EDD-Ⅰ distribution based on B sample series (POT) and A sample series annual maximum (AM) was not significant, with the difference ranging from 0.72% to 12.68%. In contrast, the variation in the proportion of the difference between the calculated for Poisson-Gumbel distribution is 26.64%–28.03% and Poisson-Weibull distribution is 20.57%–25.09%, successively, both with relatively large differences. This suggests that it is more appropriate to use the Weibull-Pareto distribution to analyze the over-threshold sample series. Meanwhile, the Poisson-EDD-Ⅰ distribution has more stable results based on both samples. This also indirectly reveals that the extreme value expansion model has stronger applicability than the classical extreme value distribution when facing the variation of the sample series.

In summary, this paper will select the multi-year water increment design value based on the Poisson-EDD-Ⅰ distribution model as the criteria for typhoon storm surge risk classification, as in Table 7.

The water increment extremes during the landfall of typhoon storm surges that made landfall in the Guangdong area between 2010 and 2020 were selected (see Table 8). Based on the above classification criteria, the risk classification of the 10 typhoon storm surges will be made. The results of the zoning are shown in Table 8 and Fig. 9.

Although there are some differences in the designing water increment obtained based on different sample series, the results of the typhoon storm surge risk classification are generally consistent when the risk level is divided into zones. Figure 9 and Table 8 show that the peak-value of surge of the five typhoon storm surges, including Chandu, were all less than the 10-year recurrence level, with a risk level of grade Ⅳ. And according to the classification criteria of the SOA (2012), their disaster level is Ⅲ. In addition, according to the division criteria of the SOA (2012), the extreme value of water increment less than 1.5 m is disaster level Ⅳ. And the water increment value of 1.5 m is less than the design value of 5-year increase based on the Poisson-EED-Ⅰ extrapolation model. This shows that the classification criteria of SOA (2012) are more conservative and safer when classifying low-risk areas. The classification results of No. 1409 Rammasun and No. 1415 Kalmaegi, based on the return period, are roughly the same as those in SOA (2012), which are Ⅱ and Ⅰ, individually. In another word: in the high-risk areas, the assessment grade results based on the return period are more consistent with those of the State Oceanic Administration. Among them, the extreme value of water increment in Kalmaegi exceeds the 200-year return level of water increment.

4 CONCLUSION

Based on the transcendental threshold theory, combined extreme value distribution theory, and compound extreme value theory, the typhoon storm surge rating evaluation system proposed in this paper is suitable for determining the intensity of the typhoon storm surge in the west Guangdong sea area and it has certain reference significance for moisture-proof mitigation in other coastal areas in China. Here are the main conclusions drawn from this study:

(1) The use of POT sampling can make fully and more effectively use of the limited measured information to obtain more extreme sample point data and data analysis information, and avoid the influence of the annual cycle on the inter-annual variation of the extreme sample.

(2) Over-threshold samples tend to have "peaks and thick tails". Compared with the classical extreme value distribution, the Weibull-Pareto distribution can better fit such samples. When different sampling methods are used and the length of the sample series varies, the extreme value expansion model is more applicable than the classical extreme value distribution, and the imputed design values are more stable.

(3) Example analysis on typhoon storm surge classification shows that the classification criteria of the SOA (2012) are more conservative when classifying low-risk areas. In contrast, in the high-risk area, the results of the return period-based classification assessment are more consistent with those of the State Oceanic Administration.

5 DATA AVAILABILITY STATEMENT

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.