Wave Boundary Layer Model based wind drag estimation for tropical storm surge modelling in the Bay of Bengal

Ocean Engineering 191 (2019) 106509

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Wave Boundary Layer Model based wind drag estimation for tropical storm surge modelling in the Bay of Bengal C. Gowri Shankar, Manasa Ranjan Behera * Department of Civil Engineering, Indian Institute of Technology Bombay, Mumbai, Maharashtra, India

A R T I C L E I N F O

A B S T R A C T

Keywords: Storm surge Storm wave Wind drag coefficient Wave Boundary Layer Model ADCIRC-SWAN Bay of Bengal

Modelling of storm surge numerically has always been strenuous as it is associated with various influencing factors. The major governing component being the wind forcing or the wind stress - that signifies the compu­ tational accuracy of simulated surge and wave parameters. The present study is aimed at analysing the most reliable wind drag evaluation method for real-time predictions of storm surge. In an attempt to examine the same, two tropical cyclones (Thane and Vardah) eventuated in the Northern Indian Ocean were modelled. Four most widely used linear empirical and quadratic relations along with the enhanced Wave Boundary Layer Model (WBLM) were executed for the estimations of wind drag coefficient. The surge was subsequently simulated (using the coupled hydrodynamic circulation and wave model: ADCIRC and SWAN, respectively), individually for each of the above wind stress methods to obtain the corresponding water levels and wave features. The modelled values were further validated with the observed water level and wave data. It was quite intuitively observed that, WBLM outperformed (i.e. significantly correlated with the observed values) its linear counterparts since, the former pragmatically includes the effects of air-sea interface at high wind speeds in the model. However, due to insufficiency of the observed parameters, water levels obtained at locations of maximum wind speeds could not be satisfactorily validated. The WBLM-based computation of significant wave heights in deep as well as shallow water nevertheless enabled efficient and reasonably-reliable estimations of the peak incidents (along with relative wave directions and energy densities). Although this study essentially includes cyclonic wind velocities, it was identified that - with given wind field that includes the ‘pre-existing’ wind-wave climate in the relevant domain (i.e. before the cyclone genesis) the modelled estimates could be more precise.

1. Introduction Cyclonic storms have proved catastrophic in multiple occurrences, thereby mandating immediate and dedicated study for precise real-time predictions. It has been observed that five out of the seven SAARC Na­ tions (except Nepal and Bhutan that are surrounded by land mass) are greatly affected due to tropical cyclones, tidal waves and their augmented effects. The boundless coastline of India (�7500 km) sur­ rounded by the two major basins of the Northern Indian Ocean (Bay of Bengal on the east and Arabian Sea on the west) has historically always been prone to severe cyclonic storm surge events. In particular, the Bay of Bengal (BoB) basin has emanated the most devastating ones that has been cataclysmic in varying degrees. Components like the shallow shores of the BoB, the low and flat coastal terrain of India and Bangladesh and the funnel shaped coastline combination tend to pro­ duce such intensive surge events. Almost 14% of the population of India

have their habitats along the coastal stretches and are highly vulnerable to such extreme sea-level events. In the recorded history of such events, the surge (~10 m) caused by Cyclone Bhola in the year 1970 on the coasts of Bangladesh (Madsen and Jakobsen, 2004) is considered the most intensive one. A Super Cyclonic Storm (SCS) that made its landfall in the Odisha coast in 1999 is observed to be the highest (of about 7.5 m) surge causing event in context of Indian scenario (Dube et al., 2000). Thus, it necessitates a quantitative understanding of these phenomena, which could consequently aid the prediction, such that the precaution­ ary actions could be adopted on time. Although, quantifying a natural phenomenon includes numerous complicated factors per se, several models have been developed to simulate and predict a storm’s behav­ iour. Storm surge modelling and prevising of such events therefore, need to be essentiated in order to safe-guard the lives and interest of the nation. Numerically modelling storm surge has always been considered an

* Corresponding author. E-mail address: [email protected] (M.R. Behera). https://doi.org/10.1016/j.oceaneng.2019.106509 Received 1 April 2019; Received in revised form 29 August 2019; Accepted 28 September 2019 Available online 7 October 2019 0029-8018/© 2019 Elsevier Ltd. All rights reserved.

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

Fig. 1. Bay of Bengal basin with the tracks of Thane and Vardah (blue and red line). The locations of the tide gauge stations and buoys are demarcated with circle and star respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. Comparison of wind drag coefficient estimation with different formulations at various wind speeds.

exacting task and the research community has constantly been working towards better and more precise solutions over the past few decades. Such modelling requires various input parameters like wind velocity, central pressure of the cyclone, local bathymetry, tidal input and wave stress gradients (in case when coupled with a wave model). These input parameters are collected and generated through various data sources and models, which is then fed to the surge model to simulate the vertical elevation of surge and inundation levels. A detailed approach on such studies towards Indian coastal regions are reported by Murty et al. (1986); Dube et al. (1997, 2009). It is a familiar fact that the ocean currents and waves are steered through the momentum fluxes across

air-sea interface. To evaluate this momentum flux at the air-sea interface at high wind speeds is still an exhausting way to accomplish through numerical methods. A simplified approach to identifying this transfer of momentum is through shear stress, which is also known as the wind stress and is expressed as

τ ¼ ρa u2* ¼ ρa Cd u210

(1)

where, ρa is the density of air, u� is the shear velocity, u10 is the hori­ zontal wind speed at an height of 10 m above the sea surface and Cd is the wind drag coefficient. Several researchers have stated about the empirical computation of 2

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

reduces at very high winds (crossing 35–40 m/s). Hence the un­ certainties revolving around this wind drag coefficient has not been completely perceived. The phenomenon of wind drag reduction at high winds has report­ edly (Makin and Kudryavtsev, 2002; Donelan et al., 2004; Kudryavtsev, 2006) been ascribed to the change in momentum transfer due to the uninterrupted breaking of waves. However, some other studies (Andreas, 2004; Kudryavtsev and Makin, 2011, etc.) assert that it is due to the sea-spray effect on the boundary layer that was consequently generated due to the same phenomenon (i.e. uninterrupted breaking of waves). Therefore, the preferable alternative is to cap the wind drag on a threshold value or even assume a law for decreasing coefficients (Bunya et al., 2010; Dietrich et al., 2010; Zijlema et al., 2012). It is hence noteworthy to mention that the study performed by Moon et al. (2009) with the WBLM coupled with wind and surge models, generated better results with realistic scenario. Further, the development of WBLM through the inclusion of sea-spray effect (Chen and Yu, 2016, 2017) has brought in a pragmatic way of representing the wind drag coefficient for extreme wind conditions (>35 m/s). The recent, enhanced WBLM (Chen and Yu, 2016) solves the dynamics through momentum and energy conservation equations of the wave boundary layer, by accounting for more physical attributes involved in an air-sea interaction process. The total stress computed through the model is a constant and is conven­ tionally said to be the sum of the wave induced stress and turbulent stress. With the above mentioned way of including the sea-spray effect, the WBLM is able to reproduce the decreasing tendency of the wind drag coefficient at very high wind speeds (Chen and Yu, 2016). The storm wave modelling study by Chen and Yu (2017) also substantiated the significance of the above-mentioned enhanced model. However, its ef­ ficacy in assessing surge caused by a tropical cyclone is yet to be authenticated. Also, Chen and Yu (2017) executed the enhanced WBLM with the wave spectrum obtained from a standalone run of a wave model (SWAVE) where the ambient flows between the tide, surge and wave is neglected. Hence, the current study directs its objective towards a comparative analysis of cyclonic surge modelling using the enhanced WBLM, with the two dimensional wave spectral input coming from a coupled circulation and wave model to account for the non-linear interaction of tide, surge and wave. Along with the enhanced WBLM, four other bulk formulae methods are also used for computing the wind drag coefficients followed by the storm surge and wave estimations. For executing the same, two cyclonic conditions namely, Cyclone Thane (2011) and Cyclone Vardah (2016) (as their tracks and landfall points are shown in Fig. 1) were considered. Both these cyclones have been characterised as Very Severe Cyclonic Storm (VSCS) by the Indian Meteorological Department (IMD), as their wind speed exceeded 35 m/s. The present study employs the coupled hydrodynamic circulation model ADCIRC and the third gener­ ation wave model, Simulating Waves Nearshore (SWAN, e.g., Booij et al., 1999) to simulate the storm surge with the wind stress coefficients obtained from WBLM. The following sections give a brief discussion about both Cyclone Thane and Cyclone Vardah, and the wind fields developed for the modelling. The execution of the coupled model (ADCIRC þ SWAN) with the empirical methods and WBLM. The results obtained from numerical modelling were subsequently validated with the observed data (tide gauge stations and buoys) collected from the Indian National Centre for Ocean Information Services (INCOIS). Further, a comparative analysis was carried out in the evaluations of storm surge and storm wave pa­ rameters in correlation with the wind drag estimation methods.

