Probable maximum tropical cyclone parameters for east and west coast of India

Storm surge simulation models require certain parameters for evaluating the worst possible event that could occur at a site with respect to a certain return period. The most significant probable maximum tropical cyclone parameters are pressure deficiency at the centre (ΔP) and maximum wind speed (Vmax) during the cyclone. In this study, the probable maximum tropical cyclone parameters that would yield the maximum probable storm surge along the Indian coasts of Arabian Sea and Bay of Bengal are estimated. Datasets are created based on various existing data archives for each basin. The datasets are subjected to extreme value analysis for determining the ΔP and Vmax parameters. The data are fitted to various probability distributions (Gumbel, Fréchet, Weibull and Log-normal) whose parameters (scale, shape, and location parameters) are estimated using graphical (least square fit) and numerical (order statistics approach) methods. A mean recurrence interval of 1000 and 10,000 years is considered for strategic structures. The best fit distribution and its parameters are obtained based on goodness of fit criteria. The resulting ΔP and Vmax are compared with theoretical maximum cyclone parameter values of each basin and revised till an optimal set of values are reached. The analysis shows that ΔP and Vmax values for Arabian Sea and Bay of Bengal are best represented by Weibull distribution. The estimated parameters are useful input to a storm surge model to determine the design basis flood level for the strategic coastal sites.


Introduction
Strategic and industrial installations located along the coast are required to be designed to withstand the adverse flooding phenomena caused by storm surges, precipitation, and other natural disasters. Hence, it is necessary to determine a design basis event for each flood inducing phenomena and provide appropriate protection to ensure the safety of a site as specified by Atomic Energy Regulatory Board (AERB 2014). The standard procedure is to calculate the probable maximum water level from a design basis flood (DBF) event. The design basis flood level (DBFL) is arrived based on the probable maximum water level. In the case of coastal locations, the maximum water level is a combination of highest astronomical tide, probable maximum storm surge and wave run-up (AERB 2002). The probable maximum storm surge is derived from past events with respect to a return period. This is not feasible in all areas because accurate storm surge data (the difference between the tide level and the total water level) is not available in many gauge locations of the site (IAEA 2011). Hence, the probable maximum storm surge is obtained from simulations of a probable maximum tropical cyclone (PMTC) (AERB 2002;USNRC 2011).
PMTC is defined as the hypothetical storm event that generates the maximum possible sustained wind speed in the considered study area. It is defined by a set of parameters out of which the most significant are (AERB 2002;IAEA 2011): 1. Pressure difference between centre of cyclone and the periphery (ΔP in millibar or hPa) and 2. Maximum sustained wind speed during the cyclone (V max in kmph).
The PMTC parameters are estimated using deterministic or probabilistic method. The deterministic method employs a mathematical approach to evaluate ΔP by balancing pressure gradient force and force of gravity under hydrostatic approximations. This method requires site-specific information like height and temperature of the tropopause, sea surface temperature and moisture distribution that are not readily available in all areas. As a result, the probabilistic approach is recommended (AERB 2008). Probabilistic method or extreme value analysis (EVA) is necessary to determine the most likely value of any variable, which can then be used in design of various infrastructures and facilities (Gupta et al. 2018). It is performed on datasets using statistical techniques for different recurrence intervals or return periods based on the sensitivity of the site. In the case of strategic installations such as nuclear power plants and nuclear research facilities situated near the coast a return period of 10,000 years and 1000 years is suggested, respectively (AERB 2014). The probability of each cyclone event in the dataset is determined using a plotting position formula. The dataset is then fitted to a set of prescribed extreme value probability distributions (AERB 2008;Goda 1990;Kudale 2010). The fitting is done by estimating the parameters of the distribution using various methods like graphical (least square fit), numerical (order statistics), method of moments and method of maximum likelihood (Lieblein 1974;Palutikof et al. 1999). The evaluated distribution parameters are applied in the linearized empirical equation of each distribution and extrapolated for higher return periods or mean recurrence interval (MRI, also known as average recurrence interval (ARI)). The distribution parameters generated by the preceding approaches are tested to reflect the real situation by comparing the linear plot to actual data. This is achieved using goodness of fit (GOF) tests with a best fit criterion to identify the best fitting probability distribution parameters and its estimation technique.

