Impact of ocean mixed‐layer depth initialization on the simulation of tropical cyclones over the Bay of Bengal using the WRF‐ARW model

The sensitivity of the simulated tropical cyclone (TC) intensity and tracks to the different ocean mixed‐layer depth (MLD) initializations is studied using coupled weather research and forecasting (WRF) and ocean mixed‐layer (OML) models. Four sets of numerical experiments are conducted for two TCs formed during the pre‐ and post‐monsoon. In the control run (CONTROL), the WRF model is initialized without coupling. In the second experiment, the WRF‐OML model is initialized by prescribing the MLD as a constant depth of 50 m (MLD‐CONST). In the third experiment, the spatial varying MLD obtained from the formulation of depth of the isothermal layer (MLD‐TEMP) is used. For the fourth experiment (MLD‐DENS), the model is initialized with the density‐based MLD obtained from ARMOR‐3D data. The results indicate that the CONTROL exhibits an early intensification phase with a faster translation movement, leading to early landfall and the production of large track deviations. The coupled OML simulations captured the deepening phase close to the observed estimates, resulting in the reduction of errors in both the vector and along the tracks of the storm. The initialization of the different estimates of the MLD in the WRF‐OML shows that the TC intensity and translation speed are sensitive to the initial representation of the MLD for the post‐monsoon storm. The gradual improvements in the intensity and translation speed of the storm with the realistic representation of the OML are mainly due to the storm‐induced cooling, which in turn alters the simulated enthalpy fluxes supplied to the TC, leading to the better representation of secondary circulation and the rapid intensification of the storm.


