Appraisal of Data Assimilation Techniques for Dynamical Downscaling of the Structure and Intensity of Tropical Cyclones

The dynamical downscaling technique is used for the understanding of physical mechanisms associated with the atmospheric phenomena. We have developed high‐resolution analysis (6 km) for three tropical cyclones (TCs), namely, Phailin (2013), Nilofar (2014), and Chapala (2015) originated over the North Indian Ocean using the dynamical downscaling approach. The study aimed at the identification of appropriate methodology for generating analysis so that it becomes useful for identifying the role of environmental and internal dynamics on intensification processes and structural changes of TCs. The simulations using Weather Research and Forecasting model and four‐dimensional variational (4DVAR), hybrid three‐dimensional ensemble‐variational (3DEnVAR), and a hybrid four‐dimensional ensemble‐variational (4DEnVAR) data assimilation (DA) techniques are compared. The impact of DA is quantified by comparing errors in position, minimum sea level pressure, and maximum wind speed with the best track data set of India Meteorological Department. The intensities of TCs simulated by three downscaling methods are validated in terms of changes in minimum sea level pressure, maximum surface winds, and boundary layer and middle tropospheric relative humidity. The skills scores, namely, equitable threat score, false alarm ratio, the probability of detection, and biases (BIAS), are calculated to identify the best suitable DA technique. It is found that the hybrid DA techniques improve the overall quality of analysis compared to those developed using only variational DA techniques. The simulation using the hybrid 4DEnVAR DA technique is found to be better for simulation of the track, intensity changes, and structural characteristics of TCs.


Introduction
Owing to the improvement in numerical weather prediction (NWP) models, availability of satellite observations especially over the regions having sparse network of in situ weather observations and development in advanced data assimilation (DA) techniques, the simulation accuracy and forecasting of tropical cyclones (TCs) features have improved to a great extent in the last decade. However, according to the author's knowledge, the simulation of changes in the intensities of TCs is still a challenge and needs to be attempted. There are major concerns about the representation of intensities of the TCs in the reanalysis data sets. Schenkel and Hart (2012) showed that some of the reanalysis uses vortex relocation, which in turn responsible for the simulation of the nonphysical life cycle. In contrast, Gao et al. (2018) indicated that a vortex relocation/bogusing, assimilation of TC observation data, and model parameterizations are essential to adequately represent warm cores in terms of strength and three-dimensional structure. Both of these studies indicated the requirements of more observations around the TC core. The high-resolution reanalysis developed using mesoscale models and dynamical downscaling techniques are being used to understand the physical mechanisms associated with the structure and intensity changes associated with TCs (Hendricks et al., 2004;Rajasree et al., 2016;and Rajasree et al., 2017). However, the conflicting results from above these studies indicated that it is essential to find out which DA methodology is most suited for the development of reanalysis, which represents the structure and intensities of TCs appropriately.
The dynamical downscaling techniques are found to be useful for the development of high-resolution reanalysis. These reanalysis data sets represent the changes in the structure and intensities of TCs with reasonable accuracy and therefore are found to be helpful for the understanding of the physical mechanisms associated with the evolution of the TCs (Hendricks et al., 2004;Rajasree et al., 2016;Rajasree et al., 2017). Many DA strategies are being used for developing high-resolution reanalysis, and it is important to test which approach is most appropriate to simulate the structure and intensity of TCs. The three-dimensional variational (3DVAR) DA technique uses nearly homogeneous, quasistatic, time-invariant flow-independent background forecast error covariance (BEC). However, in reality, the forecast error covariance may vary with the flow of the system. Due to the vortical and nongeostrophic motions in TCs, the BEC of model forecasts depends on the characteristics of flow. Therefore, quasistatic BEC used in the 3DVAR system has limitations in simulating the evolutions of TCs very accurately. Although a four-dimensional variational (4DVAR) DA scheme implicitly uses a time-evolving covariance model under linear tangent dynamics, it uses a static covariance model at the beginning of each 4DVAR cycle. These limit the accuracy of simulation of the evolution of TCs using 3DVAR and 4DVAR techniques and may limit the understanding of the physical mechanisms associated with the development of TCs.
The background-error covariances in ensemble Kalman filter (EnKF) are calculated using an ensemble of the short-range forecast. Therefore, BEC used in EnKF is flow dependent and consequently found to be better compared to 3DVAR/4DVAR techniques. Recent studies (Fairbairn et al., 2014;Hamill & Whitaker, 2005) suggested that an ensemble-based DA system can produce promising results in weather forecasts in NWP models. However, it needs a large number of ensembles to obtain excellent spread required to represent the system more realistically. In a hybrid DA approach, the coupling of ensemble and variational frameworks is used. It uses an ensemble-based, flowdependent estimate of the BEC. Hybrid techniques can be implemented easily with preexisting variational DA systems and can produce comparable forecasts as EnKF but with substantially smaller ensemble members Wang et al., 2007). Many previous studies using hybrid DA system (Buehner, 2005;Buehner et al., 2010aBuehner et al., , 2010bWang et al., 2008b;Kutty et al., 2018, Wang, 2011, Zhang & Zhang, 2012Wang et al., 2013;Zhang, Huang, & Poterjoy, 2013) have shown that hybrid technique produces better results compared to those obtained from 3DVAR and 4DVAR DA schemes. Although these studies (Lu et al., 2016. andLu et al., 2017) show the positive impact of the flow-dependent BEC estimated from different hybrid DA techniques on TC analysis and forecast compared to static BEC estimated from nonhybrid DA, there are limited studies over North Indian Ocean (NIO).
NIOs have different characteristics from other cyclone formation basins due to their unique topography and freshwater discharges from several rivers. The TCs form over NIO, mainly in premonsoon and postmonsoon seasons. The climatological frequency of occurrence of TCs over NIO is 4. The number TCs form during the postmonsoon season is more than that during premonsoon season over NIO. During the period 1990-2018, a total of 15 numbers of TCs underwent one or two periods of rapid intensification (RI) episode(s) in their lifetime over NIO. Among these 15 TCs, five (one over the Arabian Sea [AS] and four over the Bay of Bengal [BoB]) were originated in premonsoon and nine (four over the AS and five over the BoB) in postmonsoon  (Mahrt & Ek, 1984;Mahrt & Pan, 1984;Pan & Mahrt, 1987;Ek, 2003)

