Joint Risk of Rainfall and Storm Surges during Typhoons in a Coastal City of Haidian Island, China

Public health risks from urban floods are a global concern. A typhoon is a devastating natural hazard that is often accompanied by heavy rainfall and high storm surges and causes serious floods in coastal cities. Affected by the same meteorological systems, typhoons, rainfall, and storm surges are three variables with significant correlations. In the study, the joint risk of rainfall and storm surges during typhoons was investigated based on principal component analysis, copula-based probability analysis, urban flood inundation model, and flood risk model methods. First, a typhoon was characterized by principal component analysis, integrating the maximum sustained wind (MSW), center pressure, and distance between the typhoon center and the study area. Following this, the Gumbel copula was selected as the best-fit copula function for the joint probability distribution of typhoons, rainfall, and storm surges. Finally, the impact of typhoons on the joint risk of rainfall and storm surges was investigated. The results indicate the following: (1) Typhoons can be well quantified by the principal component analysis method. (2) Ignoring the dependence between these flood drivers can inappropriately underestimate the flood risk in coastal regions. (3) The co-occurrence probability of rainfall and storm surges increases by at least 200% during typhoons. Therefore, coastal urban flood management should pay more attention to the joint impact of rainfall and storm surges on flood risk when a typhoon has occurred. (4) The expected annual damage is 0.82 million dollars when there is no typhoon, and it rises to 3.27 million dollars when typhoons have occurred. This indicates that typhoons greatly increase the flood risk in coastal zones. The obtained results may provide a scientific basis for urban flood risk assessment and management in the study area.


Introduction
Typhoons are considered extremely devastating natural hazards worldwide, which have caused enormous property and human losses and severely impacted public health [1]. Under the influence of global warming and the rise in sea level [2], the frequency and intensity of natural disasters such as typhoons, rainstorms, and storm surges have increased. Furthermore, there has been an increase in the damage caused by these disasters. Affected by the same meteorological systems, typhoons, rainfall, and storm surges are three correlated variables. Based on previous reports [3][4][5][6], although the correlation between rainfall and storm surges is often weak, it has a significant impact on coastal urban flood management. Strong typhoons often bring on heavy rainfall and high surges. Floods can easily occur in coastal cities when typhoons occur as they often bring heavy rainfall and storm surges. For example, Typhoon Rammasun, with a maximum wind speed of 60 m/s, attacked Haikou in Hainan Province, China, from 17 July to 19 July 2014, resulting in heavy daily rainfall (509.2 mm) and high storm surges (3.83 m) on 18 July. The return period of both rainfall and storm surges is In this study, daily rainfall data and daily maximum storm surge data during 1974-2012 were collected from the Haikou hydrological station, which are provided by the Haikou Municipal Water Authority. The continuity of the data was checked. The typhoon data are available from the best-track dataset by the Shanghai Typhoon Institute of China Meteorological Administration (CMA). The dataset contains information on each typhoon track every 6 h, including the time, location (latitude and longitude), maximum sustained wind (MSW), and minimum pressure near the typhoon center. The CMA typhoon best-track dataset is included in the International Best Track Archive for Climate Stewardship (IBTrACS) project [38], which is an official World Meteorological Organization (WMO) global archiving and distribution resource for typhoon best-track data [39].
When the distance between the typhoon center and the study area is less than 500 km [40,41], the typhoon is considered to have potential impact on the study area. There were 128 typhoon events affecting Haidian Island from 1974 to 2012, and the tracks of those typhoon events are presented in Figure 2. As shown in the figure, most of the typhoons that affect Hainan Island originate in the South China Sea and approach the island from east to west. The data of each typhoon point include MSW, center pressure, and distance between the typhoon center and the study area. They were selected by the following steps. First, we selected the affected daily rainfall and storm surges by typhoons through the manual identification method. For example, on 1 September, daily rainfall and storm surges were affected by typhoons. Then, we determined the corresponding typhoon points (usually more than one) on 1 September and calculated the distance between typhoon points and Haikou. Finally, the typhoon point with the smallest distance was selected as the basic typhoon data on 1 September, since rainfall and storm surges are mainly affected by the typhoon position. As shown in Figure 3, typhoon points 3-5 affected Haikou on 1 September, and typhoon point 4 had the smallest distance. Therefore, the MSW, center pressure, and typhoon position of typhoon point 4 were selected as the basic typhoon data. In this study, daily rainfall data and daily maximum storm surge data during 1974-2012 were collected from the Haikou hydrological station, which are provided by the Haikou Municipal Water Authority. The continuity of the data was checked. The typhoon data are available from the best-track dataset by the Shanghai Typhoon Institute of China Meteorological Administration (CMA). The dataset contains information on each typhoon track every 6 h, including the time, location (latitude and longitude), maximum sustained wind (MSW), and minimum pressure near the typhoon center. The CMA typhoon best-track dataset is included in the International Best Track Archive for Climate Stewardship (IBTrACS) project [38], which is an official World Meteorological Organization (WMO) global archiving and distribution resource for typhoon best-track data [39].
When the distance between the typhoon center and the study area is less than 500 km [40,41], the typhoon is considered to have potential impact on the study area. There were 128 typhoon events affecting Haidian Island from 1974 to 2012, and the tracks of those typhoon events are presented in Figure 2. As shown in the figure, most of the typhoons that affect Hainan Island originate in the South China Sea and approach the island from east to west. The data of each typhoon point include MSW, center pressure, and distance between the typhoon center and the study area. They were selected by the following steps. First, we selected the affected daily rainfall and storm surges by typhoons through the manual identification method. For example, on 1 September, daily rainfall and storm surges were affected by typhoons. Then, we determined the corresponding typhoon points (usually more than one) on 1 September and calculated the distance between typhoon points and Haikou. Finally, the typhoon point with the smallest distance was selected as the basic typhoon data on 1 September, since rainfall and storm surges are mainly affected by the typhoon position. As shown in Figure 3, typhoon points 3-5 affected Haikou on 1 September, and typhoon point 4 had the smallest distance.