| INTRODUCTION
Tropical cyclones (TCs) that form over warm tropical oceans are one of the most destructive weather phenomena, with a high potential to damage coastal installations due to its extreme winds, heavy rainfall, flooding and storm surges (Emanuel, 2003;Samala et al., 2013). As the coastal regions of the North Indian Ocean (NIO) mainly of the Bay of Bengal (BOB) are highly crowded, many cities with dense populations ranging from metropolitan to cosmopolitan are vulnerable. The prediction of the movement and intensity of TCs and their associated rainfall over the BOB are of great interest, and remain a challenging problem (Raghavan and Sen Sarma, 2000;Bhaskar Rao et al., 2001;Srinivas et al., 2010). Though the prediction of the movement of TCs over the BOB has improved significantly over last few decades with the application of improved dynamic models and the assimilation of new satellite information using advanced assimilation schemes (Yesubabu et al., 2014), there is still a lag and uncertainty in the prediction of storms (e.g. Greeshma et al., 2015;Mohanty and Gopalakrishnan, 2016). One of the fundamental reasons for the failure of the existing numerical models is at recapturing the observed pattern of the TC intensities, which is possibly due to the lack of the proper representation of atmospheric ocean-coupled processes  in numerical models.
Several theories (e.g. Gray, 1975;Emanuel, 2003;Yu and McPhaden, 2011;Vissa et al., 2013) postulate that the warm oceans play a critical role in the genesis, intensification and maintenance of TCs by supplying the required energy through enthalpy fluxes. Identifying and incorporating the critical ocean-atmospheric parameters is extremely useful to improve the skill of the numerical weather models. A warmer sea surface temperature (SST), a deeper oceanic mixed-layer depth (OMLD) and the presence of anti-cyclonic eddies (positive sea level anomalies) are the prime indicators for the higher upper ocean heat (UOH) content (e.g. Vissa et al., 2012;Mohan et al., 2015) that mainly influences storm intensity. Previous researchers (Leipper and Volgenau, 1972;Sadhuram et al., 2006;Vissa et al., 2013) suggest that storms intensify through the energy fluxes at the air-sea interface over the regions of the higher UOH. Coupled global models advocate that, along with the SST, the OMLD plays a crucial role in the prediction of TC intensity (e.g. Chan et al., 2001). Mao et al. (2000) highlighted the intensification rate of the simulated storms in a coupledoceanic atmospheric model, which largely depends on the feedback of the OMLD. The OMLD acts as a source or sink in the development of a TC that depends on its deepness. The upper ocean responds during and after the passage of a TC; turbulent and vertical (diapycnal) mixing are the primary plausible mechanisms that cool the sea's surface. It can also weaken storm intensity (e.g. Vissa et al., 2012Vissa et al., , 2013. Although this negative feedback of wind-induced ocean cooling on the intensity of the storm has already been reported (e.g. Chang and Anthes, 1979;Schade and Emanuel, 1999), it has not been incorporated into storm-prediction models both regionally and globally until recently (Kim and Hong, 2010;Dare and McBride, 2011;Samala et al., 2013).
A weather research and forecasting (WRF) model is widely used for TC predictions in both research and operational applications. For the NIO, several studies (Mohanty et al., 2004;Pattanayak et al., 2008;Srinivas et al., 2012) have proven that the WRF model is capable of simulating tropical storms with reasonable accuracy. In spite of proven skill at simulating TCs across different basins, the model lacks in the implicit representation of ocean feedback as it is only supplied through the variation of SSTs supplied from lower boundaries. The oceanatmosphere coupling functionality is introduced in the WRF model (v. 3.5) with a 1D OMLD module to simulate mixed-layer depth (MLD) in an integral sense (Pollard et al., 1973) and incorporating its associate feedback into the TC. The impact of coupling the WRF with an ocean mixed-layer (WRF-OML) model is studied on the simulation of Hurricane Katrina by Davis et al. (2008), who indicates that the role of atmospheric ocean coupling successfully eliminated the spurious intensity errors found in the real-time predictions. The study by Chan et al. (2001) reveals that the impact of the MLD initialization is noticed in the case of the TC simulations using a fully 3D-coupled ocean-atmospheric model; they clearly indicate that the deeper MLD is highly favourable for the development of more intense storms at a given translation speed. The coupled model simulations of the WRFregional ocean modelling system (ROMS) by Zhao and Chan (2017) reveal that the slower translation speed and shallower MLD could negatively impact the development of TCs through the development of a symmetric cold wake and the reduction of surface energy fluxes. Most recently, Prakash and Pant (2017) showed that the coupled ocean-atmospheric model (WRF-ROMS) is capable of reproducing the observed storm-induced ocean response during the passage of Cyclone Phailin which formed over the BOB. The aforementioned studies suggested that the results of hurricanes can be improved if the realistic initial MLD is supplied to the WRF-OML.
Recently, Wang and Duan (2012) introduced a secondorder OML module in the WRF (WRF-OMLM-Noh) based on the model formulation of Mellor and Yamada (1982) and compared its performance with the existing results of the WRF-OMLM developed on a Pollard formulation (WRF-OMLM-Pollard). Their findings suggest that the performance of the OMLM-Noh is better than the existing version of a mixed-layer model in the WRF (OMLM-Pollard) in the realistic representation of a storm-induced SST, the subsequent MLD changes and the exchange of air-sea fluxes. They also attributed the improvement of the OMLM-Noh due to the realistic representation of the MLD in the model's initial conditions. The tendency of the OMLM-Pollard in the overestimation of intensities has been recognized for a while (Zhu et al., 2004), and is attributed to the underestimation of storm-induced seasurface temperature (SST) cooling and the subsequent reduction of ocean-coupling feedback in the form of fluxes in the WRF-OMLM (Yablonsky and Ginis, 2009). To replace the OML, several studies on the application of sophisticated atmosphere-ocean coupled models for the TC cases have indicated the impact of oceans on the atmosphere, and vice versa (e.g. Warner et al., 2010;Zambon et al., 2014aZambon et al., , 2014b. For the BOB, the recent TC simulation study by Prakash et al. (2018) showed the advantage of using a fully coupled atmosphere-oceanwave model (COAWST) when analysing the storminduced ocean mixing over the uncoupled ROMS model. However, the operationalization of these fully coupled three-dimensional (3D) ocean demands high computational resources. Considering these limitations, Yablonsky and Ginis (2009) and Ginis et al. (2010) suggested that a simple one-dimensional (1D) mixed-layer model is sufficient to reproduce the TC-induced SST cooling with limited computational resources. For the BOB, very few coupled ocean-atmosphere modelling studies (Mohan et al., 2015;Srinivas et al., 2016;Greeshma et al., 2019) indicate that the response of storm-induced SST cooling on the storm intensity predictions is captured by the coupled WRF-OML model. Mohan et al. (2015) attempt to analyse the impact of assimilating observations and the effect of the air-sea coupling using the WRF-OML in numerical simulation of Cyclone Nilam over the BOB. Their results also indicate that the incorporation of ocean coupling with a constant MLD initialization through the OML has an influence on the improvement of the prediction intensity of TCs, but the WRF-OML has little impact on track and rainfall predictions. Few studies in the literature study the impact of air-sea coupling in the simulation of TCs over the BOB; however, they have employed the option of a constant MLD at the time of model initialization (e.g. Mohan et al., 2015). Also, a few studies on the use of ocean-coupled models and the realistic representation of the MLD in the initial conditions over other basins (Warner et al., 2010;Zambon et al., 2014a;Prakash et al., 2018). However, the role of the initialization of realistic pattern MLDs on TC predictions has not been explored over the BOB.
The aim of the present study is to analyse the impact of the spatial variation of the MLD on the simulation of the TC track and intensity predictions over the BOB by considering two TC cases that formed during two different seasons (i.e. pre-and post-monsoon). The rationale behind considering the TCs in the pre-and post-monsoon seasons is to examine whether the changes in spatially varying the MLD can alter the results of storm simulation over the BOB. Initial versions of the WRF-OML (up to v. 3.8) work based on the initialization of a constant OMLD and deep-layer lapse rate, which is later updated (from v. 3.8) to use the climatological MLD is derived from the HYbrid Coordinated Ocean Model (HYCOM) and which works based on the depth of the isothermal layer (ILD) formulation. This performance of the WRF-OML based on the ILD formulation is reasonable over the regions (e.g. the mid-latitudes) where evaporation dominates or it equals the evaporation since the ILD coincides with the depth of the mixed layer (ML; calculated from the density). However, over the BOB where the precipitation exceeds the evaporation, the ILD exceeded the actual MLDs (Vissa et al., 2013), which prompt a search for the density-derived MLD sources rather than initializing with the ILD as the MLD in the WRF-OML model. The present study has used the MLD derived from the density criteria (ARMOR-3D level 4 data).
The paper is structured as follows. A brief history of the two cyclones selected for study is given in Section 2. Section 3 describes the model configuration, experimental design and data used for defining the SST boundary conditions. The results of the simulations are presented in Section 4. Section 5 summarizes and gives the main conclusions drawn from the study.