Cumulus parameterization
Kain-Fritsch scheme (Kain, 2004) season and remaining one (TC Gonu, 2006) as an onset vortex. It suggests that the probability of occurrence of TCs which undergo rapid changes in intensity are more in the postmonsoon season than premonsoon season. The three postmonsoonal TCs, namely, Phailin (2013), Nilofar (2014), and Chapala (2015), originated over NIO and, underwent RI, rapid weakening (RW), and the eyewall replacement cycle (ERC) during their life cycle have been selected for this study to investigate the simulation of their structure and rapid intensity changes. Figure 1d shows the track of these TCs. This paper aims to evaluate the performance of different DA techniques for accurate prediction of the (rapid) changes in structure and intensity of TCs and influence of environmental factors on these changes so that developed reanalysis can be used for understanding the physical mechanisms associated with them. Therefore, the three TCs mentioned above have been simulated using the Weather Research and Forecasting (WRF) model and three DA techniques, namely, 4DVAR, 3DEnVAR, and 4DEnVAR. WRF model has been widely used to simulate TCs and develop reanalysis. It is essential to know which parameterization and DA techniques simulate the cyclone features most accurately compared to the observations. The recent study (Lu et al., 2016) uses a different combination of DA systems (hybrid and nonhybrid) in NWP to study the structure and intensity changes of TCs and showed that while one DA scheme simulates the intensity better during RI, the other schemes provide the better track simulation (structure) simulation. Rather than choosing a random DA scheme in the WRF model, we aim to study the best suitable DA system to study TCs structure and intensity changes. The results obtained from the simulations are evaluated to identify the better DA technique in this paper. The paper is organized as follows. The following section describes the data sets used and the methodology adopted for this study. The results obtained from different DA techniques and their validation and interpretations are presented in section 3. The conclusions drawn from the study are given in section 4.

