Extreme wind return periods from tropical cyclones in Bangladesh: insights from a high-resolution convection-permitting numerical model

We use high resolution (4.4 km) numerical simulations of tropical cyclones to produce exceedance probability estimates for extreme wind (gust) speeds over Bangladesh. For the first time, we estimate equivalent return periods up to and including a 1-in-200 year event, in a spatially coherent manner over all of Bangladesh, by using generalised additive models. 10 We show that some northern provinces, up to 200 km inland, may experience conditions equal to or exceeding a very severe cyclonic storm event (maximum wind speeds in ≥ 64 knots) with a likelihood equal to coastal regions less than 50 km inland. For the most severe super cyclonic storm events (≥ 120 knots), event exceedance probabilities of 1-in-100 to 1-in-200 events remain limited to the coastlines of southern provinces only. We demonstrate how the Bayesian interpretation of the generalised additive model can facilitate a transparent decision-making framework for tropical cyclone warnings. 15

In this study we seek to improve our understanding of the historical extreme gust speed hazard associated with recent TCs. To address the lack of observation data in this region, we use the latest generation Met Office regional model over the BoB, to 65 simulate 9 versions of 12 historical tropical cyclone cases representing 1979-2019. This generates spatially and temporally consistent, counterfactual simulations (relative to observed TC cases), albeit limited by the constraints of the model configuration and computational resources. This ensemble configuration enhances our understanding of how each cyclone may evolve if a similar event were to happen again. We combine the ensemble information in a spatially coherent manner to produce hazard maps at 4.4km resolution over Bangladesh for extreme wind (gust) hazards. Using Bayesian inference, we 70 estimate gust speed exceedance intervals (return periods) across all of Bangladesh, and demonstrate how this information can be directly integrated into a decision making framework.