Fig. 3. Methodology implemented for the simulation of the coupled model (ADCIRC þ SWAN) through WBLM method.

the wind drag coefficient as a function of wind speed. Nevertheless, this approach might not be suitable for extreme events like cyclonic surge estimation, as the wind speed exceeds 35 m/s quite frequently. The empirical relations proposed by Garratt (1977); Large and Pond (1981) and Wu (1982) are the most extensively used ones. The wind drag co­ efficient (Cd) is also used to calculate the boundary layer parameters as reported through many contemporary models like Sea, Lake, and Overland Surges from Hurricanes (SLOSH, e.g., Jelesnianski et al., 1992; Phan et al., 2013), Advanced Circulation Multidimensional Hydrody­ namic Model (ADCIRC, e.g., Westerink and Luettich, 1991; Dietrich et al., 2011; Bhaskaran et al., 2014), Finite Volume Coastal Ocean Model (FVCOM, e.g., Huang et al., 2010; Zheng et al., 2013). In terms of application, these empirical methods were usually presumed to be applied for high wind conditions, whereas the recent studies by Powell et al. (2003) and Holthuijsen et al. (2012) stated that it was inappro­ priate. The wind drag coefficient attains a levelled-off nature and rather

2. Cyclone description 2.1. Cyclone Thane On the afternoon hours (1200 UTC) of 25th December 2011 IMD reported a ‘Depression’ formed over the southeast Bay of Bengal and lay 3

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

Fig. 4. Model domain of Bay of Bengal basin with bathymetric depth in meters.

centred approximately 8.5� N; 88.5� E (IMD National Bulletin Report, 2011). The system intensified into ‘Deep Depression’ after 15 h with wind speed exceeding 14 m/s and transforming into a ‘Cyclonic Storm’ (CS - wind speed crossed 17 m/s) at 1800 UTC on 26th Dec 2011. The system was named as ‘Thane’ by IMD and it continued to stay as a cyclonic storm for the next 30 h period and moved towards the Tamil Nadu coast of India. It intensified further into a Severe and then a Very Severe Cyclonic Storm in 6 h interval at 1200 UTC on 28th Dec 2011. Thane attained its maximum intensity of 46 m/s and made its landfall around 0100 to 0200 UTC on 30th Dec 2011 between the coasts of Cuddalore and Puducherry along the east coast of India. The havoc caused by Thane to the coastal structures and properties were enormous. It was observed that a maximum of 1–1.5 m height of storm surges above the astronomical tide inundated the coastal stretches of Ennore, Chen­ nai, Puducherry and Cuddalore. It was reported that there were 46 fa­ talities together in Puducherry and Cuddalore district (IMD/WMO report, 2012).

3. Materials and methods

2.2. Cyclone Vardah

∂2 ζ ∂ζ ∂~J x ∂~J y þ τ 0 þ Sp þ ∂t2 ∂t ∂x ∂y

3.1. Advanced circulation (ADCIRC) model This study uses the two dimensional depth integrated (2DDI) ADCIRC model, that solves the shallow water equations (SWE) to obtain the water levels (ζ) and vertically integrated momentum equations to ⇀

obtain the currents (U) (Kolar et al., 1994; Luettich and Westerink, 2004; Dawson et al., 2006; Westerink et al., 2008). The uninterrupted Galerkin finite-element method is used to solve the unstructured triangular ele­ ments in space, whereas finite difference method is adopted in solving the temporal progression (Dietrich et al., 2011). This numerical tech­ nique maintains the stability in the model at high-resolution spatial discretisation where the gradients tends to be the largest (Tanaka et al., 2011). The water levels are obtained through GWCE (Generalised Wave Continuity Equation), which is the merged and differentiated form of continuity and momentum equations and is given as:

Cyclone Vardah formed as a depression in the morning hours (0900 UTC) of 06th December 2016, over the southeast Bay of Bengal, and lay centred approximately 8.5� N; 91� E (IMD National Bulletin Report, 2016). The system intensified into a ‘Deep Depression’ after 33 h and transformed into a ‘Cyclonic Storm’ at 0030 UTC on 08th December 2016. Further, it was transformed into a Very Severe Cyclonic Storm at 1200 UTC on 10th December 2016. Vardah attained its maximum in­ tensity of 38 m/s and made its landfall around 0900 to 1130 UTC on 12th December 2016 near Chennai. Vardah reduced to a Severe Cyclonic Storm (SCS) just before making landfall, and moved westward by dissipating its intensity into a Depression and finally as a Low Pressure system. Storm surges of maximum 1 m height inundated the coastal stretches of Ennore, Chennai and Thiruvallur districts of Tamil Nadu and Nellore district of Andhra Pradesh (IMD/WMO report, 2017). The tracks of cyclone Thane and Vardah has been depicted in Fig. 1 with demar­ cated locations of the tide gauge stations (Ennore and Chennai) and buoys (BD11 and Wave Rider Buoy Pondicherry (WRBP)). The tide gauge stations and buoy locations were chosen based on the path of both the cyclones (Thane and Vardah) and the data availability.

Sp UH

∂τ0 ∂x

VH

∂τ0 ¼0 ∂y

(2)

where: ~ Jx ¼ þ

Sp Qx

τsx;

wind

∂U ∂x

∂U þ fQy ∂y τbx waves

Qy

þ τsy;

g ∂2 ζ Sp 2 ∂x

þ ðMx

ρ0

�

�

g ​ Sp H

Dx Þ þ U

∂ Ps ∂x gρ0

∂ζ þ τ 0 Qx ∂t

αη g ​ Sp H

∂ζ ∂x (3)

and: ~ Jy ¼ þ

Sp Qx

τsy;

wind

∂V ∂x

∂V þ fQx ∂y τby waves

Qy

þ τsy;

ρ0

g ∂2 ζ 2 ∂y

þ My

�

�

∂ Ps αη ∂ y g ρ0 � ∂ζ Dy þ V þ τ 0 Qy ∂t gH

gH

∂ζ ∂y

(4)

The currents are obtained from the vertically integrated momentum equations and are given as follows:

4

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

Fig. 5. Comparison of model (GAHM) computed wind speed (m/s) versus the observed buoy (BD11) data for both the cyclones. (a) Cyclone Thane; (b) Cyclone Vardah.

∂U ∂U ∂U þ Sp U þV ∂t ∂x ∂y

∂ Ps αη ζþ ∂x gρ0 τbx Mx Dx wind þ τsy; waves þ ρ0 H H

τsx;

�

fU ¼