Data and Methodology
The dynamical downscaling generally used for obtaining fine-scale data sets using continuous integration of a regional model with a single initialization or consecutive integrations with frequent reinitialization of the atmospheric fields and periodic updates of the lateral boundary condition. For the downscaling, the initial and boundary conditions are obtained from global reanalysis data sets or global circulation models. In this study, we have used Weather Research and Forecasting (WRF) Model Version 3.9.1 (Skamarock et al., 2008) to develop high-resolution (6km × 6km) reanalysis with the two-way nested computational domains shown in Figures 1a and 1b. WRF model is initialized with the National Centers for Environmental Prediction Global Forecast System (GFS) analysis with a resolution of 0.5°× 0.5°(~55 km × 55 km). The 6-hourly cyclic (00, 06, 12, and 18 UTC) DA approach is used to downscale GFS reanalysis using the WRF model and WRF-DA system. The simulations are initiated 2 days before the genesis of each TCs and carried out until 1 day after the landfall in continuous cyclic assimilation mode. The first initial condition for each simulation is obtained from GFS analysis using WRF Preprocessing System. After improving the initial and updating boundary conditions, 12 hr simulations are carried out using the WRF model. The consecutive initial conditions are obtained from the sixth-hour output from WRF simulations. The first 6 hr of the model run are ignored as a spin-up time. Figure 1c shows the schematic diagram of the experimental setup used for the DA cycle and model simulations. For investigating suitable DA techniques for simulation of changes in intensities and structure of TCs, the initial conditions are modified independently by using three DA techniques, namely, four-dimensional DA (4DVAR), hybrid three-dimensional variational DA and EnKF (3DEnVAR), and hybrid 4DVAR and EnKF (4DEnVAR) with the assimilation window of ±1 hr. The boundary conditions also updated using respective DA techniques used for modification of initial conditions. The cost functions of these DA techniques are given by equations (1)-(5), and the test for correctness is described by equations (6) and (7).
3. H is the observation operator. 4. R is the observation error covariance. 5. B is the background error covariance matrix, which is estimated through the National Meteorological Center method (Parrish & Derber, 1992). In this work, the climatological background error was calculated using 12 and 24 hr forecast at 00 and 12 UTC for the duration of 90 days from 1 October 2016. 6. M is the forecast model, and the whole assimilation window is split into i observation windows. 7. x 0 , the control variable is the initial state of the model at the beginning of the time interval. 8. y i is the observation vector. 9. N is the number of observational vectors distributed over time.
10. x ′ i being the ensemble perturbation for the ith ensemble, 11. α i is the effect of ensemble weight, that is, ensemble perturbation weight, and 12. C is the correlation matrix (localization applied to the ensemble perturbations).
Here β s and β e are the weights assigned for the static and ensemble covariance, respectively, where 1 The weight between static and ensemble background errors was selected based on the recommendation from Schwartz and Liu (2014) and given as 75% to the ensemble and 25% to the static BEC in both hybrid 3DEnVAR and 4DEnVAR DA.