Geospatial processing
Tropical cyclone simulations are derived from a 9-member ensemble for 12 historical events, using the latest generation Met Office Unified Model (Brown et al., 2012) convection-permitting regional atmosphere configuration, based on Bush et al. 75 (2020)hereafter referred to as RA2. The RAL2 4.4km domain avoids placing model boundaries over the Himalayas and covers Nepal, Bhutan, Myanmar, most of India, and parts of the Tibetan plateau. To ensure model stability over this mountainous terrain, the RAL2 model was run with a 30 second time-step. Further details of the regional modelling process and validation against the observations-based International Best Track Archive for Climate Stewardship (IBTrACS, Knapp et al., 2010Knapp et al., , 2018 and ERA5 reanalysis (Hersbach et al., 2020) can be found in Steptoe et al. (n.d.). In general, median peak 80 gust speeds from the RAL2 model ensemble are found to be 22 to 43 m s -1 faster compared to ERA5, but it is known that extreme gusts associated with vigorous convection in ERA5 are generally under-estimated, sometimes by a factor of two (Owens and Hewson, 2018). We use the ensemble output to first derive event 'footprints'a common method within the catastrophe modelling community to define peak hazard relating to a given event. In this case, footprints are based on the maximum wind gust speed achieved within each model run of 48 hours, that implicitly collapses the time dimension to leave 85 a 2D gust field in a longitude-latitude frame of reference. Although the original regional model data covers a significant portion of the BoB, we crop the data to approximately 87.5°E to 93.0°E and 20.5°N to 27.5°N.

Generalised Additive Modelling (GAM)
To summarise information from all 9 regional model ensemble member footprints into a coherent spatial summary of the tropical cyclone hazard, we use a generalised additive model (GAM), after Hastie & Tibshirani (1986), based on the R package 90 mgcv of Wood (2017), as a flexible spatial regression framework. GAMs are an extension of generalised linear modelling that use smooth functions of covariates to build a linear predictor and have previously been applied in similar geospatial natural hazard assessments, such as storm count data over Europe (Youngman and Economou, 2017), spatial prediction of maximum wind speed over Switzerland (Etienne et al., 2010) and return level estimation for U.S. wind gusts (Youngman, 2019). In each interaction (autocorrelation) present in the source data, and use the spatial dependence as a source of information.
For our purposes, we use a Gaussian location-scale (GLS) model family (Wood et al., 2016) to describe the natural logarithm (log) of the gust speed, where both the mean and the log of the standard deviation are smooth functions of predictorsin this case, longitude and latitude. Although other model families were trialled (such as generalized extreme value and gamma 100 distributions) the GLS family was found to have the best trade-off between computational efficiency and model fit. The general form of our GAM is: where ( ) is the response variable, namely log gust speed for each ensemble member in each grid cell = 1, . . . , N, N = 207,081. (a function of the mean) and (a function of the variance) are each defined as thin-plate regression splines (Wood, 2003) isotropic smooth functions of covariates and (longitude and latitude respectively). Each smooth 110 function requires a user-defined maximum amount of desired flexibility (wiggliness), traditionally quantified by the number of "knots". This flexibility is objectively penalised within mgcv to avoid over-fitting, while optimally explaining the trends in the data (Wood, 2003). Trial and error shows that (600) knots are required to construct thin-plate spline basis functions that avoid over smoothing given the resolution of the regional model data. Under this model formulation, the mean ( ) can be interpreted as an aggregated prediction across the ensemble members. 115 The smooth model parameters are estimated using restricted maximum likelihood (REML). However, once the model is fitted, it can be shown that it has a Bayesian interpretation. In particular, the coefficients of the smooth functions are assumed to have a multivariate Normal prior distribution, whose covariance matrix determines the wiggliness penalisation (see Wood, 2017 for further details). A Gaussian approximation of the posterior distribution for the coefficients then provides a multivariate Normal 120 distribution as the posterior (Gelman et al., 2013). In practice, once a GAM model is fitted to each named storm, under the Bayesian interpretation, we obtain 1000 simulations from the posterior distribution of the smooth function coefficients via random draws from a multivariate normal distribution (MVN). The MVN mean vectors are the REML coefficient estimates, and the MVN covariance is derived as a function of the covariance matrix of the sampling distribution of the model coefficients.
In Bayesian inference, sampling from the posterior distribution implies we can then derive samples from the posterior 125 predictive distribution of gust speed for each grid cell, ( ). The predictive distribution, a unique feature of Bayesian inference, fully quantifies estimation uncertainty and variability in gust speed across ensemble members. We take 1000 samples from the posterior predictive distribution and construct prediction intervals based on the empirical quantiles of these samples. To aggregate gust information from all ensembles of all named storms, we pool the 1000 posterior predictive simulations from each event into a total of 12,000 samples from the predictive distribution of gust speed across all 12 events. 130 Figure 1 summarises the key parts of this process.
Assessing the GAM specification for ( ) with detrended quantile-quantile (worm) plots (based on the method of Augustin et al., 2012), Figure 2 shows that generally storms are well represented. For some storms (such as Aila, BOB01, BOB07, Bulbul, Rashmi & TC01B) there is a tendency for the GAM to overestimate the tails of the distribution (positive kurtosis) 135 relative to the 4.4km data, as indicated by quantile-quantile plot points falling below the zero residual line. In these cases, the GAM will over-estimate extremes. Akash is the only storm where maximum gust speeds are likely to be underestimated in the GAM relative to the 4.4km data, but only for extreme upper-tail gust speeds. Checking for the consistency of variance over the range of predictor values, shows that the distribution of the residuals is stationary for both longitude and latitude (not shown). 140

Tropical Cyclones in Bangladesh
Aggregating the 12 historical tropical cyclones ensembles, Figure 3 shows the 50 th , 95 th and 99 th percentiles of the posterior predictive maximum gust speed distribution across Bangladesh. Based on historical cases, the provinces of Chittagong, Barisal and Khulna are most exposed to high wind speed associated with tropical cyclone gusts, whilst Sylhet and Rajshahi are least In addition to specific thresholds, exceedance probability curves ( Figure 5) summarise information for gust speeds up to 80 m s -1 (155 kn) for 18 of the most populated towns and cities in Bangladesh (grey lines) with four key cities highlighted. Coastal cities of Cox's Bazar and Chittagong are unsurprisingly the population centres most exposed to high gust speeds. Chittagong and Cox's Bazar are roughly 2.5 and 4.8 times more likely to experience tropical cyclones exceeding 'Very Severe' cyclonic storm conditions than Dhaka, for a landfalling cyclone. 165

Decision-making under uncertainty
By defining a loss function, it is possible to exploit the information in the Bayesian posterior predictive distributions to create a warning model based on decision theory (Lindley, 1991). Following Economou et al. (2016), defining a loss function L (a,x) to quantify the consequences of the various actions a (e.g. issuing warnings) that could be taken in the event of a landfalling TC of varying intensities x (see Table 1 for an example of four discrete gust categories), provides a method of mapping 170 predictive information onto an action. The optimum action a*, given some predictive information y (i.e. predictions of gust speed ( ) from the GAM), is one that minimises the loss L(a, x) taking into account the uncertainty in the predictive information, expressed as the probability of TC intensity x given predictive information y, p(x | y): a * = arg min ∫ ( , ) ( | ) 175 In practice, ( | ) can be easily computed from the predictive samples from the GAM, while the loss function L(a, x) is defined subjectively. Defining L(a,x) is a non-trivial process, as it should encapsulate the relative cost of false-positive (i.e. where action against a TC was taken, but the TC did not occur) and false-negative (i.e. where no action was taken, but the TC did occur) events. For the purposes of demonstrating the principle of this approach, we define a dummy loss function in Table  180 1, based on the four TC warning levels used in Bangladesh (WMO, 2018). Here relative loss is defined on a 100-point scale, where 0 equates to no loss associated with a given landfalling event, and 100 equates to maximum loss. Evacuation typically takes places at the 'Great Danger' level. Figure 6 illustrates the optimal warning that should be issued based on Table 1 and the range of gust speed information 185 summarised by our GAM. This can be interpreted as the default optimal action to take given an impending land-falling TC, and in reality the optimal action (a*) would be updated once new forecast information of the TC becomes available.
Interpretation is strongly conditioned by the loss function and the accuracy of the gust speed information, but our aim here is demonstrate a proof-of-concept transparent workflow that clearly translates hazards into actions. In this case, the northern extent of TC risk, as highlighted in Figures 2 and 3, is again reflected in the warning level, but in practice separate loss functions 7 could be defined for each province, or for different economic sectors of society. By understanding the exposure, vulnerability and decision-making process of each user, bespoke warnings could be issued.

Limitations
Despite the ensemble simulation framework, our analysis is still restricted to only 12 historical cases, which represent the recent 40-year period. The number of events was determined by the availability of source data (ERA5) for driving the regional 195 model (RA2), for TC events that made landfall over Bangladeshin this case limited to the period of ERA5 data availability, which at the time of analysis extended back to 1979. Given the relatively low ERA5 resolution (31 km), we selected TCs defined as at least a Category 1 event in the IBTrACS database , to be sure they would be identifiable within the low-resolution ERA5 data and could be downscaled by the RA2 model.

200
The initial conditions posed in the regional model play a significant role in determining the outcome of each event. In forecasting situations this is desirable behaviour: well-chosen initial conditions ensure the model retains a realistic representation of reality. Even though the modelling domain that produced these 4.4km data had the freedom to deviate in a physically plausible way (see Steptoe et al., n.d.), it does not have the ability to sample the full spectrum of possible BoB tropical cyclone events. Simulations driven by a wider range of initial conditions, derived from a wider range of historical 205 cases, would improve the sample size of cyclonic conditions on which this analysis is based. Note that this wouldn't necessarily reduce uncertainty in exceedance thresholds (in a frequentist paradigm), but it would update our view (i.e. our posterior estimate) of what is credible within the continuum of possible tropical cyclone events. In Bayesian parlance, our posterior view of Bangladesh tropical cyclones would become our new prior belief if subsequent simulation data became

available. 210
A different limitation is posed by the initial aggregation of the 4.4 km model over time. This removes our ability to draw inferences on annual occurrence of (or longer-term variability in) TC events. This means that our estimates of exceedance probabilities are conditional on a tropical cyclone event actually impacting Bangladesh. For the purposes of risk assessment, we do not feel this limitation is significantcurrent generation weather forecast models are capable of accurately predicting 215

Summary & Conclusions 220
Generalised additive models (GAMs) provide a useful framework for condensing spatial hazard information in an interpretable way, from multiple numerical model simulations, into a single spatially coherent hazard map. Using a restricted maximum likelihood approach to fit the GAM allows us to interpret model predictions in a Bayesian fashion that logically provides credible exceedance estimates. High-resolution convection-permitting numerical predictions of 12 historical cyclone events, in an ensemble model set-up, gives an improved sense of the plausibility and likelihood of possible extreme events without 225 being constrained by the lack of observational history in this region. Combining ensemble simulations with a GAM then allows us to robustly quantify the likelihood of maximum gust speed exceedances in a spatially coherent manner.
Our new maps of exceedance intervals show that north-western provinces of Bangladesh are relatively exposed to high wind speed hazardsin some areas the exceedance probabilities are equal to those experienced along the coast. Our hazard-to-230 decision making framework suggests that these areas may need to be considered in an equivalent manner to coastal regions, from a disaster risk reduction perspective. In coastal areas of Cox's Bazar and Chittagong we show super cyclonic conditions may occur as frequently as 1-in-20 to 1-in-100 years. We hope that these kilometre scale hazard maps facilitate one part of the risk assessment chain to improve local ability to make effective risk management and risk transfer decisions. Future work to co-produce a proper loss function, given wind speed thresholds, would facilitate a method of transparent operational decision 235 making that could be used as the basis of an operational warning system.

Code availability
Python, R and GAM model and data analysis code is available at https://doi.org/10.5281/zenodo.3953772

Data availability
The data used in this study is available at https://doi.org/10.5281/zenodo.3600201 and released under CC-BY 4.0 240

Author contribution
HS prepared the manuscript, with input from TE, and undertook the data analysis. TE and HS jointly developed and coded the GAM model. TE developed and coded the decision-making framework used in 3.1.

Figure 1
Summary of generalised additive modelling and the derivation of the posterior prediction gust speed distribution. The posterior predictive distribution is found for each grid cell of the regional model domain. Gust speed prediction intervals are found from the percentiles of the posterior predictive distribution.

Figure 2
Detrended quantile-quantile (worm) plots for each GAM model per storm. We discretise the quantiles into 50 bins (open circles). The red dashed lines represents zero deviance between data and theoretical quantiles defined in the GAM. Where model quantile deviates below (above) the zero deviance line, this implies that the model predictions are overestimated (underestimated) relative to the data: for any given theoretical model quantile, the data quantile is lower (higher). Deviance residuals respect the model family used when fitting the GAM and are calculated via the simulation method of Augustin et al. (2012).   (WMO, 2018). Event exceedance probabilities show the likelihood of a maximum tropical cyclone gust speed being greater than or equal to the corresponding classification wind threshold, conditional on a tropical cyclone making landfall over Bangladesh. An exceedance threshold of 50% (0.5%) represent a 1-in-2 (1-in-200) chance of a tropical 395 cyclone exceeding a given threshold. Areas where the exceedance probability is > 50% (< 0.5%) are shaded black (grey). Province boundaries are outlined in white, with the 18 most populated towns and cities marked by circles.   Table 1.

Loss Function
Warning Level (y)

OK Warning Disaster Great Danger
Event ( Table 1 Dummy loss function for actions associated with 4 Bangladesh TC warning levels, and their associated wind speed intensity. In this case loss is defined on a 100-point scale, where 0 = no loss, and 100 = maximum loss, associated with a given landfall TC event.