Towards an efficient storm surge and inundation forecasting system over the Bengal delta: Chasing the super-cyclone Amphan

The Bay of Bengal is a well-known breeding ground to some of the deadliest cyclones in history. Despite recent advancements, the complex morphology and hydrodynamics of this large delta and the associated modelling computational costs impede the storm surge forecasting in this highly vulnerable region. Here we present a proof of concept of a physically consistent and computationally efficient storm surge forecasting system tractable in real-time with limited resources. With a state-of-the-art wave-coupled hydrodynamic numerical modelling system, we forecast the recent super cyclone Amphan in 5 real-time. From the available observations, we assessed the quality of our modelling framework. We affirmed the evidence of the key ingredients needed for an efficient, real-time surge and inundation forecast along this active and complex coastal region. This article shows the proof of the maturity of our framework for operational implementation, which can particularly improve the quality of localized forecast for effective decision-making.

May) and another during post-monsoon (October-November) (Bhardwaj and Singh, 2019). The Bay of Bengal is a small semienclosed basin (Figure 1), which accounts for only 6% of cyclones worldwide. However, more than 70% of global casualties from the cyclones and associated storm surges over the last century occurred there (Ali, 1999). The number of fatalities concentrates in the Bengal delta across Bangladesh and India (Ali, 1999;Murty et al., 1986) where more than 150 million 25 people live below 5 m above mean sea level (MSL) (Alam and Dominey-Howes, 2014). Seo and Bakkensen (2017) noted a statistically significant correlation between storm surge height and the fatality in this region. Two of the notable cyclones that struck this delta include the 1970 cyclone Bhola and 1991 cyclone Gorky which killed about 300,000 (Frank and Husain, 1971) and 150,000 (Khalil, 1993) people respectively. In recent decades, the death toll has reduced by orders of magnitude. For example, cyclone Sidr, a category five equivalent cyclone on the Saffir-Simpson scale that made landfall in 2007, claimed 3406 30 lives (Paul, 2009). At the same time, the cost of material damages has significantly increased (Alam and Dominey-Howes, 2014;Emanuel, 2005;Schmidt et al., 2009). During recent cyclones, government and voluntary agencies of Bangladesh and India took coordinated effort to evacuate millions of people to safety at cyclone shelters before the cyclone landfall (Paul and Dutt, 2010). This kind of informed coordination benefited from improved communication, increased shelter infrastructure, and, most importantly, from the improvement of the numerical prediction of storm track and intensity. 35 Over the last decades, global weather and forecasting systems have advanced significantly. Several global models now run multiple times a day, starting from multiple initial conditions, with horizontal resolution in the order of tens of km, providing forecast with a range of hours to a week. These forecasts have proven useful during weather extremes like tropical storms (Magnusson et al., 2019). Operational hurricane forecasting systems like Hurricane Weather Research and Forecasting (HWRF) have emerged and reached a level of maturity to provide reliable cyclone forecasts several days in advance throughout the 40 tropics (Tallapragada et al., 2014).
Storm surge forecasting has been developed and implemented around the world along with the advancement of the weather and storm forecasting (Bernier and Thompson, 2015;Daniel et al., 2009;Flowerdew et al., 2012;Lane et al., 2009;Verlaan et al., 2005). Nowadays, operational surge forecasting systems typically run on high-performance computing systems, either on a scheduled basis or triggered on-demand during an event (Khalid and Ferreira, 2020). In the Bay of Bengal region, the  Yang et al., 2020), coarse-resolution modelling (Suh et al., 2015), and reduced-physics modelling (Murty et al., 2017). In the past decade, unstructured-grid modelling systems are getting more and more popular due to their efficiency in resolving the topographic features and their reduced computational cost compared to structured-grid equivalents (Ji et al., 2009;Lane et al., 2009;Melton et al., 2009). 60 The published history of storm surge modelling in the northern Bay of Bengal goes hand-to-hand with the land-falling of very severe cyclones (Das, 1972;Flather, 1994;Dube et al., 2004;Murty et al., 2014;Krien et al., 2017). From the modelling of the historical storm events, previous studies have identified several ingredients as essential for accurate modelling and forecasting of storm surge induced water level and associated flooding over the Bengal delta. The most important of these ingredients is an accurate bathymetry and topography (Krien et al., 2016(Krien et al., , 2017Murty et al., 1986). Second, it is required to have a large- Figure 2. Conceptual diagram of the involved processes that determine the water level evolution and its interaction with the controls determining the inland flooding. Idier et al., 2019). Due to this interaction, the highest surge is obtained for a storm making landfall around 2 hours before 100 the high tide. The strong dependence on tide reinforces the importance of having an accurate tidal model for the region and a reliable network of water level gauge for validation (As-Salek and Yasuda, 2001;Krien et al., 2017;Murty et al., 2016).
Wave setup is another component that has a significant impact on the nearshore sea level (Idier et al., 2019;Krien et al., 2017;Murty et al., 2014). It manifests as an increase in the sea level occurring in the nearshore zone that accompanies the dissipation of short waves by breaking (Longuet-Higgins and Stewart, 1962). During cyclone Sidr, the modelled wave-setup 105 was around 30 cm (Krien et al., 2017). This estimation would be potentially underestimated by up to 100% due to an early dissipation of wave energy by depth-induced breaking arising from usual parameterizations utilized in spectral wave models Pezerat et al. (2020).
At the seasonal scale, the mean sea level around the Bengal delta also shows considerable evolution. The amplitude of this variation can go as high as 40 cm due to freshwater influxes during monsoon from the GBM river system and offshore-ocean 110 steric variability Tazkia et al., 2017). During the wintertime cyclone-prone season, this steric variability can induce 10-15 cm variation of mean sea level of the bay.
Except for the Sundarban mangrove forest, coastal Bangladesh is mostly embanked through a network of 139 polders (Nowreen et al., 2013). These embanked areas can be flooded in three ways -(a) by overflowing/overtopping, (b) by breaching of the embankments/control structures, and (c) by heavy rainfall inside the embanked area. The height of these embankments 115 plays the most vital role in actuating the inundation from a storm surge (Krien et al., 2017). The performance of the embankments during the cyclone is another crucial factor, depending on the composition (typically, earth vs stones/concrete). During a cyclone, the wave action as well as overtopping/overflowing on the embankments can create breaches at a weak point, thus creating local inland flooding (Islam et al., 2011). Particularly for Bangladesh, the consequence of such storm surge induced flooding can be long-lasting as the topography inside polders is often below the mean sea level due to ground compaction 120 (Auerbach et al., 2015).