factor and y0 is the reference latitude; U and V are depth integrated currents in the x- and y-directions, respectively; Qx ¼ UH and Qy ¼ VH are fluxes per unit width; f is the Coriolis parameter; g is the gravita­ tional acceleration; Ps is the atmospheric pressure at the surface; ρ0 is the density of water; η is the Newtonian equilibrium tidal potential and α is the effective earth elasticity factor; τs , winds and τs , waves are surface stresses due to winds and waves, respectively; τb is the bottom stress; Mx & My are lateral stress gradients; Dx & Dy are momentum dispersion terms; and τ0 is a numerical parameter that optimizes the phase prop­ agation properties (Atkinson et al., 2004; Dietrich et al., 2012). Various studies have been reported on the wetting and drying algo­ rithms that ADCIRC uses to activate and deactivate elements based on the upturn and downturn of water depth (Blain et al., 2002; Dietrich

g ​ Sp þ

∂V ∂V ∂V þ Sp U þV ∂t ∂x ∂y

�

�

fV ¼

(5)

�

∂ Ps αη ζþ ∂y gρ0 τsy; wind þ τsy; waves þ ρ0 H g

τby

þ

My

Dy H

(6)

where H ¼ ζ þ h is the total water depth; ζ is the deviation of the water surface from the mean; h is the bathymetric depth with respect to the mean sea level; Sp ¼ cosy0 =cos y is the spherical coordinate conversion 5

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

Fig. 6. Maximum wind speed of the entire domain obtained from the model (GAHM) for both the cyclones. (a) Cyclone Thane; (b) Cyclone Vardah.

et al., 2004, 2006; Dick et al., 2013). An element is considered to be wet only when all three vertices are wet, if not then the corresponding element is taken as dry for that particular time-step. ADCIRC adopts implicit and explicit procedures to solve the GWCE. Implicit scheme operates on the principle of Jacobi Conjugate Gradient (JCG) method, whereas the explicit scheme uses lumped diagonal mass matrix which solves at a faster rate than implicit (Dietrich et al., 2004, 2012). Further, once the water levels are obtained with respect to wetting and drying algorithms, the momentum equations are solved to estimate the water currents.