| DESCRIPTION OF THE HISTORY THE TWO TCs
To test the sensitivity of the OMLD initialization under different upper ocean conditions with the coupled WRF-OML, two tropical storms (Nargis and Vardah), which formed over the BOB in two different seasons (pre-and post-monsoon), were considered. Detailed reports regarding the history of the two selected TCs are available from the India Meteorological Department (IMD), Delhi (IMD, 2009(IMD, , 2017. A brief description of the cyclones is now provided. Very severe cyclonic storm Vardah formed initially as a low-pressure system over the southeastern BOB in early December 2016. The depression stage of the storm was noticed on December 6 and then intensified into a deep depression on December 7 that crossed the Andaman and Nicobar Islands. It moved northwards and intensified as a cyclonic storm on December 8. Slightly changing its direction westwards, Vardah consolidated into a severe cyclonic storm on December 9. The very severe cyclonic storm attained its peak intensity on December 11 with a central sea level pressure (CSLP) of 982 hPa and maximum sustained winds (MSWs) of 120 kmÁhr −1 . After it weakened into a severe cyclonic storm, the system landed over Chennai on December 12 and thereafter it weakened rapidly into a depression on the following day and passed over the Karnataka region and moved out into the Arabian Sea.
Cyclone Nargis was one of the devastating storms over the NIO in the decade 2000-2009. It originated as a depression over the southeast BOB in the morning of April 27, 2008, intensified into a cyclonic storm at 0000 UTC on April 28, and then into a very severe cyclonic storm at 0300 UTC on April 29. The system initially moved in a northwesterly direction and further recurved northeast and crossed the southwest coast of Myanmar near 16 N between 1200 and 1400 UTC on May 2. It maintained the intensity of very severe cyclone for about 12 hr after landfall.

| DATA, MODEL CONFIGURATION AND EXPERIMENTAL DETAILS
In the present study, an advanced research WRF (ARW) (v. 3.9.1) and a simple OML (WRF-OML) are used to analyse the impact of mixed-layer dynamics on the simulation of TCs over the BOB. The modelling system is configured into two nested domains of 27 and 9 km horizontal resolutions covering the BOB and its adjoining regions. The high-resolution inner domain covers the BOB and the east coast of India. The first domain with 27 km covers the Indian Sub-continent and NIO basin. The spatial extents of the two model domains are shown in Figure 1a. The model physics used in the study are adopted from previous modelling studies on physics sensitivity (Srinivas et al., 2013;Hari Prasad et al., 2016;Srikanth et al., 2016;Reshmi Mohan et al., 2018;Vijaya Kumari et al., 2019). The prediction of TCs over the BOB includes: the Thompson scheme (Thompson et al., 2008) for cloud microphysics; the Dudhia scheme (Dudhia, 1989) for short wave radiation; the rapid radiative transfer model (RRTM) (Mlawer et al., 1997) for long wave radiation; the Yonsei University non-local diffusion scheme (Hong et al., 2006)   turbulence; the Kain-Fritsch mass flux scheme (Kain, 2004) for cumulus convection; and the Noah land surface model scheme for representing land surface processes (Chen and Dudhia, 2001). The initial conditions are provided by Global Forecast System (GFS) analysis available at a 0.5 × 0.5 horizontal resolution from the National Center for Environmental Prediction (NCEP). Lateral and lower boundaries are updated at 6 h intervals. Even though the GFS analysis is available at 3 h intervals, the lower initial boundary conditions (SST) are available at 6 h time resolutions. To maintain same time resolution for all the boundary conditions, the boundaries were updated with the GFS forecasts at 6 h intervals. For the WRF-OML system, the ocean initial conditions were obtained from the global reanalysed data of the HYCOM (https://www.hycom.org/data/glba0pt08/expt-91pt2). The WRF-OML is a conventional bulk OML model based on the formulation of Pollard et al. (1973). It incorporates the feedback obtained from the windinduced ocean mixing and deepening of the mixed layer to atmospheric model in the form of SST changes. The OML model is introduced by Davis et al. (2008) in the WRF, and is designed for the simulation of the mixedlayer process in an integral sense. Using the wind stress supplied to the turbulent mixed layer, the model simulates currents in the OML, which results in the mixing of colder temperature with the deeper layers. The entrainment of colder water from the deeper layers to the top changes the surface temperatures and cools the OML, and also, in turn, changes the SST. This process is carried out only through the vertical redistribution of temperatures as the OML assumes that no heat transfer takes place across the grid points. The effect of the Coriolis term is included in the OML to incorporate the rotation of inertial currents and its associated impact on mixing and SST changes. The model neglects the pressure gradients and horizontal advection terms. The absence of the horizontal advection and upwelling terms in the OML formulation induces uncertainties of about 15% in the simulated storm SST cooling (Price, 1981). For faster moving storms (translation speed of TCs > 4 mÁs −1 ), the horizontal advection term has a little impact on the simulation of storm-induced SST cooling (Yablonsky and Ginis, 2009). Since the WRF-OML model computes the changes in the SST mainly based on the initial fields of the MLD, deep-layer temperature lapse rate and surface wind stress at the top, Davis et al. (2008) pointed out that the accuracy of the WRF-OML depends greatly on the specification of the initial values of the ocean's thermal state. Also, they suggested that the spatially varying and realistic initial MLD values needed to obtain reliable estimates of storm-induced SST variations.
The OML basically computes the simulated depth of the mixed layer (dh/dt). For an arbitrary time solution (n), the rate of change in depth of mixing is given by: where γ is a constant; u is the velocity; and R i is the Richardson number: Surface ocean temperature changes can be derived from the conservation of heat: where h is the MLD; T is temperature; Γ is the temperature lapse rate; c is the heat capacity per unit mass; ρ is seawater density; g is gravity; α is the thermal co-efficient of expansion; u and v are horizontal velocities; Δu is mean velocity; and Q is heat flux down through the surface. Four sets of numerical simulations are performed for the two TC cases by initializing the model at 1200 UTC on April 28, 2008, and 0000 UTC on December 9, 2015, respectively, and integrated up to 120 hr for Nargis and Vardah, respectively. In the first experiment (CONTROL), along with GFS analysis, the lower boundary SST fields are replaced by time-varying fields of high-resolution real-time global (RTG) SST, and another three experiments are conducted by incorporating the coupling feedback of the OML physics through the WRF-OML system for each TC. This feedback in the WRF-OML configured at every model time step (i.e. 60 s) where the WRF model supplies surface wind stresses (τ) and heat flux (moisture, latent and sensible) to the OML, while the OML intern provides the storm-induced response in the form updated SSTs to the WRF (Davis et al., 2008;Mohan et al., 2015). Of the three WRF-OML experiments, the first is configured with a constant MLD of 50 m (MLD-CONST). In the second (MLD-TEMP), the coupled system is initialized with the spatial variation of the MLD computed based on the isothermal temperature profiles of the HYCOM data set and the lapse rate for the OML configured at −0.14 CÁm −1 . In the third (MLD-DENS), the initial values of the MLD are replaced with density-based MLDs obtained from the three-dimensional (3D) thermohaline field (ARMOR-3D level 4 data). Though the density-based MLD can be derived from the global HYCOM reanalysed data set available because it has a high horizontal resolution (about 0.08 ), it has a lesser number of vertical levels (32). In particular, the major limitation arises from upper ocean surface levels as the global data contain only seven layers within the upper 100 m of the ocean's surface. Therefore, the spatial MLDs obtained from the ARMOR-3D in the MLD-DENS experiment were used. The MLD product of the ARMOR-3D was derived from the density formulation by using the high-resolution density profiles of the Argo data. The ARMOR-3D MLD is computed based on the variable density criteria, that is, the depth at which density has varied significantly from the surface value by a threshold value of 0.2 C (de Boyer, 2004). For more details about the density-based MLD data set, see Buongiorno Nardelli et al. (2017). Among the experiments (MLD-CONST, the MLD-TEMP and MLD-DENS), the present authors only changed the spatially varying MLDs while maintaining the same physics and domain configurations expected in the OML, and it was assumed that the simulated changes in track and intensity were primarily observed due to the combined effect of the initial MLD and the defined initial profile of temperature and density.
To validate the model results, the intensity parameters and track information for the two cyclones were obtained from the cyclone reports of the IMD (2009,2017). The azimuthally averaged radial-height crosssections of the temperature anomaly and tangential winds, which provide a quantitative structure of the cyclone from the model simulations, are compared with the Cooperative Institute for Research on Atmosphere (CIRA) products of radial height cross-sections of Advanced Microwave Sounding Unit (AMSU) profiles and multi-platform tropical cyclone analysis products (Knaff et al., 2011). The surface wind analysis from the CIRA is available at a horizontal resolution of 10 km and at 6 hr intervals. The gridded rainfall estimates from the IMD and Tropical Rainfall Measuring Mission (TRMM) 3B42v7 satellite-merged precipitation data are used for a comparison with the model-simulated rainfall.

| RESULTS AND CONCLUSIONS
The results are organized into three subsections. The time variation of the simulated storm tracks and intensity that have already been discussed is first presented. Moving further, the response of the coupled model at simulating the TC-induced features over the upper ocean, such as simulated changes in the OML and a storminduced SST cold wake, are then discussed. Finally, the impact of ocean feedback on the secondary circulation, organization of winds and rainfall of TCs are examined. All the results for the two TCs (Vardah and Nargis) are presented from the high-resolution (9 km) model domain.

| Tracks and intensity of the simulated storms
The model-simulated tracks for the two cyclones along with the IMD best track estimates are shown in Figure 1b,c. The best track estimates of the IMD indicate a westward movement of the system in the case of the post-monsoon Cyclone Vardah, while an eastward movement is observed in the case of pre-monsoon storm Nargis. All the model simulations captured the pattern of the observed track during the initial 36-48 hr of simulation. Thereafter, the differences in the vector track are clearly visible between the coupled and control simulations. For Nargis, though the model is initialized at 1200 UTC on April 28, 2008, the best track data are available only between 1200 UTC on April 29 and 0000 UTC on May 3. The initial track deviations of the storm from the IMD in all experiments could be due to the dislocation of the initial vortex from that observed. The role of the WRF-OML coupling on track predictions is seen only after 42 hr of model simulation in both TC cases. This slow response is probably due to the limitation of the OML model discussed in several studies (Zhu et al., 2004;Halliwell et al., 2008). It is primarily due to the simulated OML response, and the resultant TC-induced cooling strongly depends on both initial MLDs as well the characteristics of the TCs (translation speeds and sizes of the storms) passing over the basin. The initial MLD values (Supporting Information Figure S1a-c) are not changing significantly (< 10 m). This could be the primary reason for the slow response in the track simulations for Nargis. The results of the coupled model simulations (Figure 1b) shifted to southwards compared with the uncoupled simulations (CONTROL) and moved closer to the observed storm positions during landfall. The track variations in the WRF-OML are due to variations in the simulated mixed layer, which induce the differences in the simulated SST cold wakes in the right forward sectors and shift the model storm to move towards the warmer SSTs (leftward forward sectors). This southward movement of the coupled model compared with uncoupled simulations has been discussed by Srinivas et al. (2017). In the present study, the translation speeds (kmÁ6 hr -1 ) are calculated for the two cyclones as a measure of distance travelled by the model storm in 6 hr intervals. The translation distances (in 6 hr) for the CONTROL, the MLD-CONST, the MLD-TEMP and MLD-DENS are 102.4, 93.5, 95.8 and 91.3 km, respectively, for Nargis, and 113.9, 111.5, 108.3 and 107.9 km, respectively, for Vardah. In the case of the translation speeds (in 6 hr) for the observed storm computed from the best estimate of the IMD are 89.7 and 86.5 km for Nargis and Vardah, respectively. Among all the experiments, the deviations in the track are great with the CONTROL and the vector track spread between the WRF-OML experiments is relatively less. Apart from the vector track deviations, the translation speed of the storm is fast in all the experiments, especially with the CONTROL in both TCs. Of the coupled model simulations, the MLD-DENS simulated tracks are closer to the observed and its translation speed improved when the model storm encountered relatively shallow MLDs.
The time series of vector track error (VTE) computed with the IMD estimates for both TCs are shown in Figure 2a,b. The track errors indicate that the deviations of the CONTROL experiments are increasing with forecast length at high rate as compared with the WRF-OML runs in both TC cases. The response of mixed-layer coupling on the VTE is seen only after 30-36 hr of simulation for both TC cases; the impact and spread of the VTE appear to be higher for the post-as compared with the pre-monsoon. This coincides with a previous discussion on the slow response of the OML (Halliwell et al. 2008) due to meagre variations in the initial condition of the MLD. The coupled model simulations consistently simulated westward movement for Vardah and eastward movement for Nargis closer to the observed storm till landfall, while uncoupled simulations (CONTROL) moved northward during the mature phase of the storm. Also, the errors (VTE) are not only due to the cross-track deviations of the model, as it arises along the track deviations of the model from the variations in the translation speed of the model storm. Thus, the translation speed and VTE clearly suggest that experiment MLD-DENS has a relatively better simulation, followed by the MLD-TEMP and MLD-CONST. The present results suggest that in a basin such as the BOB, where precipitation exceeds evaporation, the variations in the initial MLD obtained from temperature and density criteria will have a significant impact on the vector track variations of the simulated storm.
F I G U R E 2 Time variation of the vector track error (VTE) (km) (a, b), simulated central sea level pressure (CSLP) (hPa) (c, d) and maximum wind (mÁs −1 ) (e, f), for Nargis (left) and Vardah (right). The arrows on the x-axes indicate the landfall points in the respective simulation Time variation of the simulated intensity parameters' CSLP and MSW along with the best estimates of the IMD are shown in Figure 2c-f. For Nargis, the observed estimates indicate that the lowest pressure drops and MSW are about 962 hPa and 47 mÁs −1 , respectively, which occurred for a period of 6 hr. Though all the simulations underestimated the pressure drop, the CONTROL experiment predicted an early life cycle of the storm with a high intensification rate. It attained maximum intensity before 30 hr of the observed, and the lowest pressure simulated by the CONTROL experiment was 968 hPa. Even though the OML experiments underestimated the intensity by nearly 8 hPa, resulting in weaker winds compared with the IMD and CONTROL, the simulations seem to produce a gradual deepening and mature phases of the storm close to the timings of the observed storm. The significant underestimation in the pressure drop between the model and the observed, starting from the predeepening to weakening phases of the storm, is probably due to the cold start initialization of the WRF model without the assimilation of conventional and satellitederived data. Previous studies on the simulation of Nargis (Osuri et al., 2012;Srinivas et al., 2012) suggested that the assimilation of satellite-derived winds using the mesoscale modelling system is need to improve the intensity predictions of Nargis. Moreover, in the coupled ocean atmospheric model, the OML component, which is a simple slab ocean model that does not account for all the ocean-wave-atmospheric interactions such as the wind wave-induced positive feedback on storm intensity (Liu et al., 2011), was used. Therefore, considering these limitations, Osuri et al. (2013) highlighted the need for a fully coupled ocean-atmospheric assimilation system for the better prediction of TCs over the NIO. The time variations of the MSW in all the experiments of Nargis exhibit slight variations in attaining the maximum intensity (by 12 hr) compared with the CSLP. The maximum surface winds simulated in the coupled model are about 34, 36 and 35 mÁs −1 at 1200 UTC on May 1, 2008, by the MLD-CONST, the MLD-TEMP and MLD-DENS, respectively. Moreover, the simulated maximum surface winds in all the experiments are underestimated by 10-15 mÁs −1 . As discussed in the CSLP analysis, the underestimation of intensity by the model is probably due to the lack of an initial vortex at the time of initialization, and also mesoscale model tendency, as indicated by several authors (Trivedi et al., 2002;Osuri et al., 2012;Srinivas et al., 2012), which needs to be addressed by increasing model resolution, with the proper representation of the initial vortex and assimilation of all available observations. Thus, the intensity estimates of Nargis suggests that the impact of initialization of the spatial MLD variations is clearly seen on simulation of the TC's intensity. However, the intensity differences between the mixed-layer experiments are less sensitive for Nargis, particularly with the spatial variations of the MLD derived from density and temperature.
In the case of Vardah, the IMD estimates showed a gradual deepening of the storm for a period of 36 hr and it attained its lowest pressure drop of 975 hPa and an MSW in the order of 37 mÁs −1 in its mature phase, which occurred between 0000 and 0600 UTC on December 12, 2016. The CONTROL experiment rapidly intensified and attained maximum intensities of 978 hPa and 33 mÁs −1 early by 12 hr compared with the observed estimates. The intensity differences between the OML experiments were less sensitive in the first 24 hr of the model integration, and the variations in intensity started to appear after 30 hr of simulation. The differences in both the intensity and time of attaining maximum intensity were seen when the model reached its mature stage. As discussed in the track simulations, the slow response of the OML (Zhu et al., 2004;Halliwell et al., 2008) is probably the major reason for this time lag of 30 hr. Zhu et al. (2004) reported similar variations in the coupled OML simulations of mesoscale model 5 and their study attributed it as a response of the OML. Further, Zambon et al. (2014b) suggested that the coupling of the wave component (simulating waves nearshore) with the atmospheric model (WRF) can only provide the required feedback for the co-efficients of surface roughness lengths to the atmospheric model (WRF) to produce a significant difference in simulations of the maximum wind speed and pressure drop just after initialization of the WRF model (as seen in the analysis of Hurricane Sandy). Of all three OML experiments, the MLD-CONST attained maximum intensity first by 12 hr before the MLD-DENS and 18 hr earlier compared with the MLD-TEMP. Though the MLD-DENS captured the observed intensity pattern during the deepening phase of the storm, it attained the observed maximum intensity of 975 hPa by 18 hr before the observed time. The MLD-TEMP simulated a slowly intensified storm and attained the observed maximum intensity (975 hPa) after the 60 hr of model simulation. Interestingly, the time of attaining maximum intensity is close to the observed estimates by 12 hr. However, the intensity of the storm is overestimated by the MLD-TEMP during the deepening phase between 36 and 54 hr of model simulation. Early intensification of the MLD-DENS compared with the MLD-TEMP is possibly due to the specification of lower depths obtained from density profiles which suppressed the SST-induced cooling and enhanced the flux feedback from the ocean, which favoured the rapid intensification of the storm. The incorporation of the spatial MLD on TC prediction is clearly visible in the case of post-monsoon storm Vardah, as the intensity difference are seen from the deepening phase of the storm and it could be due to the association of strong winds and a deepening of the MLD, which is examined in the next section.

| Storm-induced ocean response
In this section, the simulated TC-induced ocean response from the three WRF-OML experiments is presented in terms of the simulated spatial variation of the MLD, storm-induced SST cooling and flux feedback. The WRF-OML estimates the MLD after computing the rate of change in the depth of mixing, which is directly proportional to surface wind speed. The initial MLDs (Supporting Information Figure S1) supplied to the WRF-OML indicate that there are minor differences in the spatial distribution of the MLD computed from the formulation of isothermal temperature and density in the two TC cases. The deviations of the spatial MLD are large in the case of the post-monsoon TC (Vardah) as compared with pre-monsoon storm Nargis.
To illustrate how the initial differences in the spatial distribution of the MLD between the MLD-TEMP and MLD-DENS are directly reflected in the storm-induced deepening, the spatial distribution of the simulated MLD along with the corresponding storm tracks are shown in Figure 3. In the case of the MLD-CONST, a maximum MLD deepening of about 5-10 m depth is only noticed during its peak intensification period, while in the MLD-TEMP, there is a clear signal of the presence of a deeper MLD up to 15-20 m along the track, and it gradually decreases from the right side to the left side of the track. A similar pattern is also seen in the MLD-DENS with a high magnitude of deepening up to 20 m. The results of the simulated response of the OML suggest that the deepening of the MLD clearly occurred on the right side of the storm's passage for both TCs (Supporting Information Figure S2). The results suggest that the simulated deepening of the OML is highly sensitive to the initial MLDs supplied to the WRF. Also, the simulated deepening is clearly seen for the MLD-DENS, and is found to be high for Vardah. The impact of deepening is found to be less (5-10 m) for the MLD-CONST.
Previous studies (Maeda, 1964;Gentry, 1970) showed that when an intensified cyclone passes through a shallow mixed-layer, storm-induced strong winds lead to the upwelling of the ocean waters from the subsurface to cool the surface ocean waters. It leads to a decrease in the SST during the passage of the cyclone, which in turn acts like a negative feedback of the ocean to reduce the flux feedback to suppress the intensification of the TC (Mao et al., 2000). To present the storm-induced SST cold wake, the simulated SST difference (Figure 4) computed between the SSTs at the time of model initialization and after the passage of the TC, along with the corresponding microwave SST obtained from a special sensor microwave imager (SSMI), is presented. The SSMI SSTs clearly indicating the storm-induced cooling (> 1.5 C) for both TCs on the right side of the cyclone track. While the results of the simulated cold wake in both TC cases indicate that when the spatial variation of the MLD is set to constant (50 m) in the initial conditions (MLD-CONST), though the SST cooling is observed along the right side of the storm track, the magnitude is found to be less (0.25-0.5 C) compared with the other WRF-OML simulation. The magnitude of the storm-induced SST cooling is slightly increased from 0.5 to 0.75 C when the initial MLD is obtained spatially from the isothermal layer depth (MLD-TEMP). Though the simulated MLD of the MLD-TEMP is deeper than the MLD-CONST (Figure 3), the storm-induced SST cooling is high in the case of the MLD-TEMP, and is because the initial MLD is shallower in the MLD-TEMP compared with the MLD-CONST (Supporting Information Figure S1). The initialization of spatially varying the MLD produced the significant changes in the simulated life cycle of the storm and resulted in the simulated strong surface winds (than the MLD-CONST) that resulted in mixing with deeper layers (Supporting Information Figure S2b) and might have produced the strong cooling conditions as reflected in the form of SST cooling. The magnitude of the storm cooling (> 1.25 C) is significantly increased and a well-marked pattern is observed on the right side of the cyclone track in the simulation of the MLD computed from density (MLD-DENS). Simulation of the storm-induced SST is one of the prominent parameters for the coupled model simulations for TCs. The results indicated that the simulated SST cold wake obtained from the MLD-DENS realistically captured the observed pattern of the SSMI SST product in both cyclone cases. However, there is an underestimation of storm-induced SST cooling (0.5-1 C) even in the best simulation (MLD-DENS); it is probably due to the limitation of the 1D OML model highlighted by Yablonsky and Ginis (2009). Their study shows that the OML (initialized with deeper MLDs) has a tendency to underestimate significantly (more than half) the storm-induced cooling (SST), particularly for the TCs translating at speeds of 3.5 and 5.0 mÁs −1 . Also, their study attributed the underestimation in ocean surface cooling primarily to neglecting the process upwelling mechanism. Kossin (2018) indicates that the average translation speed for NIO storms is in the range of 3.5-4 mÁs −1 . Owing to these factors, the underestimation of the simulated storm-induced cooling can be mainly attributed to the limitation of the OML model, which can be rectified by coupling with a 3D ocean model (Zambon et al., 2014b;Prakash and Pant, 2017).
The ML acts like a source enhancing or reducing the intensity of the storm passing over, and the upper ocean feedback to the TC atmosphere is primarily through the exchange of enthalpy fluxes (latent and sensible heat). To illustrate the impact of the simulated MLD changes on the surface fluxes, a time series of enthalpy fluxes averaged over the region of 4 × 4 around the TC centre are presented in Figure 5. The results show that there are significant variations in the magnitude of fluxes between the coupled and uncoupled models. Also, the variation of the sensible heat flux is large for the pre-monsoon storm, as warm summer temperatures could possibly enhance the sensible flux feedback in the coupled model simulations of Nargis, while the major differences in latent heat flux are seen after the storm reaches its highest intensity in both cases. Srinivas et al. (2017) also point out similar variations in the flux feedback in the TC simulations between the coupled and uncoupled models, and their study explained that a lack of formulation of wind-induced mixing in uncoupled models is a possible reason for these flux variations. Among the coupled model simulations, significant differences were found mainly in the mature and weakening stages of the storm. To illustrate the precise differences in the simulated fluxes, Supporting Information Figure S3 shows the time variation of area-averaged latent heat and sensible heat fluxes over the four quadrants around the centre of the TC (left forward, right forward, right rear and left rear) and in the quadrants determined based on the simulated track position of the storm for all time intervals. Maximum variations in the simulated fluxes are seen in the left quadrant of the storm as maximum winds in the left sectors of the storm, which lead to a reduction of latent heat flux, and increases in sensible heat fluxes are seen mainly in the left forward and rear sectors. Significant differences are found in the coupled simulations, mainly from the deepening to the mature stages of the storm. The CONTROL simulation exhibits minimum variations in sensible heat flux, while they are maximum in the case of the simulated latent heat fluxes in all quadrants, and this is possible due to the reduction of storm-induced cooling in the CONTROL, which enhanced surface moisture fluxes and reduced sensible flux transfer from the ocean to the atmosphere. In coupled OML F I G U R E 5 Time series of the simulated area-averaged sensible (top) and latent heat fluxes (bottom) around the cyclone grid for Nargis (left) and Vardah (right) simulations, the maximum fluxes are simulated by the MLD-CONS, while the minimum transfer of fluxes in all quadrants is seen in the case of the MLD-DENS, and it is mainly due to the simulation of enhanced storm-induced cooling in the MLD-DENS, which reduced the transfer of fluxes to the WRF.

| Surface winds simulated
In this section, the simulated surface winds are analysed with the multi-platform winds from the CIRA at the mature stage of the cyclonic storms considered in the study. The vector track positions clearly indicated that there are significant differences in translation speed between the coupled and uncoupled models. For Nargis, the CIRA winds ( Figure 6) are considered from 0000 UTC on May 2, 2008, while the simulated winds are considered at 1800 UTC on May 1. The observed estimates from the CIRA indicated strong winds of magnitude 65-80 kn, which are confined in the right rear sector of the cyclonic storm, and the surface winds of magnitude 35 kn are seen over wider regions from 13 to 18 N with a circularly symmetry around the centre of the TC positioned at 93 E and 16 N. The observed pattern of isotachs around the centre of the storm and the spatial distribution of maximum winds (> 35 kn) are well captured by the MLD-DENS and MLD-TEMP simulations, with simulated peak winds of 65 kn in the right forward sector of the TC. Moreover, the position of the TC centre and spatial extent of the maximum winds simulated by the MLD-DENS are close to the observed estimates of the CIRA. The MLD-CONST simulated slightly weaker winds with speeds around 55 kn. Also the organization of the isotachs is slightly shifted to left forward sectors; and the position of the system is close to the coast of Myanmar. Though the similar patterns of weaker winds (45 kn) confined over a smaller region is noticed in the CONTROL, the position of isotachs is shifted to rear sectors of the TC and the centre of the TC is dislocated from the observed by around 100 km. In the case of Vardah, the CIRA observation is considered at 1800 UTC on December 11, 2016, when the storm attained its mature stage and is about to make landfall, while the simulated winds are shown at 0600 UTC on December 11, 2016. The CIRA estimates clearly indicate strong winds of > 50 kn in the right rear sector of the TC. The spatial distribution of winds and magnitude are well captured in the MLD-TEMP and MLD-DENS, while the CONTROL and MLD-CONST predicted a weaker storm with a lesser magnitude of winds (< 50 kn), and the spatial extent of maximum winds (50 kn isotach) is confined over a smaller region. Although the experiments initialized with spatially varying MLD-simulated surface winds closely match the observed maximum winds' (50 kn) pattern and locations, the spatial extent of the isotachs is slightly overestimated in these experiments (MLD-TEMP and MLD-DENSE). This is probably due to the change in the simulated life cycle of the storm after the initialization with a spatially varying MLD. The IMD estimates suggest that the observed storm attained its maximum intensity between 0000 and 0600 UTC on December 12, 2016, for Vardah. However, the experiments CONTROL, the MLD-CONST, the MLD-TEMP and MLD-DENS showed early intensification and attained their mature stage, respectively, at 1800 UTC on December 10, 0000 UTC on December 11, 1200 UTC on December 11 and 1200 UTC on December 11, 2016. As the surface wind comparison time (1800 UTC on December 11) is close to the mature stage of the MLD-TEMP and MLD-DENS, the surface winds were slightly greater in magnitude in the two experiments.

| Radius height wind and temperature of Nargis
To illustrate the impact of mixed-layer coupling on the simulation of the secondary circulation of the TC, radius height distributions of azimuthally averaged tangential winds are plotted (Figure 7a-e) for Nargis at the mature stage of the TC (0000 UTC on May 1, 2008). The observed estimates of gradient wind indicate calm winds (< 10 kn) over a radius of 20 km from the centre of the TC and strong winds (40 kn) that spread vertically up to 10 km with a radial distance of 80-180 km; further the magnitude of the gradient winds decreased significantly in both vertical and horizontal directions. While the simulations show the strongest gradient winds observed when the model is initialized with the MLD-DENS, the models also overestimated the winds. The results of the coupled and uncoupled model show that the width of the calm wind region and magnitude of the maximum gradient wind are almost similar and close to the observed estimates. However, the radial distance and vertical extensions of the strong gradient wind distribution (50 kn) differ in all four simulations. The CONTROL simulations exhibited a weaker secondary circulation with maximum winds (40 kn) extended up to 4 km and radially by 60 km; similar results are also noticed in the MLD-CONST with a maximum wind confined to a very small radius of 40 km. The results of secondary circulation are improved with the MLD-TEMP with the strongest winds (40-50 kn) extended vertically and horizontally up to 5 and 80 km, respectively. The better comparisons found with the MLD-DENS indicate the presence of strongest winds up to 6 km in a radius of 100 km from the centre of the storm.
Similarly, the radial height section of the azimuthally averaged temperature anomaly for Nargis at the mature stage of the storm (0000 UTC on May 1, 2008) is presented in Figure 7f-j. The observed temperature anomaly is available directly from the CIRA, whereas the corresponding simulated temperature anomaly is computed as a difference in temperature computed from the mature stage of the TC to the temperature before the formation of the storm. The observed estimates of the temperature anomaly show a cooling at lower levels up to 6 km and a warming region (3 C) found between 10 and 15 km from the centre of the storm to a radial distance of 50-100 km. The observed estimates indicate a strong convergence and associated cooling due to evaporation at lower levels and the existence of a strong outflow region and the subsidence of air, resulting in warming in the mid-Troposphere. Though the presence of upper tropospheric warming was captured in all experiments, the model failed to simulate the lower tropospheric cooling. Also, the upper tropospheric warming is overestimated (1-2 C) by the model compared with the observed CIRA estimates. For the better pattern of upper

| Simulated rainfall
The spatial distribution of 24 hr-accumulated rainfall both for Nargis from 0000 UTC on May 1 to 0000 UTC on May 2, 2008, and Vardah between 0000 UTC on December 12 and 0000 UTC on December 13, 2016, are compared with the observed rainfall estimates of the TRMM 3B42v7 ( Figure 8). In the case of Nargis, the model overpredicted the rainfall in all experiments, with a slight displacement of the rainfall band probably due to faster translation movement as compared with the observed estimates and the maximum overestimation found in northward regions of the track. The TRMM rainfall pattern indicates a maximum rainfall distributed in southern parts of the storm track. The uncoupled model simulated the highest rainfall (> 200 mm) with maximum rainfall distributed mainly over the north-central BOB, and simulated rainfall reduced as the storm moved close to the coast of Myanmar. The simulated pattern of rainfall from the coupled model simulations indicates that the simulated rainfall with the WRF-OML is sensitive to the initial MLDs. Among the coupled model simulations, the MLD-DENS is close to the spatial distribution of the TRMM rainfall and the overestimation of rainfall significantly reduced when compared with the MLD-CONST and MLD-TEMP. The rainfall in these simulations (MLD-CONST and MLD-TEMP) was significantly overestimated in the northern sectors close to Myanmar due to faster translation movement of the simulated storm. In the case of Vardah, all the experiments slightly underestimated the rainfall and failed to capture the observed symmetric distribution of rainfall, particularly that over southern parts of the storm track seen in the TRMM. Though the spatial pattern of the simulated rainfall was close to the observed estimates, the heavy rainfall band completely shifted northwards by 100 km due to the northwesterly movement of the storm track in the CONTROL. Among the coupled model simulations, the MLD-DENS produced a better comparison with the observed estimates along both the right and left sectors of the track, and the rainfall peak over the coast of Chennai. Thus, the role of the initialization of the MLD on the simulated TC rainfall is seen in both pre-and postmonsoon storms, but the impact seems to be more significant for post-monsoon storm Vardah. Moreover, the rainfall simulation with the MLD-DENS seems to produce a realistic pattern of rainfall compared with other coupled model simulations, suggesting that the initialization of the MLD obtained from the density criteria has a positive response on the simulation of TCs' rainfall.

| SUMMARY AND CONCLUSIONS
In the present study, the response of the simulated tropical cyclone (TC) intensity and track with respect to the initialization of the variable oceanic mixed-layer depth (OMLD) analysed for the two TCs that formed during the pre-and post-monsoon using the ocean mixed-layer (OML) model coupled with the weather research and forecasting (WRF) model has been analysed. The impact of the OML processes on the prediction of the TC intensity and track that formed over the Bay of Bengal (BOB) is also analysed using the best track estimates of the India Meteorological Department (IMD).
Of the four numerical experiments, the differences in the simulated track are seen mainly between the CON-TROL and WRF-OML simulations, and the translation speed of the storm is found to be faster in the CONTROL simulation. As for the coupled WRF-OML simulations, the translation speeds and vector track error (VTE) clearly indicate that the cross-track variations among the coupled model simulations are less, and the differences are seen mainly along the track due to changes in the translation movement of the storm, particularly for post-monsoon storm Nargis. Of all the simulations, the MLD-DENS is relatively better simulated, followed by the MLD-TEMP and MLD-CONST, suggesting the dynamic initialization of the MLD obtained from density criteria seems to produce positive results for the track simulations, particularly for postmonsoon storm Vardah. The intensity differences between the mixed-layer experiments are less sensitive for premonsoon storm Nargis; however, significant differences in the intensity are seen for post-monsoon Vardah. In the case of the latter, the differences in intensity are seen from the deepening phase of the storm. Also, early intensification is noticed in the case of the MLD-DENS, which could be attributed to the specification of lower depths in the MLD-DENS that have reduced the sea surface temperature (SST) cold wake and enhanced the flux feedback from the ocean, leading to the simulation of the enhanced intensification of the storm. Moreover, saline stratification and the presence of thick barrier layer thickness in the BOB during the postmonsoon season are also attributed to the reduced TC SST cooling (e.g. Balaguru et al., 2014;Foltz et al., 2018).
To check the major factors that play a critical role in the simulated differences in track and intensity, the results of the ocean response in terms of variations in the simulated MLD and storm-induced SST cooling and its cumulative effect on the simulated ocean fluxes are analysed. The results of the simulated MLD and storm-induced SST suggest that the deepening of storm-induced MLD and SST cooling is primarily observed on the right side of the storm passage, and the storm-induced deepening is high for the MLD-DENS, in particular for Vardah. The results of surface enthalpy fluxes indicate that significant variations exist in the magnitude of fluxes, mainly between the coupled and uncoupled models, and the variation in the sensible heat flux are large for the pre-monsoon storm, while the major differences in latent heat flux are seen in post-monsoon storm Vardah. These differences in fluxes with the numerical simulation of the WRF-OML seem to produce a better pattern for the simulated surface winds and secondary circulation of the TC, while there is also consistent improvement seen from the initialization of the constant MLD to the spatially varying MLD derived from temperature and density. A comparison of Cooperative Institute for Research on Atmosphere (CIRA) surface winds and the secondary circulation of the TC are made through the radial-height sections gradient wind and temperature anomaly. The results provide evidence for the better agreement with the MLD-DENS. Though the rainfall simulations from all the experiments are not close to the observed estimates of the Tropical Rainfall Measuring Mission (TRMM), the initialization of the spatially varying MLD shows a positive impact on the simulation of TCinduced rainfall, especially for post-monsoon storms compared with the pre-monsoon. Among the three different MLD initializations, the positive impact on rainfall is observed with the initialization where the MLD is computed from density. This suggests that the dynamic initialization of the MLD obtained from density criteria better simulates the post-monsoon storms.