Earth and Space Science
The cost functions of these DA techniques are optimized to obtain the best initial condition for simulation. WRF tangent linear model (TLM) and adjoint model system (WAMS) are discussed by Xiao et al. (2008). WAMS uses adjoint sensitivity analysis and the 4DVAR system (Zhang, Zhang, & Poterjoy, 2013). WAMS uses the tangent linear and adjoint test procedure following Thepaut and Courtier (1991) and Navon et al. (1992). The correctness of TLM and adjoint relation has been tested using equations (6) and (7) where ‖·‖ denotes the norm of the vector, f(x), g f (x, g x ), and a f (x, a x ) are full physics WRF model, TLM, and adjoint model, respectively, and x, g x , and a x are the column vectors of model state variables, perturbation state variables, and adjoint state variables, respectively. WRF hybrid DA system implemented horizontal covariance localization (Wang et al., 2008a(Wang et al., , 2008b. The horizontal covariance localization was modeled by applying the recursive filter on the extended control variables during variation minimization. The ensembles covariance localization scales used for each four-weighing factor is described in Wang et al. (2008aWang et al. ( , 2008b. WRF hybrid DA also uses an inflation factor to ensemble transform Kalman filter ensemble perturbations to ensure that on the BEC estimated from the spread of different ensemble about the ensemble mean consistent with the BEC estimated from the difference between the ensemble mean and the observation (detailed described in Wang et al., 2008aWang et al., , 2008b. Xiao et al. (2008) and Zhang, Huang, and Pan (2013) showed that adjoint and tangent linear relationships are matched well till 24 hr forecast.  (Figures 2i and 2l), respectively. The figure shows that the higher records of satellite radiances are assimilated in the initial hours of cyclic assimilation. However, they showed a decline with the increase in the number of assimilation cycles after variational bias correction. The amount of data getting assimilated in the model depends on the various factors such as the number of observations available at each hour, quality control, first guess and innovation, thinning algorithm, variational bias corrections, etc. In this paper, we have used a thinning radius of 90 km. The changes in innovation with two different cyclic iterations for different channels of satellites are provided in the suporting information. As innovation are different for three different DA techniques due to the deviation of the first guess with respect to observations; therefore, obviously, the number of observations assimilated differs. In general, due to the systematic error of the reference field (background or analysis field) with respect to the observations, and the choice of Radiative Transfer Model (Community Radiative Transfer Model in present work), the satellite radiances are considered to be biased. The expression by which the bias for a given satellite channel is modeled depends on few numbers of parameters. The parameters for all sensors on a particular satellite are estimated from time series of radiance departures (observed-minus-background). Then data assimilated system assimilates these bias parameters, and the analysis is corrected accordingly in real time. Now if the bias of some channels of satellite sensor is less (radiances calculated from background fields) is not much deviating from those observed by the satellites, then the radiance data of that channel gets rejected as there is no innovation. In WRF-DA, the variational bias corrections (Derber & Wu, 1998) are inbuilt and does the bias correction for the next cycle. The same can be seen from the simulated tracks ( Figure 1d) and a number of observations assimilated (Figure 2) that the relatively large number of observations assimilated during initial cycles when there is a relatively large track error. However, the reduction in the number of radiance observations is seen when trajectory matches well with the best track. Thus, the decrease in the number of observations assimilated after few cycles can be attributed to the convergence of the cost function near to the minimum due to assimilation in the initial few cycles and the marginal improvement in the latter assimilation cycles. In other words, the reduction in the radiance data assimilated in the analysis is attributed to the improvement in the initial condition due to cyclic DA. In the cyclic DA, the 6-hourly out of the WRF model of the previous cycle provides background conditions Figure 2. Number of conventional observations and satellite radiance assimilated at each cycle in 3DEnVAR and 4DEnVAR DA after quality control and bias correction. The first column (a-f) shows in situ observations assimilated in the initial conditions using 3DEnVAR and 4DEnAR DA methods for the simulation of TCs Phailin (a, b), Nilofar (c, d), and Chapala (e, f) respectively. The second (g-i) and third (j-l) columns show the satellite radiance observations assimilated in the initial conditions for simulations of TCs Phailin (g, j), Nilofar (h, k), and Chapala (i, l) using 3DEnVAR and 4DEnVAR techniques, respectively.
for the next cycle (at same resolution). It reduces the bias in the background with respect to the observations causing smaller perturbation in the model state variables and inturns the rapid convergence of minimization of the cost function. It can also be seen from Figure 1d that the tracks of three TCs matches well with the IMD best track after two to three cycles. The detailed analysis of the rejection of radiance of different channels of satellites is out of the scope of this paper and will be presented in another paper.
For validation of the model performance and to check the performance of DA techniques, we have compared the simulated track and intensity with the IMD best track data set ( Figure 1d) and observations of reflectivity from Doppler Weather Radar (DWR) of IMD Vishakhapatnam. Three-hourly accumulated rainfall data sets from Tropical Rainfall Measuring Mission (TRMM) 3B42v7 are used to calculate skill scores of the model simulation for TCs Phailin and Nilofar. Daily accumulated rainfall retrieved using Integrated MultisatellitE Retrievals for Global Precipitation Measurement Mission (GPM) algorithm from the GPM is used to calculate skill scores of the model forecast of TC Chapala. Meteosat Second Generation satellite image and image (brightness temperature superimposed on the Meteosat-7 satellite image) developed by the U.S. Naval Research Laboratory are used to study the double eyewall structure of TC Phailin and Chapala. The quantitative verification of model rainfall is carried out using equitable threat score (ETS), bias (BIAS) score, false alarm ratio (FAR), and probability of detection (POD) using and TRMM and GPM rainfall data sets. The formula for these indices and details of the calculations are provided by Jolliffe and Stephenson (2012). For the verification of the simulated structure and intensity of the TCs various environmental parameter and such as equivalent potential temperature (ϴ e ), relative vorticity, relative humidity, minimum sea level pressure, and maximum surface wind speed are analyzed with the aim to find best assimilation technique to simulated the observed structure and intensity change of the TCs.  techniques from IMD best track data set at every 6 hr interval are shown in Figure 3. These differences for the TCs Phailin, Nilofar, and Chapala are shown in Figures 3a-3c, 3d-3f, and 3g-3i, respectively. Figure 3a shows minimum track errors estimated by 4DEnVAR simulations at the initial hours of the lifecycle of TCs. During the RI of the TC Phailin (10 to 11 October), ΔR decreases to its minimum value of~2.1 km in 4DVAR and 3DEnVAR simulations and 0 km in 4DEnVAR simulation (at 06 UTC 11 October). At this period 4DVAR simulated track shows comparatively large values of ΔR, ΔV max , and ΔMSLP. The maximum magnitude of ΔR in all three simulations are around 162 km at 18 UTC 08 October for 4DEnVAR simulations. The 4DEnVAR technique shows more values of ΔR compared to the other two methods at the time of dissipation. On average, 4DVAR and hybrid 3DEnVAR estimates comparable ΔR for Phailin. Figure 3b shows the overall values of ΔV max are minimum in 4DEnVAR simulation and maximum in 4DVAR simulations; however, during initial intensification, 3DEnVAR simulation shows minimum ΔV max . The magnitudes of ΔV max show decrease during the RI phase of TC Phailin (10-11 October), and 4DEnVAR captures the intensity of TC more appropriately compared to others. The minimum value of ΔV max is at 12 UTC, 10 October. 4DVAR shows a minimum (~1.59 kt), and 3DEnVAR shows the maximum (~8.2 kt) value. Figure 3c indicates that the ΔMSLP value is maximum (~15 hPa) at 00 UTC of 10 October, and the error is less in 4DVAR and large in 4DEnVAR simulation compared to other two simulations. The minimum ΔMSLP value (~0.6 hPa) is at 18 UTC of 11 October in 4DEnVAR simulations. The overall MSLP prediction is better in 4DVAR during the period of initial intensification, but with the intensification, MSLP prediction improves in hybrid assimilations, and 4DEnVAR simulation shows less ΔMSLP in case of TC Phailin. October. During the same time, the values of ΔV max simulated by both 4DVAR and 4DEnVAR are~31 kt, whereas that by 3DEnVAR is~27 kt. During the initial period of the TC Chapala (Figure 3i), the ΔMSLP is more in 4DEnVAR (peak value~20 hPa at 12 UTC) of 12 October compared to 3DEnVAR and 4DVAR, which estimated these values as~4 and~10 hPa, respectively. Overall, ΔMSLP value is lesser in 3DEnVAR simulations compared to the other two simulations. During the most intensified stages of Chapala, the estimated ΔMSLP values are smaller in 4DEnVAR compared to the other two simulations. From Figures 3a, 3d, and 3g, it is seen that the errors in the TCs track is more before its intensification to CS stage and reduce with the intensification to the higher-intensity stages. However, at the dissipation stage of the TCs, at depression (D) and deep depression (DD) stage, these errors are considerably less compared to these stages during the period of initial intensification. The different peaks in track, intensity, and MSLP errors (especially sudden increase in errors in TC Chapala during its RI in 4DEnVAR simulation) shown in Figure 3 from different DA initialized.

Reflectivity
The simulated radar reflectivity (dBZ) has been compared with DWR images of IMD Vishakapattanam just before landfall by the TC Phailin. Figures 4a-4d show the DWR image, 4DVAR, 3DEnVAR, and 4DEnVAR simulated dBZ, respectively, at 14 UTC of 12 October. It is seen that the position of landfall in all three simulations matched very well with the observation with a center position near latitude 19°N and longitude near 85°E. However, the values of dBZ are slightly overestimated in all three simulations. The maximum values of observed dBZ are in between 48 and 52 dBZ and that in the model simulation is in between 50 and 55 dBZ at the inner core of the TC Phailin. The spatial spread in dBZ is maximum in 4DVAR and minimum in the 4DEnVAR simulation.

Forecast Skill Scores
Forecast validation has been performed by calculating forecast skill scores such as ETS, BIAS, FAR, and POD for rainfall forecast to identify the best forecast among 4DVAR, 3DEnVAR, and 4DEnVAR for all three TCs. Skill scores were calculated using 24 hr accumulated rainfall during the RI periods of the TCs. First (a-d), second (e-h), and third column (i-l) of Figure 5 show the skill scores for the TCs Phailin, Nilofar, and Chapala, respectively. The threshold rainfall is indicated along with the x axis of the plots. In the case of Phailin (Figure 5a), ETS in all simulations increases with the increase of threshold rainfall, but after 90 mm, ETS decreases with the increase in threshold rainfall.
For lower threshold rainfall (till 60 mm), 4DVAR ETS is less, and the other two simulations show comparable ETS. For lower to moderate rainfall (60-80 mm), 4DVAR shows maximum ETS, and after 110 mm threshold, 3DEnVAR (4DEnVAR) shows the maximum (minimum) ETS. Therefore, 4DVAR forecasts well the low rainfall events and 4DEnVAR forecasts well the heavy rainfall events compared to other simulations. In the case of Nilofar (Figure 5e), ETS in all simulations decreases with the threshold rainfall, and 3DEnVAR shows the minimum ETS value, whereas 4DVAR and 4DEnVAR show almost equal ETS. In the case of Chapala (Figure 5i), 4DVAR, 3DEnVAR, and 4DEnVAR show comparable ETS for all low to high threshold rainfall, whereas for lower threshold values, ETS in 3DEnVAR and moderate threshold ETS in 4DEnVAR simulation is slightly more than that in the other two simulations. In the case of TC Phailin, BIAS ( Figure 5b) in all three simulations increase with the increase in threshold rainfall. For lower to the moderate threshold, BIAS is slightly less in 4DEnVAR than other two simulations, whereas for moderate to high threshold values of rainfall, BIAS in the 3DEnVAR simulation is less than the other two simulations. BIAS in the 4DVAR simulation is maximum for all threshold values. Similar to Phailin, BIAS in 4DVAR, 3DEnVAR, and 4DEnVAR simulations of TC Nilofar (Figure 5f) increases with the increasing threshold values of rainfall. For the almost entire range of threshold rainfall, the 4DEnVAR simulation shows the lowest BIAS, whereas 4DVAR and 3DEnVAR show approximately equal BIAS till 90 mm threshold. After this threshold, BIAS in the 4DVAR simulation is maximum. But in the case of Chapala (Figure 5j), BIAS is more in the 3DEnVAR simulation and less in 4DEnVAR simulation than other two simulations for almost the entire range of threshold rainfall. For large values of precipitation, BIAS is maximum in 4DEnVAR and minimum in 4DVAR. The BIAS increases with the rise in threshold rainfall values. In TC Phailin (Figure 5c), FAR is more in 3DEnVAR than 4DVAR and 4DEnVAR simulations till threshold rainfall equal to 70 mm. For moderate to heavy rain in Phailin, FAR is less in 3DEnVAR, whereas it is more in the 4DVAR simulation. For very high threshold rainfall, 3DEnAR and 4DEnVAR simulations show almost equal FAR. In TC Nilofar (Figure 5g), for the entire range of threshold rainfall, the 3DEnVAR simulation shows the maximum FAR, and 4DEnVAR simulation shows the minimum value of FAR. In the case of Chapala (Figure 5k), for a lower threshold of the rainfall (till 40 mm), the FAR value is maximum in 4DVAR and minimum in 4DEnVAR simulation, and for moderate thresholds (50-100 mm) maximum FAR is in the 3DEnVAR simulation. The minimum FAR is simulated in 4DEnVAR simulation till 130 mm threshold rainfall. After 130 mm threshold rainfall, the minimum FAR is in the 4DVAR simulation. The results from the simulation of TCs Phailin and Nilofar (Figures 5d and 5h) show the minimum POD in 3DEnVAR simulations, whereas that in the case of TC Chapala (Figure 5i) is found in 4DEnVAR. The 3DEnVAR simulations show minimum (maximum) POD for almost the entire range of threshold rainfall other than very high rain (above 130 mm) for TCs Phailin and Nilofar.

Intensification of TCs
The evolution intensity of TCs in their entire lifespan can be determined in terms of maximum surface wind (V max ) and MSLP ( Figure 6). Figures 6a-6c show V max and MSLP simulated in 4DVAR, 3DEnVAR, and 4DEnVAR for the entire life span of the TCs Phailin, Nilofar, and Chapala respectively. The second column of Figures 6d-6f shows the boundary layer and midtropospheric relative humidity in the core region (2°× 2°g rid box around the center) for the TCs Phailin (Figure 6d), Nilofar (Figure 6e), and Chapala (Figure 6f). Figures 6a-6c show the intensification as an increase (decrease) in V max (MSLP) with time and dissipation as an increase (decrease) in MSLP (Vmax). All three TCs underwent RI. According to IMD best track data set, the V max for the TC Phailin increased from 45 kt from 00 UTC 10 October to 110 kt on 00 UTC of 11 October, which implies the RI of the cyclone. During this period, 4DVAR, 3DEnVAR, and 4DEnVAR simulations show the increase in V max by 15, 20, and 22 kt, respectively, and~40 hPa decrease in the MSLP in 4DVAR and 3DEnVAR simulations whereas about 43 hPa decreased in 4DEnVAR. In general, RI of TC Phailin is underpredicted in all three simulations, whereas during 14 UTC 10 October to 14 UTC 11 October, 3DEnVAR and 4DEnVAR simulations captured almost~30 kt increase in V max in 24 hr, which implies the RI of Phailin during this period. Similarly, the observed RI of Nilofar (Figure 6b) from 00 UTC 26 October to 00 UTC of 27 October is underpredicted in all three simulations. However, 4DEnVAR estimates captured RI early and showed an increase in V max during 23 UTC, 25 October to 23 UTC, 26 October, by approximately 33 kt. During the same period, the 4DVAR and 3DEnVAR simulations showed an increase in V max by~27 and~28 kt, whereas both the simulation captured RI early by 2 hr with an increase in the value of V max by~31 kt. IMD reported the RI in cyclone Chapala (Figure 6c) during 00 UTC, 29 October to 00 UTC, 30 October, which is also underpredicted in 4DEnVAR simulations, whereas 3DEnVAR and 4DVAR simulations show an increase in values of V max less than 35 kt during this period. 4DEnVAR shows delayed RI of Chapala during 18 UTC, 29 October to 18 UTC, 30 October, with an increase of 30 kt in V max in 24 hr. The 4DEnVAR overpredicts the intensity of the TC Chapala. Due to the overprediction of TC intensity by 4DEnVAR technique could not capture the rapid changes in the intensity of TC Chapala.
Although the overall intensification and dissipations patterns and temporal evolution of TCs are well captured in model simulations, the intensity and MSLP of TCs are overpredicted during the initial period and underpredicted during the more intense stages of TCs in all three simulations. Therefore, model simulations can capture the rapid increase in intensity, however, unable to reproduce it at the appropriate time. The intensity and MSLP predictions are better during the dissipation of TCs in all three simulations. All simulations appropriately showed the increase in the values of relative humidity in the boundary layer and the middle troposphere with intensification and decreased with the dissipation of TCs.

Structure of the TCs
This section is aimed to investigate which simulation captures the spatial and vertical structure of TCs more appropriately. The parameter equivalent potential temperature (ϴ e ) follows the path of moist air in the free atmosphere with a constant value. In a mature TC, the ϴ e surfaces are different in the inner core of the cyclone than the outer radii. In the inner region, values of ϴ e decrease with height mostly up to 600 hPa pressure level, followed by an increase in the values of ϴ e . Therefore, the ϴ e isopleths are bending radially outward after around 3 km height. For the outer region with a radius larger than 70 km, the ϴ e at surface levels is almost horizontal. In the innermost region, near to the center, ϴ e increases radially inward, and the ϴ e isopleths are weakly slanted outward in the upper level. This upward and outward pattern of ϴ e isopleths is due to the upward and outward flow of the primary eyewall.
At the time of TC genesis, the spatial distribution of MSLP and 10 m wind are analyzed for all three TCs and shown in Figure 7. The vertical cross-sectional structure of TCs in terms of relative vorticity and ϴ e surfaces are shown in Figure 8. Spatial distribution of simulated sea level pressure (SLP) (contour) and 10 m wind (shaded) in the core region of TC Phailin at 12 UTC of 9 October from 4DVAR, 3DEnVAR, and 4DEnVAR simulations are shown in Figures 7a-7c (first column). That for TC Nilofar at 12 UTC of 26 October (when IMD declared Nilofar as a severe cyclonic storm) are shown in Figures 7d-7f (second column) and for TC Chapala at 00 UTC of 29 October are shown in (third column) Figures 7g-7i. Figure 8 is similar to Figure 7 but shows the vertical cross-sectional structure of the TCs, where contour lines represent the ϴ e surfaces, and shaded regions represent the relative vorticity (ξ × 10 −5 ) at the core of the TCs. From Figures 7a-7c, it can be seen that the spatial distribution of surface wind (10 m wind) speed is almost similar in 4DVAR, 3DEnVAR, and 4DEnVAR simulations of TC Phailin at cyclonic storm (CS) stage. The maximum 10m wind is in range of 50-60 kt in all three simulations, which are overprediction of the CS and severe cyclonic storm stages of TC Phailin.
Also, the SLP distributions are similar in all three simulations with a minimum pressure of 994 hPa. The vertical structure of Phailin (Figures 8a-8c) is more organized in 4DEnVAR than 3DEnVAR and 4DVAR simulations. Vertically upward and outward propagated ϴ e surfaces show the existence of more prominent Eyewall near the center of the cyclone in 3DEnVAR and 4DEnVAR simulations. In 4DVAR simulated Phailin, maximum ξ (200-250) is in the boundary layer (900-700 hPa) far away from the center of the storm, whereas 3DEnVAR and 4EnVAR show the maximum ξ at the center of the storm. At the middle-upper tropospheric layer of Phailin, maximum ξ is captured in 3DEnVAR simulations.
During the Severe Cyclonic Storm stage of TC Nilofar, all three simulations captured the almost similar pattern of the spatial distribution of SLP and 10 m wind speed (Figures 7d-7f), whereas the vertical structure (Figures 8d-8f) show that the most intense vertical structure is captured in the 4DEnVAR simulation. The maximum value of ξ is~400 × 10 −5 s −1 simulated near to the center of TC Nilofar by 4DEnVAR technique, however, that using 4DVAR and 3DEnVAR techniques is simulated in the eyewall region. The radially inward-upward-outward ϴ e surfaces show the structure of TC Nilofar and existence of eyewall near to the center. The spatial distribution of 10 m wind speed and SLP of TC Chapala in the CS stage (Figures 7g-7i) show that the most intense structure has been simulated in 4DEnVAR with values of MSLP and V max as Figure 9. Spatial distribution of model-simulated (a-c and h-j) radar reflectivity (dBZ) and (e-g and l-n) relative humidity (in percent) at 850 hPa using 4DVAR (first column), 3DEnVAR (second column), and 4DEnVAR ( (Figures 8g-8i) show that the maximum ξ (~400 × 10 −5 s −1 ) is at the center of 4DVAR simulation, and the minimum is in the 3DEnVAR simulation. The ϴ e surfaces show the existence of eyewall most prominently (with larger eyewall compared to other simulations) in 4DEnVAR simulation.

Double Eyewall Formation
The TCs with concentric eyewall experience significant changes in the intensity. Once concentric eyewalls form, the secondary eyewall contracts, and the primary eyewall gets replaced by the outer secondary eyewall. The previous observations showed that most of the intense TC undergoes one or multiple ERCs in their total lifespan. The mechanism of ERC consists of three parts, and these are the formation of the secondary eyewall at the outer radii, the decay of the primary eyewall, and the establishment of the large eye with a broader eyewall and a quasi-axisymmetric structure. During the contraction and intensification of the secondary eyewall, the TC stops intensifying and starts to weaken. After the replacement of the primary eyewall by the secondary eyewall, the TC can resume its intensification if the environmental condition is favorable. Therefore, to predict the intensity of TCs, it is essential to study the ERC to gain a proper understanding of the mechanism involved in the process. The TCs Phailin and Chapala underwent two and one ERC in their entire lifespan, respectively. Figures 9d and 9k show brightness temperature superimposed on Meteosat-7 Visible satellite image for TC Phailin at 00 UTC 12 October 2013 and for TC Chapala at 00 UTC of 11 November 2015, respectively. These images show the existence of the concentric double eyewall structure of the TCs at that time. The model simulated radar reflectivity (Figures 9a-9c) and 850 hPa relative humidity (Figures 9e-9g) at 00 UTC 12 October 2013 show two maxima near the center of the Phailin in all three simulations. Therefore, the simulations are successful in simulating the concentric eyewall structure and ERC of the TC Phailin ( Figure 9d). Also, the model simulated radar reflectivity (Figures 9h-9j) and 850 hPa relative humidity (Figures 9l-9n) at 00 UTC 12 October 2015 show two maxima near the center of the Phailin in all three simulations indicating the primary and secondary eyewall of the TC Chapala. Therefore, the simulations are successful in simulating the concentric eyewall structure and ERC of the TC Chapala (Figure 9k).

Conclusions
In this paper, we have compared the three DA techniques, namely, 4DVAR, 3DEnVAR, and 4DEnVAR methods to generate high-resolution reanalysis for the three TCs occurred over NIO, namely, Phailin (2013), Nilofar (2014), and Chapala (2015). The mesoscale model WRF used for the simulations is initialized with GFS analysis and carried out in cyclic mode to reduce the errors in the developed analysis. At each cycle to the initial conditions were improved by assimilating the GDAS radiance data set is of the weather observations. The developed reanalysis is compared with the IMD best track data set for validation of track, maximum surface winds (intensity), and mean sea level pressure (MSLP) prediction. Also, the simulated double eyewall structure of TC Phailin and Chapala is validated with brightness temperature superimposed on the Meteosat-7 image (produced by NRL). The rainfall forecast along the track of these TCs s are verified with TRMM and GPM rainfall by calculating forecast skill scores, namely, ETS, BIAS, FAR, and POD. It is found that hybrid DA techniques perform better than the nonhybrid 4DVAR method. Recently, Hamill, Whitaker, Kleist, et al. (2011), Kutty et al. (2018), Lu et al. (2016Lu et al. ( , 2017, Schwartz et al. (2013), and Wang (2011) used hybrid DA techniques to study TC and found that the track and intensity prediction of TC initialized with hybrid DA was significantly better than the track forecast generated using only variational method. The result from the present study also matches well with these studies.
The generated reanalysis data set reasonably matches with the best track data set provided by IMD. The errors in the predicted track reduce with the intensification of TCs; however, the error in the simulated maximum surface wind and MSLP decrease with the intensification. The intensities of the TCs are overestimated TCs at the initial intensification stages in all three simulations. In general, 4DEnVAR techniques provide better estimates with comparatively minimum errors. The simulated position and pattern of the spiral band of TC Phailin shown in radar reflectivity match well with the IMD observed radar reflectivity image with an error of 10 dBZ in its value. The simulations using hybrid DA techniques captured the observed rapid intensity changes in TCs, however, not able to predict the exact timing of the RI of the TCs. The temporal error in RI prediction is less than ±6 hr. The study shows that the 4DEnVAR simulations estimated the most intense and organized vertical structure of TCs. The double eyewall structure of TC Chapala and Phailin are simulated appropriately by all three simulations. In general, it can be concluded that the 4DEnVAR technique is most suitable for the generation of useful reanalysis for the understanding of the physical process involved in structure evolutions and rapid intensity changes of the TCs over NIO.