3
PMTC parameters have been estimated using statistical analysis (probabilistic method) for strategic sites (nuclear power plants or associated research labs) around the world (Tai et al. 2014;Wang 1982;Yin et al. 1995;Yin and Wang 1991;Zhao 2009). Statistical analysis of India Meteorological Department (IMD) cyclone data (from 1971 to 2007, 36 years) collected along the Tamilnadu coast, was carried out by employing least square fit technique on Gumbel distribution to identify 50-year return period ΔP as 66 hPa (Rao et al. 2010). The analysis does not provide any distribution parameters to be extrapolated for 1000 and 10,000 years. The data used here is also insufficient to make any reasonable estimate for such high return periods or low frequency events. Hence, a separate analysis is required to define the best fit distribution and its associated parameters so that they could be used in determining PMTC parameters for any return period. Atomic Energy Regulatory Board (AERB) suggests the usage of least square fit and order statistics approach to determine best fit distribution and PMTC parameters. The most adverse 100-year return period cyclone event for Bay of Bengal (based on 127 year observed cyclone data) is ascertained to have a ΔP and V max of 98 hPa and 265 kmph (Pandey 2020). Similarly, the observed maximum ΔP for Arabian Sea is observed to be 66 hPa (1975 Porbandar cyclone) based on 100 years data from 1917 to 2016 (BMPTC 2018). The maximum observed 100-year return period sustained wind speed for Arabian Sea is determined as 216 kmph (GSDMA 2014;Poulose et al. 2020). The analysis relies on observed data which is not sufficient for PMTC parameter estimation of 1000 and 10,000 years return period. This makes EVA based estimation of PMTC parameters the most viable method as they do not require extensive data records.
In this study, EVA is employed to determine the PMTC parameters (ΔP and V max ) for Arabian Sea and Bay of Bengal by fitting the peak event values to extreme value distributions (AERB 2008;Goda 1990;Gumbel 1958). The Best track and cyclone e-Atlas data from IMD are compiled into a database to derive the probability distribution parameters using various estimation techniques. The estimated PMTC parameters are compared with the theoretical maximum cyclone parameter of each basin. The resultant values could be used as input in storm surge models to yield the DBFL of a site located along the east and west coast of India.

Methodology for probabilistic method
The general procedure for estimation of extreme value for meteorological parameters comprise of the following steps: • The dataset is analysed in site representativeness, data reliability and completeness. • There are two approaches in creating the dataset namely, Block maxima (BM) and Peak-over-threshold (POT). In BM approach, the peak value of each year for a prescribed record period is subjected to EVA. In POT approach, a threshold value is selected and any event of magnitude higher than the threshold is considered for EVA (Coles et al. 2001;Subramanya 2008). • The most appropriate or best fitting statistical distribution and its associated distribution parameters is determined for each PMTC parameter based on GOF tests.
• The estimated PMTC parameter (from best fit distribution) for suggested return period is compared with the theoretical maximum cyclone parameters (Sect. 2.1.6). The distribution parameters are revised to determine the optimal PMTC parameters.
In order to perform the above tasks, a database with sufficient data is required. The flowchart in Fig. 1 details the step-by-step procedure to determine the PMTC parameters by probabilistic method.
(1) where x is the considered PMTC parameter value; a is location parameter; g is shape parameter and b is scale parameter. The value of g for Weibull distribution could be 0.75, 1.0, 1.4 or 2.0 (Goda 1990). The non-exceedance probability is linearized to yield the reduced variate ( Y ) of the distribution. In the case of Gumbel, Fréchet and Weibull distributions, this is achieved by taking negative natural logarithm two times on both sides of the empirical distribution equation. The equations could be called as linearized empirical distribution equation and are given in Eqs. (5), (6) and (7) for Gumbel, Fréchet and Weibull distributions, respectively.
The Gumbel and Weibull distribution parameters are calculated directly from the observed values, but the Fréchet distribution parameters are evaluated using logarithmic transformation on the variables. The Log-normal distribution parameters are evaluated by simultaneously solving below Eqs. (8), (9), (10) and (11) (Abramowitz and Irene 1964).
The probability of each cyclone event in the database is estimated using the plotting position formula and is applied in the above linearized empirical distribution equations.