Dingemans, 1997; Booij et al., 1999; Zijlema, 2010). Stot is the source or sink term representing the sum (Stot ¼ Sin þ Swc þ Snl4 þ Sbot þ Sbrk þ Snl3 ) of wave growth by wind (Sin Þ, action lost due to whitecapping (Swc Þ, surf breaking (Sbrk ), bottom friction (Sbot ), and the action exchanged between non-linear wave interaction components (Snl3 and Snl4 ). SWAN operates on the basis of action density spectrum N (σ, θ) instead of energy density spectrum E (σ, θ), since, in the presence of ambient currents the action density is conserved whereas energy density is not (Whitham, 1974). Hence the phase resolving is achieved steadily even at finer scales. Moreover, to maintain the stability on unstructured grids, SWAN adopts an analog label to the Gauss Siedel iteration method (Zijlema, 2010). Due to the above-mentioned key-characterizations, the circulation model (ADCIRC) and wave model (SWAN) are ‘coupled’ to run on a parallel environment with the same unstructured triangular mesh for the computation of storm wave and surge characteristics. In the fully coupled model, ADCIRC first solves the shallow water equations to evaluate the wind velocities, water levels and current velocities, which is then passed on to SWAN to force its computations. SWAN follows with the calculation of wave parameters and their corresponding radiation stresses, on the following time-step, which is then fed as an input to ADCIRC. The exchange of data between these two models on the parallel computing environment happens at each time step internally through local memory (Dietrich et al., 2012). The coupling time interval between these two models are based on SWAN’s time-step. For the current model simulations ADCIRC time-step was considered to be 1 s as it is limited to Courant number stabilization for its explicit method whereas, SWAN’s time-step is prescribed as 1800 s (owing to greater stability of the model).

3.2. Simulating Waves Nearshore (SWAN) model SWAN, the third generation wave model used in this study solves the action balance equation through the method of Gauss Siedel Algorithm to obtain the wave parameters. The advantage of SWAN over the other computational models is that, it can be applied over large domains to resolve the length and time scales of wind generated waves (Booij et al., 1999; Dietrich et al., 2011). The gradual development of the wind generated waves can be described through the governing wave action balance equation: � ⇀ ⇀� � ∂N ∂cθ N ∂cσ N Stot ¼ þ r⇀x : Cg þ U N þ þ ∂t ∂θ ∂σ σ

(7)

where, N (x, t, σ, θ) is the wave action density spectrum; change of wave action in time t; σ is the relative frequency; θ is the wave direction. The other terms of the equation represents, r⇀x the gradient operator in ⇀

⇀

geographic space, cg the wave group velocity and U the ambient current vector, turning rate ðcθ Þ, and shifting rate (cσ ). The relative terms for the expressions are from the linear wave theory (Whitham, 1974; ⇀

6

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

Fig. 7. Total water level elevation time series of Cyclone Thane obtained through all five methods compared with the observed tide gauge stations data. (a) Ennore; (b) Nagapattinam. The purple line indicates the landfall time of Thane (0200 UTC 30 Dec 2011).

3.3. Wind drag estimation methods

increase of wind speed (up to 41 m/s) and is capped by the numerical model (ADCIRC) at 0.0035 for high wind speeds (beyond 41 m/s) to restrict unrealistic water level outputs.

The most important driving factor of numerical surge and wave modelling is the wind stress that is calculated by the corresponding models. The wind stress which is directly proportional to its wind drag formulations can improve the accuracy of any surge or wave computa­ tion to a great extent when represented in a more realistic way. In the present work, five wind drag estimation methods namely, the WBLM (Wave Boundary Layer Model) and four most widely used empirical formulations developed by Garratt (1977), Large and Pond (1981), Wu (1982) and Zijlema et al. (2012) are considered for the investigation of wind stress computation and its corresponding storm surge estimation.

3.3.2. Large and Pond (L-P) method The relation developed by Large and Pond (1981) can be written as: 8 1:2 u10 < 11 m=s < 3 11 � u10 � 25 m=s (9) Cd � 10 ¼ 0:49 þ 0:065u10 : 2:115 u10 > 25 m=s This formulation is used in the Finite Volume Coastal Ocean Model (FVCOM) to simulate coupled surge and wave characteristics. It is cap­ ped at a constant value beyond 25 m/s wind speed. The major difference between equations (8) and (9) is the upper limit at which Cd is capped, which in turn implies that there is still no solidarity in estimating the wind drag coefficient at high wind speeds.

3.3.1. Garratt’s method The bulk-formulae proposed by Garratt (1977) is the most widely used in terms of coupled ADCIRC and SWAN modelling (Bunya et al., 2010; Dietrich et al., 2011, 2012) and it is given as: � u10 < 41 m=s 0:75 þ 0:067u10 Cd � 103 ¼ (8) 3:5 u10 > 41 m=s

3.3.3. Wu’s method Wu in 1982 formulated another empirical relation for estimating the wind drag and it is denoted as:

The Cd calculated through equation (8) increases linearly with the 7

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

Fig. 8. Total water level elevation time series of Cyclone Vardah obtained through all five methods compared with the observed tide gauge stations data. (a) Ennore; (b) Nagapattinam. The purple line indicates the landfall time of Vardah (1100 UTC 12 Dec 2016).

� Cd � 103 ¼

u10 < 7:5 m=s u10 � 7:5 m=s

1:2875 0:8 þ 0:065u10

3.3.5. Wave Boundary Layer Model (WBLM) The hindmost atmospheric layer over the ocean surface is termed to be the boundary layer, which is significantly affected due to the wave motion under extreme wind conditions. Over the past few decades’ re­ searchers have advocated to evaluate the wind stress that affects the wave motion through the correlation of momentum and energy con­ servation methods (Longuet-Higgins, 1957; Phillips, 1985; Janssen, 1989; Makin and Kudryavtsev, 1999; Hara and Belcher, 2002, 2004; Moon et al., 2004a, 2004b). The current study utilizes the enhanced Wave Boundary layer model developed by Chen and Yu (2016), where the wind shear stress ðτÞ is obtained by solving the conservation of mass and momentum equations for the atmospheric boundary layer:

(10)

The linear slope in the formulation corresponding to the wind speed is identical in both Wu and L-P’s methods, whereas the constants used with respect to boundary layer characteristic differs. The method is directly adopted in SWAVE model (Qi et al., 2009) which is the modified version of SWAN to couple with FVCOM and other circulation models. 3.3.4. Zijlema’s method Zijlema et al. (2012) had developed a quadratic method that fits a polynomial regression function using different sets of field data (re­ ported through various studies) with a decreasing trend in the drag coefficient at high wind speeds. The quadratic function is expressed as: 3

Cd � 10 ¼ 0:55 þ 2:97~u

1:49~u

2

d dΠw ðu ⋅ τtot Þ þ dz dz

(11)

τtot ¼ τw ðzÞ þ τt ðzÞ

~ ¼ u10 =31:5, and u10 is measured in terms of m/s. Equation (11) where u gives a quadratic approach to estimate the wind drag coefficient at high wind speeds (in particular during cyclonic conditions > 31.5 m/s), which is 30% less than that of Wu’s method (as shown in Fig. 2). However, Zijlema’s method fails to give a positive value at wind speeds >68 m/s, and thus a minimal constraint has to be added.

ρa ε þ gðρw

ρa Þ α s k

ds ¼0 dz

(12) (13)

where, τw and τt are the respective wave induced stress and turbulent stress occur locally; u is the wind vector at the height z above the sea surface; Πw is the wave-induced motion’s vertical flux energy; ε is the viscous dissipation rate; αs is the reciprocal of turbulent Schmidt num­ ber; k is the turbulent eddy viscosity; and s is the sea spray’s volumetric 8

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

Fig. 9. Maximum water level elevation of the entire domain obtained from the model through WBLM method for both the cyclones. (a) Cyclone Thane; (b) Cyclone Vardah. The dotted line in (a) and (b) represents the cyclone’s track respectively.

concentration at height z above the sea surface (Chen and Yu, 2017). The wave-induced stress is computed through the spectrum obtained from wave model through the bulk formulae. It is to be noted that there are different formulations used by models to determine the high-frequency spectrum fed into WBLM (Donelan et al., 2012; Reichl et al., 2014). Fig. 2 represents the comparison of wind drag estimations

at various wind speeds for WBLM and three other bulk formulae’s used in this study. The wind drag coefficients calculated through WBLM shows a lowering trend instead of levelling-off at very high wind speeds similar to the field measurements reported by Powell et al. (2003) (as shown in Fig. 2). The limitation to the current modelling study is that it requires a pre-run output of the wave spectrum from a numerical wave 9

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

the viscous Stress (τv ), wave-Induced Stress (τw ), roughness length and wind drag coefficient. Eventually, the wind shear stress ðτÞ is computed and given back as an input to the coupled circulation and wave model to finally obtain the water levels and wave attributes for both the cyclonic conditions. The model domain with bathymetric depth is shown in Fig. 4, which also illustrates the spatial extent that covers a significant part of the BoB basin. The domain was discretised with an unstructured triangular mesh of 463,752 elements and 233,802 nodes, the bathymetric data from GEBCO (General Bathymetric Chart of the Oceans) with 30 arc second data interval along with a high resolution NHO (National Hydrographic Office, Govt. of India) chart was used to define the water depth at all the nodal points. Mesh size of approximately 200–1000 m was discretised along the near shore regions, while the deep ocean was discretised with 15–20 km resolution roughly. As mentioned earlier, the simulations were performed for two cyclonic scenarios namely Thane and Vardah, which occurred in 2011 and 2016, respectively. The simulation time frame for Thane was from 25th December 2011 12:00 h to 30th December 12:00 h UTC (5 Days), whereas for Vardah it was executed from 6th December 2016 6:00 h to 13th December 2016 6:00 h UTC (7 days). The wind and pressure inputs given to the model is through the best track information obtained from the Joint Typhoon Warning Centre (JTWC) using which the wind and pressure fields are constructed through the Generalised Asymmetric Holland Model (GAHM, an inbuilt module in ADCIRC) for the simulations. GAHM interpolates wind speeds to all the nodal points of the domain through the information specified in the best track data by JTWC. The bottom friction was assigned to be hybrid with parameters like Cb, Hb, θ, and λ are set to 0.0026, 1, 10 and 1/3 respectively (Bhaskaran et al., 2014). Time-step for the circulation model on both the cyclonic cases was considered to be 1 s, which sat­ isfies the stability (CFL) criteria (Westerink et al., 1994) with respect to the domain resolution. On the other hand, SWAN’s time-step was posed to 1800 s and the coupling time interval was matched between ADCIRC

Table 1 Statistical comparison of mean and standard deviation for the estimates of water level elevations computed using different wind drag methods and the observed data. Cyclone Thane – WLE (meters) Ennore Wu L-P Garratt Zijlema WBLM Observed

Nagapattinam

Mean

Standard Deviation

Mean

0.1494 0.0724 0.0926 0.0564 0.0494 0.0083

0.3497 0.2847 0.3033 0.2789 0.2602 0.3308

0.0615 0.0021 0.0255 0.0205 0.0187 0.0008

Standard Deviation 0.2200 0.2105 0.2215 0.2190 0.2230 0.2331

Cyclone Vardah – WLE (meters) Ennore Wu L-P Garratt Zijlema WBLM Observed

Nagapattinam

Mean

Standard Deviation

Mean

Standard Deviation

0.0840 0.0519 0.0675 0.0486 0.0408 0.0007

0.2773 0.2242 0.2644 0.2409 0.2383 0.2388

0.0232 0.0137 0.0177 0.0107 0.0096 0.0017

0.1950 0.1826 0.1924 0.1755 0.1777 0.2092

model (which is computationally expensive and time consuming). However, it is negligible as the current day wave models can be coupled with WBLM to execute the same. 3.4. Methodology The simulations of the present study (as represented with the methodology flow diagram shown in Fig. 3) were carried out with the coupled ADCIRC (circulation model) and SWAN (wave model) on an unstructured triangular mesh, whose two dimensional directional wave spectrum output is taken as an input forcing into the WBLM to calculate Cyclone Thane - WLE (Ennore)

1.2

Statistical Estimates

1

0.86

0.94

0.89

Wu

L-P

Garatt

Zijlema

1

0.6

0.2

0.29

0.26 0.11

0.16

0.11

0.06 0.07

0.07 0.09 0.05

0

-0.2

RMSE

Cyclone Thane - WLE (Nagapattinam) 1.2

0.97 0.98

0.8

0.4

b.

WBLM

R

0.14 0.17 0.11

0.03

0.09

Statistical Estimates

a.

0.89

Statistical Estimates

1 0.81

0.8

0.86

0.83

L-P

PE (x100)

Garatt

Zijlema

0.2

0.31

0.4 0.2

0.20

0.24 0.10 0.13 0.08 0.08

0.09

0.36 0.24 0.25 0.21

0.19

RMSE

R

0.25 0.27 0.24 0.22

MBE Statistical Methods

Cyclone Vardah - WLE (Nagapattinam) 1.2

0 -0.2

0.11

0.07 0.05 0.04

0.16 0.07 0.07

-0.03

1

0.88 0.88

0.29

0.24 0.27 0.22 0.21

WBLM

0.6

d.

WBLM

0.6 0.4

Zijlema

0.8

-0.2

Statistical Estimates

1.2

Wu

Garatt

0

MBE

Cyclone Vardah - WLE (Ennore)

L-P

0.93 0.90 0.96 0.97

Statistical Methods

c.

Wu

0.88

Wu

L-P

PE (x100)

Garatt

Zijlema

WBLM

0.92 0.90 0.94 0.95

0.8 0.6 0.4 0.2

0.28

0.23 0.13

0.18

0.18

0.11 0.11

0.11

0.15

0.21 0.24 0.20 0.19

0.08 0.09

0 RMSE

R

MBE

-0.2

PE (x100)

Statistical Methods

RMSE

R

MBE

PE (x100)

Statistical Methods

Fig. 10. Statistical error analysis comparison between the observed and different wind drag methods for the computation of water level elevations of cyclones’ Thane (a and b) and Vardah (c and d), respectively. 10

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

Fig. 11. Significant wave height time series of Cyclone Thane obtained through all five methods compared with the observed buoy data. (a) BD11; (b) Wave Rider Buoy Pondicherry.

and SWAN (i.e. the exchange of wind velocities, water levels and current velocities from ADCIRC and wave radiation stresses from SWAN will happen every 1800 s). The computational space was discretised to 72 wave directional intervals with 5 � increment and the wave frequency intervals are sorted out to 40 with a logarithmic increase across the range 0.0285–1.1726. Energy changes due to white-capping, depth-in­ duced wave breaking and wave-wave interaction were also activated. All five wind stress evaluation methods (Garratt, Wu, L-P, Zijlema and WBLM) were implemented separately to calculate the surge and wave attributes.

statistical performance indicators such as the Mean, Standard Deviation, Root Mean Square Error (RMSE), Correlation Coefficient (R), Mean Bias Error (MBE) and Mean Relative Error Accuracy or Percentage Error (PE) were used. The primary assumption made to compute these statistical parameters were that the observed (measured) data are error free. While, the generalised mean and standard deviation methods were calculated, the RMSE was computed using the following relationship, rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1X ðM OÞ2 RMSE ¼ (14) n Where ‘M’ represents the model predicted and ‘O’ represents the observed values for a set of n. RMSE derives the line of best fit between the model predicted and the observed values (Umesh et al., 2017). It is inferred that the lower value of RMSE gives the best fit between model

3.5. Statistical error analysis To evaluate and compare the performance of the model generated results using various wind drag methods and the observed values, 11

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

Fig. 12. Significant wave height time series of Cyclone Vardah obtained through all five methods compared with the observed buoy data. (a) BD11; (b) Wave Rider Buoy Pondicherry.

and observed data. Correlation coefficient is another often used method to establish the relation between the model and observed values. It signifies the strength of the relationship between the model and observed parameters with a numerical value (higher value if strong correlation persists and lower value to zero if they are independent). It is defined as follows, P ðM MÞ ðO OÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ¼ qffiP (15) ðM MÞ2 ðO OÞ2

where M and O represents the mean of the model predicted and observed values respectively. If the value of the correlation coefficient is close to one then it indicates that there is a high positive correlation between model and observed data with low random error, however, the systematic error is not induced into it. Mean Bias Error (MBE) is the average magnitude of the set of errors between prediction and observed data that does not include the direc­ tion (Padilla-Hernandez et al., 2007). If the signs of the errors are included, then it is called as the Mean Absolute Error (MAE). The MBE is represented as, 12

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

Fig. 13. Maximum significant wave height of the entire domain obtained from the model through WBLM method for both the cyclones. (a) Cyclone Thane; (b) Cyclone Vardah. The dotted line in (a) and (b) represents the cyclone’s track respectively.

MBE ¼

1X ðM n

OÞ

modelling viewpoint, the Percentage Error (PE) was computed (Brown, 2010). Percentage Error that is also known as the mean relative error accuracy is defined as the difference between the modelled and measured parameters with respect to the percentage of the measured field data.

(16)

The mean bias error signifies the trends in the error, i.e. whether the model consistently over-predicts or under-predicts with respect to the field measurement (Remya et al., 2012; Umesh et al., 2017). In order to better evaluate the performance of the model in terms of numerical 13

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

4. Results and discussion

Table 2 Statistical comparison of mean and standard deviation for the estimates of sig­ nificant wave heights (Hs) computed using different wind drag methods and the observed data.

4.1. Estimation of water levels The coupled model (ADCIRC þ SWAN) simulations were performed with all five wind stress evaluation methods (Wu’s method, Garratt’s method, L-P’s method, Zijlema’s method and WBLM) to compute the water levels, current velocities, significant wave heights at all the nodal points in the domain. The wind field for the entire domain was con­ structed using the GAHM model based on the JTWC cyclone track in­ formation and it was compared with the field data recorded from the deep-water buoy BD11. Fig. 5 (a and b) represents the comparison of the model computed wind speeds with the observed data for both the cyclones (Thane and Vardah), which clearly indicates that the model has captured the wind speeds and profiles very well along with the maximum wind intensities (as the cyclones progressed towards the coast). The maximum wind speeds (m/s) obtained at all the nodal lo­ cations in the domain is represented in Fig. 6, which shows the extent of high wind speeds caused by cyclone Thane and Vardah. At the tide gauge stations, Ennore and Nagapattinam, wind intensities of 19 m/s were noticed for cyclone Thane, while for cyclone Vardah it was 37 m/s in Ennore and 14 m/s in Nagapattinam. The peak wind speeds at the deep-water buoy location (BD11) was found to be 14 m/s and 22 m/s respectively for cyclone Thane and Vardah, and in the shallow-water buoy point (WRBP) it was 31.5 m/s and 15 m/s for cyclone Thane and Vardah respectively. The model generated wind field was then forced in to the coupled model to estimate the water levels and wave parameters. The water levels computed at the coastal locations (Ennore and Nagapattinam) using all five wind stress evaluation methods for both the cyclones were compared with the available observed values obtained from the Ennore and Nagapattinam tide gauge stations. The observed data from both the tide gauge stations were at 1 min frequency with severe noised values, henceforth, the moving average technique was adopted to smoothen out the measured data and the smoothened data

Cyclone Thane - Hs (meters) BD11 Wu L-P Garratt Zijlema WBLM Observed

WRBP

Mean

Standard Deviation

Mean

Standard Deviation

2.8620 1.5995 2.3934 2.0648 2.0181 3.6793

1.3027 2.1674 1.3669 2.1086 2.0307 0.9137

2.3764 1.4151 1.9330 1.7375 1.6922 1.8793

1.2257 0.8541 0.9866 0.6987 0.7354 0.2808

Cyclone Vardah - Hs (meters) BD11 Wu L-P Garratt Zijlema WBLM Observed

PE ¼ 100

Mean

Standard Deviation

Mean

Standard Deviation

3.2649 2.4797 4.1473 3.3057 3.3619 3.3700

2.4586 1.9333 2.7402 1.6925 1.6244 1.5195

0.9166 0.8268 1.3015 1.2268 1.2387 1.2718

0.9736 0.7127 0.9826 0.6199 0.5828 0.3561

n X Mi i¼1

WRBP

Oi Oi

(17)

All the above mentioned statistical methods of error analysis were used in this study to evaluate the performance of wind drag estimation methods for numerically computing the storm wave and surge parameters.

Fig. 14. Statistical error analysis comparison between the observed and different wind drag methods for the computations of significant wave heights of cyclones’ Thane (a and b) and Vardah (c and d), respectively. 14

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

was used for the comparison of model results (as represented in Fig. 7 and Fig. 8). It is observed that the water levels computed based on the WBLM and Zijlema’s method (as both these methods accounts the decreasing trend of drag coefficient at high wind speed) for cyclone Thane shows a very good comparison with the measurements (as shown in Fig. 7 a and b). Similarly, the WBLM and Zijlema’s method based water levels estimated for Vardah cyclone also shows a closer match with the available observations as represented in Fig. 8 (a and b). As the tide gauge station at Ennore was located along the Vardah cyclone track, the sensors stopped recording after 0600 h UTC on 12th December 2016 due to extreme weather conditions. This restricted the comparison of maximum water levels during the cyclone’s landfall. Fig. 8 (a) distinctly explains the influence of the cyclonic surge component in the change of total water level elevations with respect to its mean at the time of cyclone Vardah’s landfall (1100 UTC 12th Dec 2011). Considering the fact that the winds are stronger on the right side of a cyclone’s track due to its anti-clock wise forward movement on the Northern Hemisphere, higher surge levels are observed at Ennore (located on the right side of the track, Fig. 8a) than Nagapattinam (Fig. 8b). The maximum water levels were evaluated at all the nodal points in the domain using the WBLM method for both the cyclones and are displayed in Fig. 9 (a and b). Fig. 9 clearly explains the change in total water levels (in the entire domain) with respect to its natural tide and surge shoved variation. It was also observed that the maximum water levels are higher along the cyclone track and the near shore regions. The maximum water level elevation for cyclones Thane and Vardah at the time of their corre­ sponding landfall (landfall point and its nearby locations) was estimated to be around 1–1.5 m by the WBLM-based model, which also corre­ sponds well with the IMD reported maximum water level of 1.2 m. Statistical error analysis was carried out for the estimated water levels by computing the Mean, Standard Deviation (as represented in Table 1), Root Mean Square Error (RMSE), Correlation Coefficient (R), Mean Bias Error (MBE) and Percentage Error (PE) for both the cyclones and are shown in Fig. 10. Fig. 10 (a and b) and Fig. 10 (c and d) shows the RMSE, R, MBE and PE for cyclones Thane and Vardah, respectively. It can be noticed through the error analysis results that in most cases the Zijlema and WBLM-based water level computations consistently per­ forms better than other linear wind drag methods (Wu’s method, Gar­ ratt’s method, L-P’s method) with lower RMSE, PE and higher R. Although, Zijlema’s method correlates well with the observed water level data (as the wind speed of Thane and Vardah did not exceed 42 m/ s), the method may not hold good for calculating the drag coefficient when the wind speed is greater than 55 m/s (as the quadratic function gives a higher decreasing slope for wind speeds greater than 55 m/s, which may lead to underestimation of surge). The comparison and analysis of peak water levels estimated by the models at the time of landfall could not be performed due to the inadequacy of observed data. In the case of Cyclone Thane neither of the two tidal stations (Ennore and Nagapattinam) were in near vicinity region of the cyclone’s landfall. Whereas in case of Cyclone Vardah, Ennore station was close by but unable to record data at the landfall time due to extreme weather con­ ditions. However, the present investigation of the storm surge compu­ tations using five different wind stress methods compared proves the fact that, WBLM is the best suitable method to estimate the wind stress factor, as it can even be applicable for cyclones having wind speeds above 55 m/s.

using all the five wind stress methods were compared with the buoy data fetched from stations situated closer to the cyclone’s track and its associated landfall regions. The model computed significant wave heights (Hs) were compared with the in-situ buoy data at both deepwater buoy (BD11) and shallow-water buoy (Wave Rider Buoy Pondi­ cherry (WRBP)) locations for cyclone Thane and Vardah, and are shown in Fig. 11 (a and b) and Fig. 12 (a and b), respectively. The analysis of Figs. 11 and 12 shows that, Wu and Garratt’s methods over-estimates the wave heights due to the implementation of linear formulation of Cd (with higher slope and higher cut-off values) whereas, L-P’s method significantly underestimates Hs due to the implementation of Cd formulae with reduced slope and lower cut-off values. Zijlema’s and WBLM-based Hs results demonstrated a rather succinct and balanced performance, despite a certain (minor) degree of underestimation of the significant wave heights, thereby rendering higher preference to its re­ sults over that of the other three bulk-formulae’s. As the cyclonic wind field is given as an input without accounting for the pre- and postexisting regular wind field, there is an apparent artefact of underesti­ mation observed in both Zijlema’s and WBLM-based model estimates. The Hs computed using Zijlema and WBLM methods at WRBP due to cyclone Thane (Fig. 11 (b)) and BD11 due to cyclone Vardah (Fig. 12 (a)) are found to agree well with the in-situ data at respective locations. This agreement was achieved, as the buoy locations were present very close to the cyclones’ track, where the influence of regular (pre and post existing) wind conditions were minimal. Thus, affirmatively, the wind field generated from its best track was good enough to capture the peak storm wave height by the model (through Zijlema’s and WBLM method). Although, SWAN was relatively known to be a nearshore wave model (as it quite comfortably solves the shallow water dynamics), it also accomplished effective results in the deep-water locations for extreme weather conditions. The activation terms like white-capping, depth induced breaking and triad wave-wave interaction processes that was incited in the wave model were well executed in the deep-water regions too. As mentioned earlier, for better estimation of storm surge time se­ ries, the pre-existing wind and wave conditions in the entire area needs to be included in the numerical modelling. Fig. 13 delineates the spatial variation of the maximum Hs obtained in the entire domain using WBLM method, for the complete duration of both the cyclones. The maximum Hs shown in Fig. 13 also enumerates that the winds were quite stronger on the right side of these cyclones, which produced larger wave heights on that side compared to left side of the track. Cyclones Thane and Vardah generated maximum significant wave heights of 11.4 m and 11.2 m on their respective forward move­ ments from the deep water towards the shorelines. The deep-water storm waves created by these cyclones reduce their energy towards to the near-shore areas due to shoaling and wave breaking, could also be visualized in Fig. 13. The model computed storm wave parameters were also evaluated through various statistical error analysis methods with respect to the observed data. The mean and standard deviation, calcu­ lated and tabulated in Table 2, indicates the consistency of the signifi­ cant wave heights computed by the coupled model using the WBLM method with respect to the observed data. Fig. 14 (a and b) and Fig. 14 (c and d) represents the RMSE, R, MBE and PE for cyclones’ Thane and Vardah, respectively. The significant wave heights estimated through the WBLM method also attained maximum R and minimum RMSE, PE (as shown in Fig. 14) in comparison to the other three linear wind drag estimation methods. The over and under estimation of the model results with respect to the observed data was also clearly indicated through the MBE analysis of Fig. 14. In Fig. 14 (a), the WBLM and Zijlema’s method based results of cyclone Thane is under predicted, as the errors are attributed due to considering only the cyclone best track information and not the pre-existing wind wave climate (i.e. the corresponding buoy location BD11 is almost 250 km away from cyclone Thane’s track and the wind model considers only the best track wind information). The error analysis results portrayed through Fig. 14, demonstrates the effi­ cacy of WBLM method based computation of storm wave heights. The

4.2. Estimation of significant wave heights It is necessary to estimate the significant wave heights in the domain and validate the model capability, as it greatly influences the accuracy of any storm surge modelling study. The model (SWAN) computed signif­ icant wave heights were compared with in-situ data to evaluate the performance of the wind drag estimation methods. The in-situ data (buoy data) was obtained from the Indian National Centre for Ocean Information Services (INCOIS). The simulated significant wave heights 15

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509

WBLM and Zijlema’s quadratic function was found to generate the realistic representation of wind stress and storm wave attributes for wind speeds greater than 30 m/s, to better simulate the cyclonic storm surge. However, as aforementioned in section 4.1, Zijlema’s method can be used for cyclones that has wind speeds up to 55 m/s, as the numerical models may tend to underestimate the surge and wave components if the quadratic function is adopted to estimate the wind drag coefficient for cyclones that exceeds 55 m/s wind speed (category 4 and 5 hurricanes, that are very common in the Atlantic Ocean). Thus, it is recommended to compute the wind stress using the WBLM, as it can be used for cyclones/ hurricanes with extreme wind conditions (up to 80 m/s) to better predict the storm surge and wave attributes.

for Ocean Information Services (INCOIS) for providing the field obser­ vation data from the tide gauge and buoy stations. Authors also thank the Indian Meteorological Department (IMD) for the free access to available cyclone atlas and reports. References Andreas, E.L., 2004. Spray stress revisited. J. Phys. Oceanogr. 34 (6), 1429–1440. Atkinson, J.H., Westerink, J.J., Hervouet, J.M., 2004. Similarities between the quasibubble and the generalized wave continuity equation solutions to the shallow water equations. Int. J. Numer. Methods Fluids 45 (7), 689–714. Bhaskaran, P.K., Gayathri, R., Murty, P.L.N., Bonthu, S., Sen, D., 2014. A numerical study of coastal inundation and its validation for Thane cyclone in the Bay of Bengal. Coast. Eng. 83, 108–118. Blain, C.A., Preller, R.H., Rivera, A.P., 2002. Tidal prediction using the advanced circulation model (ADCIRC) and a relocatable PC-based system. Oceanography 15 (1), 77–87. Booij, N.R.R.C., Ris, R.C., Holthuijsen, L.H., 1999. A third-generation wave model for coastal regions: 1. Model description and validation. J. Geophys. Res.: Oceans 104 (C4), 7649–7666. Brown, J.M., 2010. A case study of combined wave and water levels under storm conditions using WAM and SWAN in a shallow water application. Ocean Model. 35 (3), 215–229. Bunya, S., Dietrich, J.C., Westerink, J.J., Ebersole, B.A., Smith, J.M., Atkinson, J.H., Jensen, R., Resio, D.T., Luettich, R.A., Dawson, C., Cardone, V.J., 2010. A highresolution coupled riverine flow, tide, wind, wind wave, and storm surge model for southern Louisiana and Mississippi. Part I: model development and validation. Mon. Weather Rev. 138 (2), 345–377. Chen, Y., Yu, X., 2016. Enhancement of wind stress evaluation method under storm conditions. Clim. Dyn. 47 (12), 3833–3843. Chen, Y., Yu, X., 2017. Sensitivity of storm wave modeling to wind stress evaluation methods. J. Adv. Model. Earth Syst. 9 (2), 893–907. Dawson, C., Westerink, J.J., Feyen, J.C., Pothina, D., 2006. Continuous, discontinuous and coupled discontinuous–continuous Galerkin finite element methods for the shallow water equations. Int. J. Numer. Methods Fluids 52 (1), 63–88. Dick, C.G., Tromble, E.M., Dresback, K.M., Kolar, R.L., 2013. Implementation and analysis of a partial-element wetting and drying framework for generalized wave continuity equation-based hydrodynamic models. Int. J. Numer. Methods Fluids 72 (10), 1015–1033. Dietrich, J.C., Bunya, S., Westerink, J.J., Ebersole, B.A., Smith, J.M., Atkinson, J.H., Jensen, R., Resio, D.T., Luettich, R.A., Dawson, C., Cardone, V.J., 2010. A highresolution coupled riverine flow, tide, wind, wind wave, and storm surge model for southern Louisiana and Mississippi. Part II: synoptic description and analysis of Hurricanes Katrina and Rita. Mon. Weather Rev. 138 (2), 378–404. Dietrich, J.C., Kolar, R.L., Luettich, R.A., 2004. Assessment of ADCIRC’s wetting and drying algorithm. In: Miller, C.T., Farthing, M.W., Gray, W.G., Pinder, G.F. (Eds.), Proceedings of Computational Methods in Water Resources, vol. 2, pp. 1767–1778. Dietrich, J.C., Kolar, R.L., Westerink, J.J., 2006. Refinements in continuous Galerkin wetting and drying algorithms. In: Estuarine and Coastal Modeling, pp. 637–656. Dietrich, J.C., Tanaka, S., Westerink, J.J., Dawson, C.N., Luettich, R.A., Zijlema, M., Holthuijsen, L.H., Smith, J.M., Westerink, L.G., Westerink, H.J., 2012. Performance of the unstructured-mesh, SWANþ ADCIRC model in computing hurricane waves and surge. J. Sci. Comput. 52 (2), 468–497. Dietrich, J.C., Zijlema, M., Westerink, J.J., Holthuijsen, L.H., Dawson, C., Luettich Jr., R. A., Jensen, R.E., Smith, J.M., Stelling, G.S., Stone, G.W., 2011. Modeling hurricane waves and storm surge using integrally-coupled, scalable computations. Coast. Eng. 58 (1), 45–65. Dingemans, M.W., 1997. Water Wave Propagation over Uneven Bottoms, vol. 13. World Scientific. Donelan, M.A., Curcic, M., Chen, 1.S., Magnusson, A.K., 2012. Modeling waves and wind stress. J. Geophys. Res.: Oceans 117 (C11). Donelan, M.A., Haus, B.K., Reul, N., Plant, W.J., Stiassnie, M., Graber, H.C., Brown, O.B., Saltzman, E.S., 2004. On the limiting aerodynamic roughness of the ocean in very strong winds. Geophys. Res. Lett. https://doi.org/10.1029/2004GL019460. Dube, S.K., Chittibabu, P., Rao, A.D., Sinha, P.C., Murty, T.S., 2000. Extreme sea levels associated with severe tropical cyclones hitting Orissa coast of India. Mar. Geod. 23 (2), 75–90. Dube, S.K., Jain, I., Rao, A.D., Murty, T.S., 2009. Storm surge modelling for the Bay of Bengal and Arabian Sea. Nat. Hazards 51 (1), 3–27. Dube, S.K., Rao, A.D., Sinha, P.C., Murty, T.S., Bahulayan, N., 1997. Storm surge in the Bay of Bengal and Arabian Sea the problem and its prediction. Mausam 48 (2), 283–304. Garratt, J.R., 1977. Review of drag coefficients over oceans and continents. Mon. Weather Rev. 105 (7), 915–929. Hara, T., Belcher, S.E., 2002. Wind forcing in the equilibrium range of wind-wave spectra. J. Fluid Mech. 470, 223–245. Hara, T., Belcher, S.E., 2004. Wind profile and drag coefficient over mature ocean surface wave spectra. J. Phys. Oceanogr. 34 (11), 2345–2358. Holthuijsen, L.H., Powell, M.D., Pietrzak, J.D., 2012. Wind and waves in extreme hurricanes. J. Geophys. Res.: Oceans 117 (C9). Huang, Y., Weisberg, R.H., Zheng, L., 2010. Coupling of surge and waves for an Ivan-like hurricane impacting the Tampa Bay, Florida region. J. Geophys. Res.: Oceans 115 (C12).

5. Summary and conclusions Numerical modelling of cyclonic storm surge with higher precision is important to enhance the prediction and mitigation planning systems to maintain a better coastal resilience. The accuracy of numerical storm wave and storm surge estimation primarily depends on the wind forcing (wind stress or wind drag) represented by the wind drag coefficient. The present study aimed at modelling of near-accurate storm waves and surge using five different wind drag estimation methods to ascertain the most suitable method for tropical cyclones (i.e. Thane and Vardah) in the Bay of Bengal basin. The vital conclusions drawn from the research work are, 1. The water levels estimated using Zijlema’s and WBLM method closely corresponded with the observed values and found to be better compared to the other three linear empirical formulations (Garratt, L-P and Wu). The WBLM-based water levels rendered the maximum correlation coefficient (R ~ 0.88–0.98) with respect to the observed field data which is higher compared to that of the other linear and quadratic methods (Garratt (R ~ 0.83–0.90); L-P (R ~ 0.86–0.94); Wu (R ~ 0.81–0.89); Zijlema (R ~ 0.88–0.97)). The spatial distri­ bution of maximum water levels propelled by the cyclones to the shoreline, obtained based on the WBLM method shows a credible agreement with the IMD reported values. 2. The comparison of model computed significant wave heights with the two in-situ buoy (BD11, deep-water buoy and WRBP, shallowwater buoy) data, also clearly shows the efficiency of the WBLM method over the other wind stress methods in extreme wind condi­ tions. The simulated storm wave heights using WBLM at both deepwater and shallow-water regions were observed to have high corre­ lation coefficients ranging from 0.80 to 0.97 for capturing the peak significant wave height. However, the building-up and receding wave heights could be better simulated by accounting for the preand post-existing wind and wave conditions in the entire area. 3. Storm surge simulations should be carried out using a coupled cir­ culation and wave model using a robust wind stress evaluation method to include the realistic scenario of air-sea interaction process in extreme cyclonic events. The present study recommends the WBLM method for the evaluation of wind stress, as it is potent enough and includes the pragmatic energetic exchange at the air-sea interface in computation of wind drag coefficient for high wind conditions (exceeding 35 m/s). Although Zijlema’s method has simulated better results in the present study (as the maximum wind speeds of Thane and Vardah were limited to 42 m/s), the quadratic function may not be applicable for wind speeds greater than 55 m/s. Conclusively, it could be stated that WBLM method should be adopted in the numerical computations of near real-time storm surge prediction studies working towards asserting modelled estimates with enhanced precision. Acknowledgements The authors would like to acknowledge the Indian National Centre 16

C.G. Shankar and M.R. Behera

Ocean Engineering 191 (2019) 106509 Murty, T.S., Flather, R.A., Henry, R.F., 1986. The storm surge problem in the Bay of Bengal. Prog. Oceanogr. 16 (4), 195–233. Padilla-Hernandez, R., Perrie, W., Toulany, B., Smith, P.C., 2007. Modeling of two northwest Atlantic storms with third-generation wave models. Weather Forecast. 22 (6), 1229–1242. Phan, L.T., Slinn, D.N., Kline, S.W., 2013. Wave effects on hurricane storm surge simulation. In: Advances in Hurricane Engineering: Learning from Our Past, pp. 753–764. Phillips, O.M., 1985. Spectral and statistical properties of the equilibrium range in windgenerated gravity waves. J. Fluid Mech. 156, 505–531. Powell, M.D., Vickery, P.J., Reinhold, T.A., 2003. Reduced drag coefficient for high wind speeds in tropical cyclones. Nature 422 (6929), 279. Qi, J., Chen, C., Beardsley, R.C., Perrie, W., Cowles, G.W., Lai, Z., 2009. An unstructuredgrid finite-volume surface wave model (FVCOM-SWAVE): implementation, validations and applications. Ocean Model. 28 (1–3), 153–166. Remya, P.G., Kumar, R., Basu, S., Sarkar, A., 2012. Wave hindcast experiments in the Indian Ocean using MIKE 21 SW model. J. Earth Syst. Sci. 121 (2), 385–392. Reichl, B.G., Hara, T., Ginis, I., 2014. Sea state dependence of the wind stress over the ocean under hurricane winds. J. Geophys. Res.: Oceans 119 (1), 30–51. Tanaka, S., Bunya, S., Westerink, J.J., Dawson, C., Luettich, R.A., 2011. Scalability of an unstructured grid continuous Galerkin based hurricane storm surge model. J. Sci. Comput. 46 (3), 329–358. Umesh, P.A., Bhaskaran, P.K., Sandhya, K.G., Nair, T.B., 2017. An assessment on the impact of wind forcing on simulation and validation of wave spectra at coastal Puducherry, east coast of India. Ocean Eng. 139, 14–32. Westerink, J.J., Luettich, R.A., 1991. Tide and Storm Surge Prediction in the Gulf of Mexico Using Model ADCIRC-2D, Report to the US Army Engineers Waterways Experiment Station, 112, Vicksburg, Miss., July. Westerink, J.J., Blain, C.C., Luettich, R.A., Scheffner, N.W., 1994. ADCIRC — an advanced three dimensional circulation model for shelves, coasts and estuaries. Report 2: user’s Manual for ADCIRC-2DDI. Tech. Report DRP-92-6. U.S. Army Corps of Engineers, Washington D.C. Westerink, J.J., Luettich, R.A., Feyen, J.C., Atkinson, J.H., Dawson, C., Roberts, H.J., Powell, M.D., Dunion, J.P., Kubatko, E.J., Pourtaheri, H., 2008. A basin-to channelscale unstructured grid hurricane storm surge model applied to southern Louisiana. Mon. Weather Rev. 136 (3), 833–864. Whitham, G.B., 1974. Linear and Nonlinear Waves. John Wiley, New York, p. 636. Wu, J., 1982. Wind-stress coefficients over sea surface from breeze to hurricane. J. Geophys. Res.: Oceans 87 (C12), 9704–9706. Zheng, L., Weisberg, R.H., Huang, Y., Luettich, R.A., Westerink, J.J., Kerr, P.C., Donahue, A.S., Crane, G., Akli, L., 2013. Implications from the comparisons between two-and three-dimensional model simulations of the Hurricane Ike storm surge. J. Geophys. Res.: Oceans 118 (7), 3350–3369. Zijlema, M., 2010. Computation of wind-wave spectra in coastal waters with SWAN on unstructured grids. Coast. Eng. 57 (3), 267–277. Zijlema, M., Van Vledder, G.P., Holthuijsen, L.H., 2012. Bottom friction and wind drag for wave models. Coast. Eng. 65, 19–26.

Indian Meteorological Department IN, 2011. National Bulletin for Cyclone Warning. http://www.rsmcnewdelhi.imd.gov.in/images/pdf/archive/bulletins/2011/bulleti n_for_indian_coast1.pdf/. (Accessed 8 March 2018). Indian Meteorological Department IN, 2012. IMD/WMO Report on Cyclonic Disturbances over North Indian Ocean during 2011. http://rsmcnewdelhi.imd.gov. in/images/pdf/publications/annual-rsmc-report/rsmc-2011.pdf. (Accessed 8 March 2018). Indian Meteorological Department IN, 2016. National Bulletin for Cyclone Warning. http://www.rsmcnewdelhi.imd.gov.in/images/pdf/archive/bulletins/2016/NROA. pdf/. (Accessed 30 March 2018). Indian Meteorological Department IN, 2017. IMD/WMO Report on Cyclonic Disturbances over North Indian Ocean during 2016. http://rsmcnewdelhi.imd.gov. in/images/pdf/publications/annual-rsmc-report/rsmc-2016.pdf. (Accessed 30 March 2018). Janssen, P.A., 1989. Wave-induced stress and the drag of air flow over sea waves. J. Phys. Oceanogr. 19 (6), 745–754. Jelesnianski, C.P., Chen, J., Shaffer, W.A., 1992. SLOSH: sea, lake, and overland surges from hurricanes. In: NOAA Tech. Rep. NWS 48, 77, Natl. Weather Serv., National Oceanic and Atmospheric Administration, Silver Spring, Md. Kolar, R.L., Westerink, J.J., Cantekin, M.E., Blain, C.A., 1994. Aspects of nonlinear simulations using shallow-water models based on the wave continuity equation. Comput. Fluid 23 (3), 523–538. Kudryavtsev, V.N., 2006. On the effect of sea drops on the atmospheric boundary layer. J. Geophys. Res.: Oceans 111 (C7). Kudryavtsev, V.N., Makin, V.K., 2011. Impact of ocean spray on the dynamics of the marine atmospheric boundary layer. Boundary-Layer Meteorol. 140 (3), 383–410. Large, W.G., Pond, S., 1981. Open ocean momentum flux measurements in moderate to strong winds. J. Phys. Oceanogr. 11 (3), 324–336. Longuet-Higgins, M.S., 1957. The statistical analysis of a random, moving surface. Phil. Trans. R. Soc. Lond. A 249 (966), 321–387. Luettich, R.A., Westerink, J.J., 2004. Formulation and numerical implementation of the 2D/3D ADCIRC finite element model version 44. XX (74). http://adcirc.org/adcirc_ theory_2004_12_08.pdf. Madsen, H., Jakobsen, F., 2004. Cyclone induced storm surge and flood forecasting in the northern Bay of Bengal. Coast. Eng. 51 (4), 277–296. Makin, V.K., Kudryavtsev, V.N., 1999. Coupled sea surface-atmosphere model: 1. Wind over waves coupling. J. Geophys. Res.: Oceans 104 (C4), 7613–7623. Makin, V.K., Kudryavtsev, V.N., 2002. Impact of dominant waves on sea drag. BoundaryLayer Meteorol. 103 (1), 83–99. Moon, I.J., Ginis, I., Hara, T., 2004. Effect of surface waves on air–sea momentum exchange. Part II: behavior of drag coefficient under tropical cyclones. J. Atmos. Sci. 61 (19), 2334–2348. Moon, I.J., Hara, T., Ginis, I., Belcher, S.E., Tolman, H.L., 2004. Effect of surface waves on air–sea momentum exchange. Part I: effect of mature and growing seas. J. Atmos. Sci. 61 (19), 2321–2333. Moon, I.J., Kwon, J.I., Lee, J.C., Shim, J.S., Kang, S.K., Kwon, S.J., 2009. Effect of the surface wind stress parameterization on the storm surge modeling. Ocean Model. 29 (2), 115–127.

17