Methods
The overall framework of the study is shown in Figure 4. First, to quantify typhoons, principal component analysis was used to characterize typhoons by integrating MSW, center pressure, and distance between the typhoon center and the study area. Then, a multivariate joint probability distribution of typhoons, rainfall, and storm surges was constructed by trivariate copula functions. Third, the impact of typhoons on flood risk caused by rainfall and storm surges was evaluated based on an urban flood inundation model and a flood risk model. Finally, the joint risk of rainfall and storm surges during typhoons was estimated from two aspects: (1) the impact of typhoons on joint probability of rainfall and storm surges, and (2) the impact of typhoons on flood risk caused by rainfall and storm surges.

Methods
The overall framework of the study is shown in Figure 4. First, to quantify typhoons, principal component analysis was used to characterize typhoons by integrating MSW, center pressure, and distance between the typhoon center and the study area. Then, a multivariate joint probability distribution of typhoons, rainfall, and storm surges was constructed by trivariate copula functions. Third, the impact of typhoons on flood risk caused by rainfall and storm surges was evaluated based on an urban flood inundation model and a flood risk model. Finally, the joint risk of rainfall and storm surges during typhoons was estimated from two aspects: (1) the impact of typhoons on joint probability of rainfall and storm surges, and (2) the impact of typhoons on flood risk caused by rainfall and storm surges.