Plotting position formulation
Weibull formula is the most common and best suited plotting position formula for EVA (AERB 2008;Benjamin and Cornell 2014;Gupta et al. 2018;Makkonen 2006;Ologhadien 2021;Weibull 1939Weibull , 1951. The data x i are ranked in ascending order and the non-exceedance probability P x i for each event is evaluated using Eq. (12): where i is the rank of the data and N is the total number of data points. This formula is suggested by AERB for both Gumbel and Fréchet distribution. While the Gringorten plotting position formula as given in Eq. (13) is also suggested for Gumbel distribution (Goda 1990;Gringorten 1963;Kudale 2010) The different plotting position formulations suggested by above literatures could affect the resultant distribution parameters. Hence, both Weibull and Gringorten plotting position formula were applied in the present study to assess the most suitable one. The Gumbel distribution using AERB prescribed Weibull plotting position formula is referenced as Gumbel-AERB. While the Gumbel distribution using Goda (1990) prescribed Gringorten plotting position formula is referred as Gumbel-Goda. The plotting position formulas for Weibull and Log-normal distribution are given by Petrauskas and Aagaard (Eq. 14) and Blom (Eq. 15) formulations, respectively (Goda 1990;Kudale 2010;Petruaskas and Aagaard 1970).

Parameter estimation techniques
Parameter estimation is required to examine the fitness of a distribution to a parameter. There are multiple techniques for parameter estimation, namely method of moments, maximum likelihood, least squares, probability plot correlation coefficient (PPCC) plot and probability plot (Bury 1999;Filliben 1975). AERB recommends two approaches namely graphical (or least square fit) and numerical (or Lieblein) technique for determining the parameters of a distribution (AERB 2008).

Least square fit method
The least square fit is a technique used to determine the parameters of a distribution by minimizing the sum of the squares of the residuals between the observed and expected values. The probability of each event is estimated using plotting position formula of each distribution as given in Sect. 2.1.2. The estimated probability is then applied in the linearized empirical distribution equation. The reduced variate is evaluated and plotted against the considered variable. A trendline is fitted to the data and a linear equation (with slope m and intercept c) is obtained. The slope and intercept values of the equation refers to the scale and location parameters of the Gumbel, Weibull and Log-normal distributions. In Fréchet distribution, the shape and scale parameters are determined using Eqs. (16) and (17) as shown below: The distribution parameters are also estimated using Lieblein technique and compared with least square fit results.

Lieblein technique
It is also known as the Order Statistics Approach and is suggested especially for EVA (AERB 2008;Coles et al. 2001;Gupta et al. 2018;Lieblein 1974). The data are categorised into a main subgroup (set of 6) and the remaining data is brought under remaining subgroup. The main subgroup is treated like a matrix where the row elements are arranged in an ascending order and summed up along each column. The summed values are normalised using appropriate weights to determine the component of distribution parameters for the main sub-group. Similarly, the component of distribution parameters for the remaining sub-group is obtained. Then both components are multiplied with proportionality factors and summed to yield the distribution parameters. The detailed procedure, formulae and weights required for this method are provided in the AERB guidelines (AERB 2008).

Goodness of fit tests
The trendline generated values of least square fit and Lieblein technique methods are compared with actual data to determine the best fitting parameters using GOF tests. The GOF test describes how well a pre-defined empirical distribution fits to a set of observations. It typically measures the discrepancy between observed values and the values expected from the distribution. In this analysis, sum of squared error (SSE) (Hosmer et al. 1997) and correlation coefficient (R) (Kessler and Neas 1994;Vivekanandan 2017;Vivekanandan and Ramesh 2017) are used as a measure of goodness of fit. The correlation coefficient is used to determine the best fit probability distribution. Each distribution has two sets of parameters yielded by least square fit and Lieblein techniques. The best fitting parameters of the best fit distribution are selected based on the SSE value. The distribution parameter that corresponds to the lowest value of SSE is adopted. The selected distribution parameters are applied in the linearized empirical distribution equation (given by Eq. (5), (6), (7), (8)) while the reduced variate is computed with respect to an MRI.

Mean recurrence interval
The non-exceedance probability P(x) in the case of limited data availability for a design basis event having T years MRI is given by Eq. (18): where is the average number of storms per year; M is the number of storms in the time period and K is the time period for which data is collected (AERB 1998(AERB , 2014Kudale 2010).

Theoretical maximum cyclone parameters
The MRI-based ΔP and V max values have to be validated with a peak value (based on the theoretical understanding of a cyclone) called theoretical maximum cyclone parameters or maximum potential intensity (MPI). The theoretical maximum cyclone parameters were determined using the sea surface temperature (SST), relative humidity of boundary layer and atmospheric thermodynamic profiles (Emanuel 1987(Emanuel , 1988(Emanuel , 1995Miller 1958).
Empirical relationships have been established between V max of MPI and SST for North Atlantic Ocean (Demaria and Kaplan 1994a), Eastern North Pacific Ocean (Whitney and Hobgood 1997) and Western North Pacific Ocean (Zeng et al. 2007). Each method has been site specifically developed based on its SST. V max of MPI also depends on duration of storms (Demaria and Kaplan 1994a) and vertical wind shear (Demaria and Kaplan 1994b). Cyclones attain maximum strength when the origin is on warmer SSTs and have longer duration over the sea (Kotal et al. 2009). The values are not with respect to any return period but provides the maximum potential cyclone intensity. These theoretically derived values aid in providing a threshold to prevent over-estimation of PMTC parameters. The ΔP and V max obtained for 10,000-year return period are compared with the maximum cyclone parameters. If the estimated value of a PMTC parameter is greater than the theoretical maximum value, the next best distribution parameters are used to re-evaluate the 10,000-year return period PMTC parameter. It is again compared with the theoretical maximum value based on the aforementioned condition. This process is carried out till the estimated 10,000-year return period PMTC parameter is within the possible maximum value. Once the 10,000-year return period PMTC parameter is determined, the corresponding distribution parameters are used to evaluate the PMTC parameter for 1000-year return period.

Study area
Majority of nuclear power plants and associated research facilities in India are situated along the coasts of the Arabian Sea and Bay of Bengal. Arabian Sea and Bay of Bengal are part of the Northern Indian Ocean that is bordered on three sides: north, west and east, with the Indian Ocean on the south. The Arabian Sea has a surface size of 3,862,000 sq km and a maximum depth of 4652 m. Similarly, the Bay of Bengal has a surface size of 2,600,000 sq km and a maximum depth of 4694 m. According to IMD cyclone data from 1891 to 2018, the number of atmospheric anomalies along the Arabian Sea and Bay of Bengal were 222 and 1211, respectively (RSMC 2018). This corresponds to approximately 2 cyclones per year in the Arabian Sea and 10 cyclones per year along the Bay of Bengal. This proves that Arabian Sea and Bay of Bengal differ in terms of meteorological conditions; and hence, PMTC parameters would also vary for each coast requiring independent analysis of each basin.

Data
The data required for this analysis are the ΔP in millibar and V max in kmph for each cyclone event retrieved from an archive. The India Meteorological Department (IMD) is responsible for monitoring and tracking atmospheric anomalies (like cyclones) in the Bay of Bengal and Arabian Sea. IMD classifies and provides the data associated to anomalies (cyclones) based on their sustained wind speed (IMD 2003). Best track, cyclone e-Atlas data and synoptic charts provided by IMD are the data sources considered in this investigation.

Best track data
The Regional Specialised Meteorological Centre (RSMC), New Delhi of IMD publishes six hourly best track data comprising of various cyclonic disturbances over the considered study area since 1982 (RSMC 2018). The individual PMTC parameter values for cyclonic episodes from 1982 to 2018 were considered in this analysis. The total number of cyclonic disturbances in the Arabian Sea and Bay of Bengal is 75 and 242, respectively. The data presented in the best track archive of IMD are available only for a small period which is not sufficient to determine the PMTC parameters for 1000 or 10,000 years return period with reasonable accuracy. Hence, cyclone e-Atlas data are also considered in addition to the best track data.

Cyclone e-Atlas
Cyclone e-Atlas contains cyclone and depression tracks over the Northern Indian Ocean since 1891(RSMC 2011a). Figure 2 shows the archived cyclone tracks for Arabian Sea. Similar archive is also available for Bay of Bengal. Cyclone e-Atlas tracks are classified based on Beaufort wind scale as depression (D), cyclonic storm (CS) and severe cyclonic storm (SCS) whose wind speed varies up to 61 kmph, 62-88 kmph and > 88 kmph, respectively (RSMC 2011b). This system of classification is not in line with the current IMD cyclone classification. Hence, the cyclone e-Atlas data are reclassified, in this study, based on best track sustained wind speed to determine a generalised maximum ΔP for each cyclone (obtained from e-Atlas) category. The V max parameter for each cyclone event in cyclone e-Atlas data is obtained using the correlation between ΔP and V max in the current IMD classification. The ratio of V max to ΔP is evaluated for best track data under each category and normalised to match cyclone e-Atlas categories. The values of ΔP, ratio of V max to ΔP and computed V max are provided in Table 1 for each category in cyclone e-Atlas.

Datasets
The data obtained from IMD's cyclone e-Atlas and best track are compiled to create five datasets. The first dataset comprises of the best track data values for ΔP and V max from 1982 to 2018 in the Arabian Sea. The second dataset is created by combining cyclone e-Atlas data (generic value of ΔP and V max is applied to each category) from 1891 to 1981 with best track data for Arabian Sea. This dataset has 46 severe cyclonic events with very high wind speeds that might lead to over-estimation. As per AERB recommendations, synoptic charts can be used for storm data hindcasting (AERB 2002). The Indian daily weather report of IMD has included synoptic charts depicting cyclone pressure distribution since 1891. Based on these synoptic charts, ΔP and V max are derived for the Arabian Sea (CWPRS 2007). The ΔP and V max of SCS events obtained from synoptic charts were compared with the 46 SCS events of Dataset 2. The third dataset was generated by replacing the 46 SCS events with values based on synoptic charts for Arabian Sea. The fourth dataset is made of best track data for Bay of Bengal. The fifth dataset is developed combining cyclone e-Atlas data from 1891 to 1981 with best track data for Bay of Bengal. The study is carried out using five datasets, which are referred to as dataset 1-5 as shown in Table 2.

Reduced variate
The reduced variate for different extreme value distributions and datasets are calculated for 1000 and 10,000-year return period and tabulated in Tables 3 and 4, respectively. The reduced variates of Gumbel-AERB; Gumbel-Goda; Fréchet and Weibull (g = 1.0)   distributions are same but the distribution parameters are different. Hence, the PMTC parameters obtained from the analysis will vary.

PMTC parameters based on GOF
The datasets are subjected to probabilistic analysis as per the methodology described in Sect. 2.1. The distribution parameters were obtained for both PMTC parameters using least square fit and Lieblein technique in Bay of Bengal and Arabian Sea. It was observed that dataset 3 provided better correlation coefficient among all datasets considered for Arabian Sea. Similarly, dataset 4 performed best for Bay of Bengal. Tables 5 and 6 show the estimated 10,000-year return period values of ΔP and V max using least square fit and Lieblein technique on dataset 3. Tables 7 and 8 shows the 10,000-year return period values of ΔP and V max using dataset 4. The correlation coefficients are found to be same for a probability distribution with parameters estimated by both least square fit and Lieblein technique. This is because correlation coefficient is not expected to vary due to any linear transformation on the data. However, it differs for different probability distributions. The best fit distribution has two sets of distribution parameters from which the best one is determined using the lowest SSE value. When comparing the distribution parameters produced using least square fit method with that of Lieblein technique, the former consistently has a lower SSE value. The best fitting distribution in both basins based on GOF tests is found to be Weibull distribution with shape parameter 0.75 and 1.0 for ΔP and V max , respectively. The ΔP and V max obtained for

Theoretical maximum cyclone parameters for Arabian Sea and Bay of Bengal
The cyclone potential of each basin was determined to define the theoretical maximum cyclone parameters. They were established for Bay of Bengal based on existing linear empirical relationships between mean SST in ºC and V max in knots. The impact of vertical wind shear was found to be detrimental while duration of storm increased the theoretical maximum cyclone parameters in the Bay of Bengal (Kotal et al. 2009). The maximum mean recorded SST of 29.19ºC occurred during 2003-2005(Khokiattiwong and Yu 2012. The SST over Bay of Bengal was projected to increase by 2.0 °C to 3.5 °C by the end of this century (Rajalakshmi and Achyuthan 2021;Vivekanandan et al. 2016). Applying the mean maximum projected SST (of 31.94 °C) in the empirical relationship developed by Kotal et al. (2009), the V max of MPI is evaluated as 400kmph. Similarly, the ΔP is estimated as 102 hPa. As cyclone ferocities are identical in both basins, the same theoretical maximum cyclone parameters are applied to Arabian Sea.

Revised PMTC parameter based on theoretical maximum cyclone parameter
The theoretical maximum cyclone parameters (Sect. 3.3) were compared with the 10,000year return period PMTC parameters (Sect. 3.2). The ΔP and V max of Fréchet distribution (Tables 5, 6, 7, 8) are very high compared to the theoretical maximum cyclone parameter and are physically impossible to occur in the considered study area. Hence, values obtained from Fréchet distribution are neglected. The next best fit distribution parameters for Arabian Sea corresponds to Weibull distribution with shape parameter 1.0 and 1.4 for ΔP and V max , respectively. The values are found to be 132 millibar and 317 kmph ( Table 5). The ΔP value was observed to be still higher than the theoretical maximum cyclone parameter value. Hence, the next best distribution is opted for ΔP. The process is repeated until the 10,000-year return period values becomes comparable to the theoretical maximum cyclone parameter values, while adhering to GOF criteria.

Recommendation of PMTC parameters
The final distribution for ΔP and V max in the Arabian Sea is given by Weibull distribution having shape parameter 1.4. Similarly, the best fit distribution for ΔP and V max in the Bay of Bengal is given by Weibull distribution having shape parameter 2.0. Figure 3a shows the relationship between ΔP and its reduced variate for best fit distribution of Arabian Sea with parameters estimated using graphical (GM) and numerical (NM) methods. Similarly, Fig. 3b provides the relationship between V max and its reduced variate for best fit distribution of Arabian Sea. Figure 4a describes the relationship between ΔP and its reduced variate for best fit distribution of Bay of Bengal. Figure 4b indicates the relationship between V max and its reduced variate for best fit distribution of Bay of Bengal. The linear equation of graphical method provided in Figs. 3 and 4 is used to determine the finalised PMTC parameter for 1000-and 10,000-year return periods. The finalised PMTC parameters obtained for Arabian Sea and Bay of Bengal are tabulated in Table 9. The results obtained by this method is expected to improve, based on addition of data in the future. Hence, it is necessary that PMTC parameters must be evaluated for every strategic site using latest data available.

Conclusion
The PMTC parameters (pressure difference between the cyclone centre and its periphery (ΔP) and the maximum sustained cyclonic wind velocity (V max )) for the Arabian Sea and Bay of Bengal were evaluated in this study. Though previous estimations exist, a new analysis was required to define the most possible PMTC parameters for MRI of 1000 and 10,000 years. The study compared PMTC parameters from datasets with various record durations, probability distributions, plotting position formulations and parameter estimation technique with the MPI of each basin. Datasets were compiled from IMD sources (cyclone e-atlas and best track data) and subject to PMTC estimation by probabilistic method. It is evident from this analysis that a larger data record would undoubtedly aid in finetuning the final values for such higher MRI. Gumbel distribution is generally best suited for this type of extreme value analysis. The consideration of other distributions like Weibull and Lognormal has assisted in selecting the best fit distribution. This investigation proves that PMTC parameters for Arabian Sea and Bay of Bengal were best represented by Weibull distribution which contradicts previous studies. The significance of employing multiple plotting position formulations for Gumbel distribution was also examined. The Gringorten formulation better represented this data than the AERB suggested Weibull formulation. Conventionally, the probability distribution parameters estimated using Lieblein technique have better fitted to random natural events than least square fit estimates. The PMTC parameters for Arabian Sea and Bay of Bengal obtained using Lieblein technique does not follow this trend. The Lieblein technique derives the probability distribution parameters from the data record yielding only one set of parameters for all distributions. The application of a single set of scale and location parameters for all distributions reduces the goodness of fit. Hence, these parameters have been consistently underperforming when compared with their least square fit counterparts. The best-fit probability distribution for PMTC parameters (ΔP and V max ) in the Arabian Sea and Bay of Bengal is determined to be the Weibull distribution with a shape parameter of 1.4 and 2.0, respectively. Similar methodology could be applied to determine PMTC parameters elsewhere. However, the best fit distribution and its associated parameters could vary based on the data considered. This study presents a basic methodology to determine the most possible PMTC parameters. Future research could be focussed on defining a methodology to determine the probable