Atmospheric evolution of cyclone Amphan
On May 13th of 2020, a low-pressure area was persisting over northern Bay of Bengal about 300 km east of Sri Lanka. By the end of 15th, the Joint Typhoon Warning Center (JTWC) upgraded the low-pressure system to a tropical depression. Indian Meteorological Department (IMD) also reported the same development on the next day. The tropical depression continued to 125 move northward and intensified into the cyclonic storm Amphan by May 16th 18:00 UTC. During the following 12 hours, the intensification of the system was limited. However, starting from 12:00 UTC on May 17th, Amphan started to intensify very rapidly. In just 6 hours around 18:00 UTC, the maximum wind speed increased from 140 km/h to 215 km/h, making it an extremely severe cyclonic storm (equivalent to Category 4 on the Saffir-Simpson scale). Over the next twelve hours, Amphan continued to intensify reaching a maximum of 260 km/h wind-speed and 907 millibars central pressure, making it 130 the most intense event on record in the Bay of Bengal. During this time, it appeared to form two distinct eye-walls, which is typical of intense cyclones. Over the next 24 hours, Amphan lost most of its strength due to the eye-wall replacement cycle in the presence of moderate vertical wind shear (30-40 km/h). The system continued to decay due to easterly shear and midlevel dry air. The vertical wind shear remained moderate during this period. The cyclone made landfall between 08:00UTC and 10:00UTC around (88.35°N,21.65°E) (Figure 1), at mid-tide. At landfall, the reported central pressure amounted to 965 135 millibars with a maximum wind speed of 150-160 km/h. Afterwards, over its journey inland, the system further eroded and disappeared by May 21st. The black lines in Figure 3 illustrates the evolution of Amphan from JTWC best track. The extended range outlook by IMD published on May 7th predicted cyclogenesis to occur during May 8th to May 14th with low probability RSMC (2020). Global and regional models started to predict a significant storm to happen in the Bay of Bengal as early as May 12th, eight days before the cyclone landfall and four days before the actual formation of the tropical depression.
The formation of the low-pressure system triggered the operational HWRF system of NOAA (https://www.emc.ncep.noaa.gov/ gc_wmb/vxt/HWRF/) on May 14th. Similarly, the operational system of IMD (https://nwp.imd.gov.in/hwrf/IMDHWRFv3.6) was also triggered a day later. For the rest of the paper, we confine our analysis to the forecast disseminated by NOAA-HWRF only, relying on the Automated Tropical Cyclone Forecasting System (ATCF) text output (Miller et al., 1990). Figure   3 illustrates a few of the selected forecast cyclone tracks sequentially issued from NOAA-HWRF.

145
The forecasts illustrate the convergence of the location of the landfall. As early as May 17th, three days before the landfall, the forecast tracks had converged towards the observed track. One day later, on May 18th, the forecast error on the storm track location had reduced to around 50 km. The forecast wind speed and central pressure did not show as much accuracy as the forecast trajectory. The initial forecasts captured the magnitude of the intensification. However, they failed to forecast the rapid intensification that occurred during the 17-18th of May. On the other hand, the subsequent forecasts initiated during May 17th 150 and May 18th failed to capture the rapid weakening of the system that occurred over May 19th-20th, before landfall. Overall, in terms of wind speed, the forecast evolution was accurate (within 1 m/s) only at 24 hours range.  Zhang and Baptista (2008). This model solves the shallow-water equations using finite-element and finite-volume schemes in an unstructured grid that can combine triand quad-elements. The model is applicable in baroclinic as well as barotropic ocean circulation problems, for a broad range of spatial scales from the creek scale to the ocean basins (Ye et al., 2020;Zhang et al., 2020Zhang et al., , 2016a. It has already been shown to have an excellent performance in reproducing shallow water processes over the Bengal shoreline and elsewhere, including 160 the coastal tide (Krien et al., 2016), wave setup (Guérin et al., 2018) and storm surge flooding (Bertin et al., 2014;Krien et al., 2017;Fernández-Montblanc et al., 2019).

Model implementation
We have implemented our model on an updated bathymetry of Krien et al. (2016) with the additional inclusion of 77,000 digitized sounding points from a set of 34 nautical charts published by Bangladesh Navy . Our bathymet-165 ric dataset is a blend of two digitized sounding datasets in the nearshore zone -one from navigational charts produced by Bangladesh Navy, and another being a bathymetry of the Hooghly estuary provided by IWAI (Inland Waterways Authority To account for the effect of short waves on the mean circulation, Wind Wave Model III (WWMIII), a third-generation spectral wave model, is coupled online with SCHISM in our modelling framework (Roland et al., 2012). In our configuration, WWMIII solves wave action equation on the native SCHISM grid, using a fully implicit scheme. The source terms in our model comprise the energy input due to wind (Ardhuin et al., 2010), the non-linear interaction in deep and shallow water, 180 energy dissipation in deep and shallow water due to white capping and wave breaking, and energy dissipation due to bottom friction. We run WWMIII with a 12 directional and 12 frequency bins. Water level and current are exchanged among the two models every 30 minutes. Wave breaking is modelled, according to Battjes and Janssen (1978). As the nearshore region of Bengal delta has a very mild slope, the value for the breaking coefficient α, which controls the rate of dissipation, was reduced from 1 to 0.1 to avoid over-dissipation (Pezerat et al., 2020). WWMIII was also forced along the model southern open  190 Our model is forced over the whole domain by the astronomical tidal potential for the 12 main constituents (2N2, K1, K2, M2, Mu2, N2, Nu2, O1, P1, Q1, S2, T2). At the southern boundary, we have prescribed a boundary tidal water level from 26 harmonic constituents (M2, M3, M4, M6, M8, Mf, Mm, MN4, MS4, Msf, Mu2, N2, Nu2, O1, P1, Q1, R2, S1, S2, S4, Ssa, T2, K2, K1, J1, and 2N2) extracted from FES2012 global tide model (Carrère et al., 2013). At each of the upstream river open boundaries of Ganges, Brahmaputra, Hooghly, and Karnaphuli, we implemented a discharge boundary condition.

195
At Meghna and Rupnarayan river open boundaries, we implemented a radiating Flather boundary condition. From a year-long tidal simulation, we found that the updated version of the bathymetry performs 2-4 times better at key locations compared to the global tidal models -which includes FES (Carrère et al., 2013), GOT (Ray, 1999), and TPXO (Egbert and Erofeeva, 2002) (see Table A1).
We derived the atmospheric wind and pressure fields for SCHISM from a blending of analytical wind field inferred from a 200 parametric wind model (close to the cyclone track) and background atmospheric field from the Global Forecast System (GFS, https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs) reanalysis (further away), following Krien et al. (2017Krien et al. ( , 2018. We used the best-track data of JTWC provided at 6-hour time step as input for the analytical wind and pressure field. Here we choose the analytical wind profile of Emanuel and Rotunno (2011) and the analytical pressure field profile of Holland (1980). To be consistent with our forecast described in Section 5, the background wind field was gener-205 8 https://doi.org/10.5194/nhess-2020-340 Preprint. Discussion started: 30 October 2020 c Author(s) 2020. CC BY 4.0 License. ated incrementally from an accumulative merging of GFS forecasts for each 6-hour forecast window, with a 1-hour time step.
The analytical and background wind fields were first temporally interpolated every 15 minutes and overlaid on the background GFS fields using a distance-varying weighting coefficient to the cyclone centre to ensure a smooth transition. We took into account the asymmetry of the cyclonic wind field following Lin and Chavas (2012).

Predictive skills 210
With the above-described model setup, we ran our model to calculate the water level and sea state evolution during Amphan.
The simulation was done a-posteriori, once the cyclone had passed, and JTWC and GFS hindcasts being already available throughout the cyclone lifetime. Figure 4 shows comparisons between observed and modelled significant wave heights (SWH), peak wave period and water level. First, we compare the significant wave height derived from the altimetric estimate of Sentinel-3B processed by the Wave Thematic Assembly Center from Copernicus Marine Environment Monitoring Service (https:// 215 scihub.copernicus.eu/), which overpass on May 18th 2020 at 16:03GMT. For each altimetric data point inside our domain, we interpolated the model output using a nearest-neighbour approach, both spatially and temporally. This comparison shows that SWH is reasonably reproduced, with RMSE typically within 1m (within 18% of the mean value). However, we observed an overestimation of SWH of about 2.5 m near the cyclone centre. One of the reasons for such difference is a limitation of analytical wind field models, as discussed by Krien et al. (2017). 220 We present an in situ observed time-series of SWH and mean wave period during the life cycle of Amphan in Figure 4 for BD08 (Figure 4, b-c) and BD11 (Figure 4, d-e) oceanographic buoys (the dataset is available online at https://incois. gov.in/portal/datainfo/mb.jsp). The two buoys are located on either side of Amphan track, 250 km to its east / 230 km to its west, respectively. We can see that at both locations, the model reproduces the SWH well. The match is particularly good at BD11, consistently within 1 m of the observed evolution. At BD08, SWH is found to be overestimated by about 2 m at its 225 peak, with a 12-h lead shift in time. Such difference can come from positional and timing uncertainties in JTWC best track data, from our assumption of a constant translational velocity in-between each 6-hourly JTWC time step, and finally from the intrinsic uncertainty of the analytical wind field itself, as shown on Figure 3b. Despite multiple sources of potential errors and uncertainties, our model appears capable of capturing the general features of the sea state, both significant wave height (SWH) and period, relatively well. The overall predictive skills of the model in terms of short waves are only slightly below those from 230 hindcast exercises typically conducted over other ocean basins using similar modelling strategies (Bertin et al., 2014;Krien et al., 2018).
The real-time availability of observed water level time series in this region is relatively scanty. We were able to access only a handful of coastal water level records, most of them located eastward of the landfall location. The right panels of Figure   4  kept unchanged. We also plotted the reconstructed tidal water level from the tidal atlas derived from our model. The overall 240 tidal propagation is well simulated in Galachipa and Kuakata, while at Angtihara and Tajumuddin the simulation tidal range is found to be underestimated. The local bathymetric error and friction parameterization might be the source of the discrepancy.
The tide gauge at Angtihara is located in a data-scarce location inside Sundarbans mangrove forest. Tide gauge at Tajumuddin is situated at the mouth of Meghna, where bathymetry is rapidly changing . However, our model correctly reproduces the peak of the water level at all locations. Depending on the site, the landfall occurred 2-4 h after high tide. The 245 10 https://doi.org/10.5194/nhess-2020-340 Preprint. Discussion started: 30 October 2020 c Author(s) 2020. CC BY 4.0 License. maximum recorded water level at these stations is around 2-2.5 m above mean sea level, with a storm surge (water level -tide) of the same order of magnitude.
To spatially assess the storm surge generated by cyclone Amphan, we first take a look into the temporal maximum water level in Figure 5(a). We can see that the whole coastal region experienced a high-water level ranging from 2-5 m for Amphan.
To quantify the contribution of the cyclone, we have looked in the surge component, defined as the difference between the total 250 modelled water level and a tide-only simulation. Figure 5(b) illustrates the temporal maximum of surge over the delta region.
The highest surge level is 5 m around the location of the landfall (88.3°E, 21.6°N). The maximum non-linear interaction between tide and surge amounts to about 30 cm in the nearshore domain (Figure 5d).
To quantify the contribution of the wave setup, we re-calculated the storm surge without the coupling with the wave model. Nevertheless, the spatial resolution at the nearshore region employed in this study (250 m) is still too coarse for capturing the maximum wave setup that develops along the shoreline, where a resolution of a few tens of m should be employed (Guérin et al., 2018). Nonetheless, this comparison shows that wave setup is not only developed along the shoreline exposed to waves but affects the whole delta up to far upstream, a process already observed over large estuaries (Bertin et al., 2015;Fortunato 260 et al., 2017).
Our findings reaffirm the necessity of a proper coupling between tide, surge and wave to forecast the water level evolution.
For Amphan, the magnitude of tide-surge interaction amounts to about 10% of the maximum total water level with spatial variation. Similarly, the magnitude of wave-setup is also spatially varying with a typical amplitude of 10-15% of the maximum total water level. This non-negligible non-linear dependency shows the importance of a fully coupled hydrodynamic-wave 265 modelling system for a successful forecast of the water levels in the hydrodynamic setting of the Bengal delta. Our modelling framework showed reasonable accuracy in reproducing observed storm surge water level along the Bengal coastline. The accuracy is in line with the typical level of performance of similar systems applied elsewhere in the world ocean (Bertin et al., 2015;Fernández-Montblanc et al., 2019;Suh et al., 2015). The hindcast experiment thus shows the viability of our model for a proof-of-concept of real-time forecasting exercise. The current state of maturity of atmospheric real-time cyclone forecast products allows their application to real-time surge forecasting, as discussed in Section 3. As explained in Section 2, the storm surge generated from cyclones depends on the atmospheric pressure, wind, and background (typically: tidal) water level. A realistic surge forecast issued three, or even two days before landfall can usefully predict the order of magnitude of the expected maximum water level. Particularly a spatial 275 prediction of expected maximum water level could be of utmost utility to determine the locations of impact and associated damages. This two to three-day timeframe can help to prepare the embankments or critical systems over the vulnerable areas, and prepare for the evacuation of the population living in low-lying coastal zones. Storm surge and water level evolution Figure 5. Hindcast of (a) maximum water level, (b) maximum surge (c) wave setup/setdown (d) maximum non-linear interaction between tide and surge. For (a), for the areas above mean sea level, the water level is converted to water level above the topography for consistency.
The inset maps show a close-up (75 km × 45 km) of the landfall region. estimated 24 hours to 12 hours before landfall, once the cyclone forecasts from the operational agencies have converged, can further help to elaborate location-specific evacuation orders.

Forecast strategy and forcing products
During cyclone Amphan, we performed near-realtime storm surge forecasts based on our model. We communicated the results to Bangladesh local government authority through personal communications, as well as to the scientific community through social media. In our forecasts, we relied on the outputs of global models (GFS, and HWRF) and advisories (JTWC) for the estimates of the atmospheric forcing. These advisories and global models are updated every 6 hours. Similar to the GFS forecast 285 cycle, we updated our forecasts at every 00:00, 06:00, 12:00, 18:00 UTC.
For each forecast, we derived the atmospheric forcing from a blend of JTWC, HWRF, and GFS data (Table A2). JTWC publishes storm advisories at 03:00, 09:00, 15:00, 21:00 UTC. The advisory includes an analysis of the storm intensity and position from a satellite fix 3-hours prior, followed by a forecast for the following 72 hours, with 6 to 12 hours time-steps.
NOAA published their HWRF model advisory for each GFS forecast cycle at 3-hour time steps. GFS provides the wind and 290 pressure fields in hourly timesteps at 0.25°regular spatial grid. The GFS model is initialized 6-hourly at 00:00, 06:00, 12:00, 18:00 UTC, and the data is available 3 hours after the initialized period.
The wind and pressure fields around the centre of the storm are derived analytically from a concatenated JTWC and HWRF.
We used GFS data as the background wind and pressure on the outer region of the analytical fields.
We initiated our forecast cycles on March 16th 2020 at 06:00 UTC -utilizing the first advisory from JTWC published at 295 03:00 UTC, merged with the forecast from HWRF issued at 18:00 UTC of the preceding day. We took the background wind and pressure field from GFS forecast published at 00:00 UTC. In the subsequent forecast cycles, we updated the previous track file by first replacing the past time-steps with the analysis from the latest JTWC advisory, then appending the forecasts from HWRF again. For the background wind and pressure fields from GFS, we retained the first 6 hours of forecasts from the previous cycles and updated the remainder of the time series from the latest forecast with new initialization. For each forecast 300 cycle, we kept the starting time of our model the same, on May 16th at 00:00 UTC, and ran till one-day after the landfall for a consistent initial condition. This approach is feasible as the storms in the Bay of Bengal typically form, grow and dissipate in about a week. A graphical representation of our strategy to temporally concatenate JTWC, GFS, and HWRF forecasts over a forecast cycle is shown in Figure 6. Number of MPI processes 20 Wall clock time (pre-processing) 10min Wall clock time (model integration) 1h45min Wall clock time (post-processing) 5min Disk storage (each simulation) 8GB Memory usage 18GB

305
Since our forecast is contingent to the JTWC advisory of three hours before, we had only a three-hour-long window to achieve the pre-processing of the forcing files, to run the model, and to complete the necessary post-processing of the model outputs.
The numerical efficiency of our model made it possible to cope with this time constraint of the forecast, using a modest computing resource. The computing environment we used is detailed in Table 1, which essentially amounts to a single desktop computer fitted with a high-end consumer-grade processor. 310 Figure 7 illustrates the evolution of the maximum total water level and surge for the corresponding forecasts issued at T-60 hours, T-36 hours, and T-12 hours, where T is the time of landfall (May 20th 2020 at 08:00 to 10:00 UTC). For the forecast issued at T-60 hours, the cyclone landfall is located near 90°E, associated with the strong surge simulated over the surrounding area. As the atmospheric cyclone forecasts gradually converge towards the actual landfall location more than 100 km westward, so do our storm surge forecasts at T-36 hours and T-12 hours. The dependency of the timing of the landfall to the forecast range 315 is notable here. At T-36 hours, the maximum water level pattern is similar to the hindcast shown in Figure 5. By T-12 hours, the magnitude of the forecast maximum water level corresponds well with the hindcast estimate.

Coupled tide-surge-wave dynamics
The water level dynamics over the Bay of Bengal and particularly around the Bengal delta is complex, and the various compo-320 nents of water level have different relative contributions depending on the location. Due to tide-surge interaction, the maximum water level can be lower or higher depending on the tidal phase in a fully non-linear tide-surge coupled simulation, compared to its linear counterpart. During Amphan the tide-surge interaction was mostly negative, except within about 10 km of the landfall location. The positive interaction near the landfall location means that the coupled tide-surge estimate was higher than the linearly added tide and surge. In this case, the typical magnitude of tide-surge interaction ranged from -30cm to +30cm.

325
The contribution from the wave-induced setup along the shoreline increases the maximum water level estimation throughout the coast by about 10-15% during Amphan. Fully coupled tide-surge-wave models have been recommended for operational forecasting in the Bay of Bengal for years (Johns and Ali, 1980;Murty et al., 1986;As-Salek and Yasuda, 2001;Deb and Ferreira, 2016;Krien et al., 2017). Besides the development of the wave-ocean coupled hydrodynamic modelling tool itself, the main challenge in implementing such a system sums up to the high computational requirements for the proposed modelling 330 systems (Murty et al., 2014(Murty et al., , 2016. Due to this constraints, operational systems are commonly run with either only surge, or only coupled tide-surge model, but without waves (Murty et al., 2017) and with a much coarser spatial resolution.

Performance of storm surge forecasting
Aside from the necessary inclusion of relevant physical processes -tide, surge, waves and their non-linear interactions -the performance of storm surge forecasting depends on the wind and pressure forcings. The errors in the wind forcing and in the 335 atmospheric pressure forcing are typically the most important ones among the various sources of error in storm surge modelling (Krien et al., 2017). In our forecast environment, the source of this error is two-fold -the analytical model used to synthesize the fields from the storm parameters, and the numerical weather model used to forecast the storm wind and pressure fields to derive the storm parameters themselves. There have been attempts to find a more accurate formulation from recent satellite scatterometer missions (Krien et al., 2018). However, not one formulation works best in all radial distances from the storm 340 centre. Despite the well-known error associated with the parametric wind field, they are widely used in applications due to their computational simplicity and lightweight data requirements (Lin and Chavas, 2012). The best way to avoid the error from the analytical wind field might be not using these formulations and rely on the full-fledged atmospheric forecasts. However, atmospheric models are costly in terms of computation. It is particularly true for cyclonic storms where high spatial resolution (typically kilometric) is required (Tallapragada et al., 2014).

345
From the validation of water level at a few tide gauge locations for our hindcast simulation, it appears that the modelling framework proposed here could capture the maximum surge level successfully throughout the coast. Our tide-gauge sites were limited to the eastern side of Amphan landfall location and were operational throughout the life cycle of the cyclone. However, it was not possible to confirm the ability of the model to reproduce the water level on the western side of the landfall location due to data unavailability. The proper reproduction of the maximum water level by the model gives strong confidence in the 350 formulation as well as in the coupling strategy of our modelling framework.

Inundation
One of the potential, but challenging, outcomes of storm surge forecasting is the prediction of inundation. Particularly in the northern Bay of Bengal, the existing modelling studies do not generally consider the inundation process (Dube et al., 2004;Murty et al., 2014). Some modelling systems take advantage of simplified inundation modelling schemes to tackle this problem 355 (Lewis et al., 2013). While the recent improvements in models bathymetry and the detailed accounting of dikes and coastal defences improved the overall modelling of the inundation, the associated error remains large (Krien et al., 2017). As an update to Krien et al. (2017), in our model we have incorporated a novel dike heights dataset, bearing varying dike crests heights, for the numerous polders scattered around the Bengal delta. The assessment of the impact of the updated embankment heights is, however, hard to quantify and validate.

360
In the absence of any operational network of inundation monitoring, to understand and better characterize the inundation dynamics, we systematically skimmed the Bengal delta local newspapers on the day of Amphan landfall and on the few subsequent days, to achieve an as-comprehensive-as-possible mapping of the inundation extent observed in situ. We digitized the reported flooded locations through © Google Earth and categorized the inundation mechanism. We digitized 88 such locations over the Bangladeshi part of Bengal delta, as shown in Figure 8 (Khan, 2020). Over the Indian side of the delta, the news reporting of inland flooding from dike breaching was not accurate enough for us to be able to geotag the locations. We overlaid the inundated locations over the inundation predicted by the model on a false-colour image derived from Sentinel-2 satellite using © Google Earth Engine (https://earthengine.google.com). Three categories of inundation mechanism are typically observed -inundation by breaching of the dikes (labelled as "Embankment breach"), inundation of unprotected low land by increased water level ("Unembanked lowland"), and flooding by overflowing of the dikes ("Embankment overflow").

370
From the reported news, it is clear that the major part of the inundation results from the breaching of dikes. We also note two breaching instances in the south-eastern corner of our domain, very far from the cyclone landfall. On the other hand, the flooding from merely increased water level in non-diked areas is mostly concentrated along the Meghna estuary. Along the estuary, designed embankments probably do not exist to protect many populated low-lying char areas (land/island formed though accretion inundated area, which is known not to be polderized. While we could find some reports of flooding around that region from our newspaper survey, the widespread inundation predicted by our model in this region probably did not occur in reality due to the presence of levees and city protection embankments. Indeed, this kind of small-scale geographic information, to the best of our knowledge, is not systematically available publicly and thus could not be incorporated in our model topography, despite being 380 of utmost importance for localized inundation forecast. The forecast flooded vs non-flooded areas show a wealth of spatial scales, demonstrating that our unstructured modelling framework has the intrinsic capability to model the small-scale hydrodynamic gradients. In our modelled forecasts, the contrast between the well-predicted water level temporal evolution, and the comparatively limited predictive skills of inundated locations points out the necessity of accurate knowledge and incorporation of reliable topographic and embankment information, at the local scale. Our modelling system has the ability and appropriate 385 resolution to ingest such topographic information efficiently due to its unstructured nature. The coastal polder management has long been a boon and a bane (Warner et al., 2018). On one side, the polders have protected the fertile areas from daily tidal flooding of brackish waters, which favoured agricultural production by reducing the soil salinity (Nowreen et al., 2013). On the other hand, these embankments have restricted freshwater and sediment supply, causing significant subsidence inside the embanked areas due to sediment compaction (Auerbach et al., 2015). Over time, the 390 land use inside these polders has also changed from agriculture to aquaculture (namely shrimp and crab farming). Although the embankments were not built to protect from the cyclonic surges, in many instances they functioned as such. The presence of these protective structures also contributed to a false sense of security in the residents living inside (Paul and Dutt, 2010).
Embankments have become part and parcel of survival in this low-lying region. It has now become essential to monitor the condition and topography of the embankments for efficient management and maintenance. From a forecasting point of view, 395 embankment breaching is extremely hard to model in the state-of-the-art hydrodynamic modelling frameworks. However, consistent periodic monitoring of the embankment conditions can improve the forecast by providing a more objective view of the associated inundation hazard.

Prospects in operational forecasting
Throughout our forecast and hindcast experiments, we have relied on freely available forecast products as the forcing for our 400 storm surge model. These forecasts data are publicly available within 3 hours of their initialization time through online portals.
The primary input to generate the forcing is the forecast storm track distributed as a text file, which is very lightweight even for a limited internet connection. We have implemented our storm surge model on a freely available open-source framework. The model code is already operationally used in regional forecasting systems and services (Fortunato et al., 2017;Oliveira et al., 2020).

405
As we have demonstrated, the modelling setup presented here is tractable in an operational, real-time forecasting scenario.
Thanks to advanced and efficient numerics, one can deploy such a system in a consumer-grade computing environment. The computing setup we have deliberately used (see Table 1) essentially amounts to a high-end x86-64 desktop personal computer (e.g. Intel® Core™i9-10900 or similar). This computing requirement takes only a portion of the available computing capacity of an institute such as the BMD (15x4 cores@2.8GHz) . Furthermore, the model is easily deployable in cloud 410 computing infrastructure (Oliveira et al., 2020), which could provide additional reliability and cost-savings for an event-driven operational forecasting system.

Conclusions
In this study, we present a retrospective evaluation of storm surge prediction during super-cyclone Amphan using an efficient, state-of-the-art coupled hydrodynamic-wave numerical model. During cyclone Amphan, the predicted maximum water level From our proactive forecast initiative during cyclone Amphan, we showed that with publicly available storm forecast products and easily accessible computing resources, it is feasible to forecast the evolution of water level throughout the vast coastline of the Bengal delta in real-time. During Amphan, a sufficiently skilful storm surge forecast was achieved as early as 36-48 hours before landfall. The forecast water level with 36-hour lead time seemed quite similar to our best (hindcast) simulation, in terms of maximum water level as well as the spatial pattern.

425
From a secondary post-disaster news-survey, we have identified the limiting factors for location-specific inundation forecasts. The main limiting factor is due to limited and outdated topographic information of the existing coastal defences. Also, dike-breaching, which was the prevalent process of inland inundation during Amphan, is not explicitly modelled in our hydrodynamic framework. These two factors call for routine monitoring of embankment topography and condition.
Cyclone and storm surge warning has always been a communication and trust issue in the Bay of Bengal (Paul and Dutt,430 2010; Roy et al., 2015). It is thus necessary to communicate well-grounded storm and surge forecast to the community for coordinated and informed decision-making during a storm surge (Morss et al., 2018). The forecasting system we implemented and assessed in the present case study provides the proof-of-feasibility and opens short-term operational prospects to fill a gap in the existing disaster management tools in this part of the world.
Code and data availability. The instructions to download and install the model used in this study can be accessed freely at https://github.com/schismdev/schism. The sources of the data used in this study have been described in the article. The data should be requested through the mentioned institutions or downloaded from the provided websites.  (Mayet et al., 2013). The harmonic analysis is done using the Tidal toolbox developed at LEGOS (Allain, 2016). The modulus of the complex difference defines the complex error for a tidal constituent.
(m) and observation (o). The total error of all the constituent at one location is calculated as the squared root of half of the squared sum.
Along the coast of Bengal delta, only four of the constituents -M2, S2, K1, and O1 are found to contribute significantly to the tidal energy (Sindhu and Unnikrishnan, 2013). As in many cases, information for other tidal harmonics is not available, 450 only these four constituents are considered for calculating the total complex error at a location.
A comparison of the complex error between the global models and the model presented here is shown in Table A1. Amplitudes (A) and errors are in centimeter, phase (φ) is in degrees. Hooghly River, Diamond Harbour, Garden Reach and Chandpur are not represented in global tidal models (FES, GOT, and TPXO) due to their location in far upstream.
A2 Forcing data source 455 The source of the forcing data used in this study for hindcast and forecast is listed in Table A2 20 https://doi.org/10.5194/nhess-2020-340 Preprint. Discussion started: 30 October 2020 c Author(s) 2020. CC BY 4.0 License.