Principal Component Analysis
The existing literature has focused on the dependence between typhoons, rainfall, and storm surges used to characterize typhoons by MSW [9][10][11]. However, typhoons have many indices, such as MSW, center pressure, distance between typhoon points and the study area, SiR34 (the radius, in km, of 34 kt, i.e., tropical cyclone (TC) size), and so on. MSW cannot fully reflect the impact of typhoons on rainfall and storm surges. In this work, to analyze the multivariate joint probability distribution of typhoons, rainfall, and storm surges, typhoons were quantified by a synthetic parameter integrating MSW, center pressure, and distance between the typhoon center and the study area, since these variables all have significant dependence with rainfall and storm surges ( Table 1). As shown in Table 1, since SiR34 has a lower correlation with rainfall and storm surges, it was not adopted in the base parameter set to perform the principal component analysis. Therefore, in this study, principal component analysis was used to characterize typhoons by integrating MSW, center pressure, and distance between the typhoon center and the study area.

Principal Component Analysis
The existing literature has focused on the dependence between typhoons, rainfall, and storm surges used to characterize typhoons by MSW [9][10][11]. However, typhoons have many indices, such as MSW, center pressure, distance between typhoon points and the study area, SiR34 (the radius, in km, of 34 kt, i.e., tropical cyclone (TC) size), and so on. MSW cannot fully reflect the impact of typhoons on rainfall and storm surges. In this work, to analyze the multivariate joint probability distribution of typhoons, rainfall, and storm surges, typhoons were quantified by a synthetic parameter integrating MSW, center pressure, and distance between the typhoon center and the study area, since these variables all have significant dependence with rainfall and storm surges ( Table 1). As shown in Table 1, since SiR34 has a lower correlation with rainfall and storm surges, it was not adopted in the base parameter set to perform the principal component analysis. Therefore, in this study, principal component analysis was used to characterize typhoons by integrating MSW, center pressure, and distance between the typhoon center and the study area. Principal component analysis is one of the most widely applied tools for simplifying a dataset by reducing multidimensional datasets to a smaller number of dimensions than the original representation without changing the characteristics of the original data [42]. For an n-dimensional observation sample matrix S with p eigenvalues, the principal component can be calculated with the following steps: Step 1. Establish observation sample matrix S: Step 2. Convert sample matrix S to normalized matrix Y: 1. Normalized matrix for positive indices: 2. Normalized matrix for negative indices: Step 3. Calculate the correlation coefficient matrix of the normalized matrix: Step 4. Calculate eigenvalues λ of correlation coefficient matrix R: where I is the p by p identity matrix, |·| is the determinant operator, and λ represents the p eigenvalues ranked in decreasing order, λ 1 ≥ λ 2 ≥ · · · ≥ λ p ≥ 0.
When the variance cumulative contribution rate of m principal components to 1 (generally greater than 85%), the factor variables U 1 , U 2 , · · · , U m are known as the first, second, ..., mth principal components of the original variables, respectively, which are expressed as: where y is the normalized matrix mentioned in Equations (2) and (3), and b j is the eigenvector of eigenvalue λ j . It can be calculated by the following formula.
Step 6. Characterize samples by integrating principal components: The sample values can be integrated by the following formulas [43][44][45].
where ω j is the weighting coefficient of principal component U j , and λ j is the variance contribution rate of principal component U j .

Copula Function
F 1 , F 2 , · · · , F n are marginal distributions of X 1 , X 2 , · · · , X n , respectively. According to Sklar's theorem [46], if F 1 , F 2 , · · · , F n are continuous, there exists an n-copula C as follows: In this study, two elliptical copulas (Gaussian copula and Student's t copula) and three Archimedean copulas (Gumbel copula, Clayton copula, and Frank copula) were employed to build joint distributions of typhoon-rainfall, typhoon-storm surge, and rainfall-storm surge. Trivariate Gaussian copula, Gumbel copula, and Student's t copula were employed to construct the joint distribution of typhoons, rainfall, and storm surges. The above copula functions are presented in Appendix A.
The joint probability, co-occurrence probability, and conditional probability of typhoons, rainfall, and storm surges can be calculated by copulas. For T representing the typhoon, H representing rainfall, and Z representing the storm surge, the joint probability of at least one variable (typhoon T, rainfall H, or storm surge Z) exceeding its extreme values is denoted as P ∪ (t, h, z). The expression of P ∪ (t, h, z) is as follows [47,48]: The co-occurrence probability of T, H, and Z all exceeding certain extremes is denoted as P ∩ (t, h, z). The expression of P ∩ (t, h, z) is as follows [47,48]: The conditional probability that H and Z exceed a certain extreme when T has exceeded the extreme is denoted as P(H > h, Z > z|T > t) . The expression of P(H > h, Z > z|T > t) is as follows:

Expected Annual Damage Evaluation
In this study, flood risk is defined as the product of probability and damages [49]. The expected annual damage (EAD) is then found to express the flood risk, which can be calculated by integrating the flood damage function with the probability function. An approximation for calculating EAD is described in [31]: where D P is the damage caused by a flood of exceedance probability P, m is the number of probability increments, and D i is the average flood damage (mean of D pi and D pi−1 ) during probability increment ∆P i for the ith interval. Flood damages are calculated by unit flood damage and maximum inundation depth. The inundation depth is calculated by PCSWMM, which is introduced in Appendix B.

Quantification of Typhoons
In order to establish the multivariate joint probability of typhoons, rainfall, and storm surges, typhoons should be quantified first. To quantify typhoons, characterization by principal component analysis was performed, integrating MSW, center pressure, and distance between the typhoon center and the study area. The results reveal that, on the basis of eigenvector loadings, the first principal component (PC1) with an eigenvalue of 0.093 is able to explain 59.1% of the total variation, whereas the second principal component (PC2) with an eigenvalue of 0.062 explains 39.3% of the variation, and both PC1 and PC2 explain 98.4% of the total variation (Table 2). Therefore, the first two principal components are sufficient to replace the original data information.
From the score coefficient matrix in Table 3, the score coefficients of MSW, center pressure, and distance in PC1 are 0.745, 0.589, and 0.312, respectively. MSW and center pressure have higher score coefficients than distance, indicating that PC1 is the comprehensive reaction of MSW and center pressure. The second PC, negatively loaded with MSW (−0.247) and center pressure (−0.19), is positively loaded with the distance between the typhoon center and the study area (0.95). The proportion of distance is the highest in PC2, which indicates that PC2 is a description of distance. Therefore, PC1 and PC2 can fully reflect MSW, center pressure, and the distance, indicating that the principal component analysis results are reasonable. Furthermore, the principal components can be described as follows.
where S 1 is MSW, S 2 is center pressure, and S 3 denotes the distance between the typhoon center and the study area. From Equations (8) and (9), typhoon T can be quantified by the following equation [43][44][45]: The quantitative correlation among variables was analyzed using Pearson's correlation coefficient r and two nonparametric dependence measures, Kendall's τ and Spearman's ρ. Table 4 presents corresponding correlation measures between typhoons, rainfall, and storm surges. The correlation between variables was found to be statistically significant at 1% significance level, as checked by a standard two-tailed t-test. The correlation is largest between typhoons and rainfall (Spearman's ρ = 0.368 at 1% significance level), while it is the smallest between rainfall and storm surges (Kendall's τ = 0.142 at 1% significance level). The reason may be that rainfall and storm surges are affected by various weather systems, and typhoon is only one factor. Additionally, the angle of approach of a typhoon may also cause the low correlation between rainfall and storm surges (e.g., being on the left side of the typhoon can reduce the storm surge, while it produces maximum rainfall on the left side). Even though the correlation between variables is low, it can have significant implications for flood risk estimates, and there are still many studies about their joint probability analysis [3][4][5][6][7][8]. Furthermore, to analyze the correlation between extreme events, upper-tail correlation coefficient λ u was also evaluated, and the introduction of λ u is described in the literature in detail [50]. As shown in Table 4, upper-tail correlation coefficients are higher than Pearson's correlation coefficients, Kendall's τ, and Spearman's ρ among typhoons, rainfall, and storm surges, indicating that the correlation among extreme events is stronger.

Trivariate Joint Probability Distribution of Typhoon-Rainfall-Storm Surge
The nonparametric kernel density estimation [51,52] was used to establish the marginal distribution of typhoons, rainfall, and storm surges. As shown in Table 5, all of the computed values of the Kolmogorov-Smirnov (K-S) statistic D are lower than the critical values (D 0.01 = 0.067). Furthermore, a comparison of kernel density estimations and empirical distributions of typhoons, rainfall, and storm surges are presented in Figures 5-7, indicating that the kernel density estimation can properly estimate the distribution functions of typhoons, rainfall, and storm surges. As for bivariate joint distributions, typhoon-rainfall, typhoon-storm surge, and rainfall-storm surge are all best fitted by the Gumbel copula due to minimal Akaike information criterion (AIC) statistics being found for these bivariate distributions (Table 6). Table 7 shows the results of goodness of fit of the trivariate distribution of typhoon-rainfall-storm surge. It can be seen that typhoon-rainfall-storm surge is best fitted by the Gumbel copula, as it has the lowest AIC statistics and passes the K-S test. Figure 8 illustrates the probability-probability (P-P) plot of the joint distribution of typhoon-rainfall-storm surge. The coefficient of determination between empirical distribution and copula distribution is above 0.99, which indicates that the selected copula distribution is reasonable, and the selected parameters are adoptable. The trivariate joint probability distribution of typhoons, rainfall, and storm surges is expressed in Figure 9. With an increase in typhoons, rainfall, and storm surges, their joint probability distribution increases.  The values of the distributions in bold were selected. The values of the distributions in bold were selected.
K-S statistic D 0.017 0.066 0.011 Table 6. Fitting results of bivariate distribution functions. AIC, Akaike information criterion. The values of the distributions in bold were selected. The values of the distributions in bold were selected.

Joint Probability of Typhoon-Rainfall-Storm Surge Analysis
The joint probability and co-occurrence probability of typhoons, rainfall, and storm surges are calculated from Equations (11) and (12). From Table 8, we can conclude the following: (1) With a decrease in the return period (RP), the joint probability of typhoons, rainfall, and storm surges increases. (2) The trivariate joint probability is always greater than the co-occurrence probability in all conditions. Hence, the simultaneous occurrence of typhoons, rainfall, and storm surges all exceeding certain threshold values is less frequent. However, one variable (typhoon, rainfall, or storm surge) more frequently exceeds its threshold value. (3) The joint probability P is nearly two  Figure 9. Some bivariate views of the trivariate probability distribution of typhoons, rainfall, and storm surges. Joint probability distribution of rainfall and storm surges when typhoon value is (a) 0.

Joint Probability of Typhoon-Rainfall-Storm Surge Analysis
The joint probability and co-occurrence probability of typhoons, rainfall, and storm surges are calculated from Equations (11) and (12). From Table 8, we can conclude the following: (1) With a decrease in the return period (RP), the joint probability of typhoons, rainfall, and storm surges increases.
(2) The trivariate joint probability is always greater than the co-occurrence probability in all conditions. Hence, the simultaneous occurrence of typhoons, rainfall, and storm surges all exceeding certain threshold values is less frequent. However, one variable (typhoon, rainfall, or storm surge) more frequently exceeds its threshold value. (3) The joint probability P ∪ is nearly two times the univariate occurrence probability of typhoons, rainfall, and storm surges, and flooding would easily occur when either of them exceeds the threshold. However, the flood design standard in China is only determined by univariate analysis (i.e., rainfall frequency analysis), which would highly underestimate the flood risk in coastal regions. Table 8. Joint probability (P ∪ ) and co-occurrence probability (P ∩ ) of typhoons, rainfall, and storm surges.

Return Period (Years
) without considering the dependence between these flood drivers. Table 8 also shows a comparison of flood risk with and without considering the dependence between typhoons, rainfall, and storm surges. From Table 8, we can conclude that P ∩ (t, h, z) is much higher than P * ∩ (t, h, z) in different return periods. Here, P ∩ (t, h, z) denotes the co-occurrence probability that typhoon T, rainfall H, and storm surge Z all exceed a certain magnitude. P * ∩ (t, h, z) denotes P((T > t) ∩ (H > h) ∩ (Z > z)) without considering the dependence between these flood drivers. This indicates that ignoring the dependence between these flood drivers may inappropriately characterize the flood risk in coastal regions and can lead to underestimating it. These findings are in agreement with the results found by Salvadori et al. [53]. Figure 10 and Table 9 show the co-occurrence probability of rainfall and storm surges under different typhoon RP conditions. The co-occurrence probability of rainfall and storm surges increases by at least 200% when a typhoon occurs (see Table 9). Furthermore, with an increase in typhoon RP, the conditional probability of rainfall and storm surges increases rapidly. For example, the co-occurrence probability of a 50-year RP rainfall and 50-year RP storm surge is only 0.004 when there is no typhoon, and it increases to 0.036 under 5-year RP typhoon conditions. When a 50-year RP typhoon occurs, it increases to 0.225. Furthermore, the probabilities P(H > h) and P(Z > z) significantly increase when a typhoon has occurred (see Table 9). Therefore, the coastal urban flood management strategy should pay more attention to the joint impact of rainfall and storm surges on flood risk when a typhoon has occurred. Table 9. Probabilities P(H > h|T > t) , P(Z > z|T > t) , P(H > h, Z > z), and P(H > h, Z > z|T > t) . typhoon RP, the conditional probability of rainfall and storm surges increases rapidly. For example, the co-occurrence probability of a 50-year RP rainfall and 50-year RP storm surge is only 0.004 when there is no typhoon, and it increases to 0.036 under 5-year RP typhoon conditions. When a 50-year RP typhoon occurs, it increases to 0.225. Furthermore, the probabilities  ( ) P H h and  ( ) P Z z significantly increase when a typhoon has occurred (see Table 9). Therefore, the coastal urban flood management strategy should pay more attention to the joint impact of rainfall and storm surges on flood risk when a typhoon has occurred.

Impact of Typhoons on Flood Risk Caused by Rainfall and Storm Surges
The impact of typhoons on flood risk caused by rainfall and storm surges was estimated based on an urban flood inundation model and an expected annual damage model. In this study, expected annual damage (EAD) was then found to express the flood risk, which can be calculated by integrating the flood damage function with the probability function. Flood extent and depth are the most important indicators of flood damage, often denoted as depth-damage curves [54]. However, a credible regional depth-damage curve is difficult to obtain due to the complexity of urban contexts [55]. Therefore, we used unit flood damage and inundation depth for cost estimation [32]. The unit flood damage on Haidian Island is from Lian et al. [32], which was selected as the average economic loss in unit inundation depth of the four rainfall events, RE-1 (12 October 2008), RE-2 (5 August 2009), RE-3 (28 September 2009), and RE-4 (12 October 2009). The inundation depth and economic loss in the four rainfall events are shown in Table 10. Table 10. Inundation depth and economic loss for different rainfall events [32]. Inundation depth was calculated by PCSWMM. The calculation steps of inundation depth are introduced in detail in Appendix B. The maximum inundation depths in different RPs of rainfall and storm surges are presented in Figure 11. Flood damage was calculated by maximum inundation depths and unit flood damage. Table 11 shows that flood damage increases quickly with the increase of return period, and P(H,Z|T) is much more than P(H,Z). After calculation, EAD is 0.82 million dollars when there is no typhoon, and it rises to 3.27 million dollars when a typhoon has occurred. This indicates that typhoons greatly increase the flood risk in coastal zones. The main reason is the increase of the occurrence probability of rainfall and storm surges in typhoon conditions.

Conclusions
In this study, the joint risk of rainfall and storm surges during typhoons was investigated by integrating principal component analysis, copula-based probability analysis, urban flood inundation model, and flood risk model methods. First, principal component analysis was used to quantify typhoons by integrating MSW, center pressure, and distance between the typhoon center and the study area. Following this, the Gumbel copula was found to be a robust and proper function for the joint probability of typhoons, rainfall, and storm surges. Then, the joint probability, co-occurrence probability, and conditional probability were indicated by the Gumbel copula. The results of joint probability indicate that ignoring the dependence between these flood drivers may inappropriately characterize the flood risk in coastal regions and can lead to underestimating it. For conditional probability, the co-occurrence probability of rainfall and storm surges increases by at least 200% during typhoons. Therefore, coastal urban flood management should pay more attention to the joint impact of rainfall and storm surges on flood risk when a typhoon has occurred. Furthermore, when a typhoon has occurred, The expected annual damage (EAD) increases from

Conclusions
In this study, the joint risk of rainfall and storm surges during typhoons was investigated by integrating principal component analysis, copula-based probability analysis, urban flood inundation model, and flood risk model methods. First, principal component analysis was used to quantify typhoons by integrating MSW, center pressure, and distance between the typhoon center and the study area. Following this, the Gumbel copula was found to be a robust and proper function for the joint probability of typhoons, rainfall, and storm surges. Then, the joint probability, co-occurrence probability, and conditional probability were indicated by the Gumbel copula. The results of joint probability indicate that ignoring the dependence between these flood drivers may inappropriately characterize the flood risk in coastal regions and can lead to underestimating it. For conditional probability, the co-occurrence probability of rainfall and storm surges increases by at least 200% during typhoons. Therefore, coastal urban flood management should pay more attention to the joint impact of rainfall and storm surges on flood risk when a typhoon has occurred. Furthermore, when a typhoon has occurred, The expected annual damage (EAD) increases from 0.82 million dollars to 3.27 million dollars. This indicates that typhoons greatly increase the flood risk in coastal cities.
These messages are useful for practical design and planning, since all flood hazards, such as typhoons, rainfall, and storm surges, are considered. However, the study also has certain limitations. First, since the main objective of this study is to evaluate the impact of typhoons on the joint risk of rainfall and storm surges, the impact of typhoon size and typhoon conditions on rainfall distribution was not taken into consideration and needs to be explored in our future work. Second, the limited socioeconomic data in the study area restricted monetizing flood damage. Referring to Lian et al. [32], flood damage was defined as a function of inundation depth and flood unit cost in this study. Future research work could focus on improving the accurate quantification of flood damage. Furthermore, the impact of the uncertainty of copulas on flood risk estimation is also an important research focus in our future work. Acknowledgments: Our cordial gratitude should be extended to the editor and three anonymous reviewers for their professional and pertinent comments and suggestions, which were greatly helpful for further quality improvement of this manuscript.

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

Appendix A. Introduction of Copula Functions
In this study, two elliptical copulas (Gaussian copula and Student's t copula) and three Archimedean copulas (Gumbel copula, Clayton copula, and Frank copula) were employed to build joint distributions of typhoon-rainfall, typhoon-storm surge, and rainfall-storm surge. Common bivariate copula functions are presented in Table A1. Table A1. Common bivariate copula functions.

Copulas C(u,v)
Gumbel The trivariate Gaussian copula, Student's t copula, and Gumbel copula were used to establish the joint distribution of typhoons, rainfall, and storm surges. Trivariate copula functions are as follows.
2. Trivariate Student's t-copula: where T v −1 (·) is the quintile function of Student's t distribution function with v degrees of freedom, m is the variable dimensions, and w = [w 1 , w 2 , · · · , w m ] T is the integrand variable matrix.
3. Trivariate Gumbel copula: where θ is a parameter of the copula function. The parameters of the above distributions were estimated by the maximum likelihood method. The Kolmogorov-Smirnov (K-S) test was used to assess goodness of fit. The best-fit copulas were selected by the Akaike information criterion (AIC).
The Kolmogorov-Smirnov statistic D is defined as follows: where C k is the theoretical distribution of the measured samples, i is the serial number of the measured samples in ascending order, and n is the number of samples. When the statistic D is less than the critical value D α , the test is accepted. The formula of AIC is defined as follows: where N is the number of samples, RMSE is the root-mean-square error of samples, and k is the number of variables. A smaller AIC value indicates a better fit.

Appendix B. Urban Flood Inundation Model on Haidian Island
PCSWMM, developed by CHI, Canada, is widely applied in the simulation of surface hydrological processes and drainage network flows [33][34][35][36]. Unlike SWMM, which can only simulate 1D conduit flow and river flow, PCSWMM can accurately simulate unsteady 2D surface flow above ground through a 2D floodplain model. In this study, the urban flood inundation model on Haidian Island was established by PCSWMM. The calculation steps of inundation depth are as follows: (1) Prepare data. The basic data of the urban flood inundation model include digital elevation model (DEM) data, river data (e.g., river cross-section and distance between cross-sections), drainage data (e.g., junction depth, conduit size, etc.), and obstruction data (i.e., building data). DEM data ( Figure 1) were accessed from the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences (http://www.resdc.cn/Default.aspx), and the DEM resolution is 25 m. River data and drainage data were provided by the Haikou Municipal Water Authority. The obstruction data were extracted from satellite remote sensing images through Environment for Visualizing Images (ENVI) software. (2) Establish 1D model. The dynamic wave method was used to calculate hydraulic simulations of drainage conduits. Based on the specific circumstances of the study area, for example, the study area is a small catchment and there is not enough soil data; the Horton model was used to calculate infiltration, since it is suitable for small catchments and areas where soil data are lacking. The parameters of the PCSWMM model were calibrated by trial and error and the recommended ranges from relevant references [37,56]. (3) Establish 2D model. The size of the 2D grid is 25 m × 25 m, and the study area was divided into 23,578 grids. The 1D conduit model and 2D floodplain model were integrated by the orifice connection method [57]. (4) Calibrate urban flood inundation model. Calibration data for storm events was based on actual flood inundation data during the Rammasun typhoon in July 2014, which was obtained through field investigation. Figure A1 shows observations and calculation inundation depths in different points of the study area. The observation points are indicated in Figure 1. The error values between observations and calculations are smaller than 0.1 m, which indicates that the model is reliable.
follows: (1) Prepare data. The basic data of the urban flood inundation model include digital elevation model (DEM) data, river data (e.g., river cross-section and distance between cross-sections), drainage data (e.g., junction depth, conduit size, etc.), and obstruction data (i.e., building data). DEM data ( Figure 1) were accessed from the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences (http://www.resdc.cn/Default.aspx), and the DEM resolution is 25 m. River data and drainage data were provided by the Haikou Municipal Water Authority. The obstruction data were extracted from satellite remote sensing images through Environment for Visualizing Images (ENVI) software. (2) Establish 1D model. The dynamic wave method was used to calculate hydraulic simulations of drainage conduits. Based on the specific circumstances of the study area, for example, the study area is a small catchment and there is not enough soil data; the Horton model was used to calculate infiltration, since it is suitable for small catchments and areas where soil data are lacking. The parameters of the PCSWMM model were calibrated by trial and error and the recommended ranges from relevant references [37,56]. (3) Establish 2D model. The size of the 2D grid is 25 m × 25 m, and the study area was divided into 23,578 grids. The 1D conduit model and 2D floodplain model were integrated by the orifice connection method [57]. (4) Calibrate urban flood inundation model. Calibration data for storm events was based on actual flood inundation data during the Rammasun typhoon in July 2014, which was obtained through field investigation. Figure A1 shows observations and calculation inundation depths in different points of the study area. The observation points are indicated in Figure 1. The error values between observations and calculations are smaller than 0.1 m, which indicates that the model is reliable. Figure A1. Calibration of the PCSWMM model.