Sensitivity of simulating Typhoon Haiyan (2013) using WRF: the role of cumulus convection, surface flux parameterizations, spectral nudging, and initial and boundary conditions

. Typhoon (TY) Haiyan was one of the most intense and highly destructive tropical cyclones (TCs) to affect the Philippines. As such, it is regarded as a baseline for extreme TC hazards. Improving the simulation of such TCs will not only improve the forecasting of intense TCs but will also be essential in understanding the potential sensitivity of future intense TCs with climate change. In this study, we investigate the effects of model conﬁguration in simulating TY Haiyan using the Weather Research Forecasting (WRF) Model. Sensitivity experiments were conducted by systematically altering the choice of cumulus schemes, surface ﬂux options, and spectral nudging. In addition to using the European Centre for Medium-Range Weather Forecasts Reanalysis ﬁfth-generation (ERA5) single high-resolution realization as initial and boundary conditions, we also used 4 of the 10 lower-resolution ERA5 data assimilation system (EDA) ensemble members as initial and boundary conditions. Results indicate a high level of sensitivity to cumulus schemes, with a trade-off between using Kain–Fritsch and Tiedtke schemes that have not been mentioned in past studies of TCs in the Philip-pines. The Tiedtke scheme simulates the track better (with a lower mean direct positional error, DPE, of 33 km), while the Kain–Fritsch scheme produces stronger intensities (by 15 hPa minimum sea level pressure). Spectral nudging also resulted in a reduction in the mean DPE by 20 km, and varying the surface ﬂux options resulted in the improvement of the simulated maximum sustained winds by up to 10 m s − 1 . Simulations using the EDA members initial and boundary conditions revealed low sensitivity to the initial and boundary conditions, having less spread than the simulations using different parameterization schemes. We highlight the advantage of using an ensemble of cumulus parameterizations to take into account the uncertainty in the track and intensity of simulating intense tropical cyclones

of cumulus schemes, surface flux options, and spectral nudging. In addition to using the European Centre for Medium-Range Weather Forecasts Re-analysis 5th Generation (ERA5) single high resolution realization as initial and boundary conditions, we also used four of the ten lower resolution ERA5 Data Assimilation System (EDA) ensemble members as initial and boundary conditions. Results indicate a high level of sensitivity to cumulus schemes, with a trade-off between using  Fritsch and Tiedtke schemes that have not been mentioned in past studies of TCs in the Philippines. The Tiedtke scheme simulates the track better (with a lower mean Direct Positional Error (DPE) of 33 km), while the Kain-Fritsch scheme produces stronger intensities (by 15 hPa minimum sea level pressure). Spectral nudging also resulted in a reduction in the mean DPE by 20 km and varying the surface flux options resulted in the improvement of the simulated maximum sustained winds by up to 10 ms -1 . Simulations using the EDA members initial and boundary conditions revealed low sensitivity to the initial and 25 boundary conditions, having less spread than the simulations using different parameterization schemes. We highlight the advantage of using an ensemble of cumulus parameterizations to take into account the uncertainty in the track and intensity of simulating intense tropical cyclones.

Introduction
As a country of 109 million people over more than 7,000 islands, the Philippines is considered one of the most natural hazard-30 prone countries in the world (Brucal et al., 2020) and is ranked in the top five of all countries in terms of exposure to climaterelated risks (Eckstein et al., 2020). One of the most important hazards the Philippines is exposed to is tropical cyclones (TCs).
TCs bring intense winds, extreme precipitation, and storm surges that affect a large portion of the Philippine population (Bagtasa, 2017;Lyon and Camargo, 2009). Due to its location in the western North Pacific Ocean, where TC formation is conducive all year, the Philippines is exposed to an average of 10 landfalling TCs annually (Cinco et al., 2016). Since 1990, 35 TCs in the Philippines have resulted in up to half of the total losses from all natural disasters amounting to about USD 20 billion in damages (Brucal et al., 2020) and an annual average death toll of 885 with estimated accumulated deaths due to TCs of approximately 30,000 from 1980 to 2013 (Yonson et al. 2016). It is estimated that about 5 million people are affected annually or over 570,000 are affected on average per destructive TC (Brucal et al., 2020). 40 https://doi.org/10.5194/nhess-2021-400 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.
One of the strongest typhoons that made landfall in the Philippines in recent history is Typhoon (TY) Haiyan (locally named "Yolanda"), which is considered the second costliest Philippine TC since 1990 (EM-DAT, 2020) and one of the deadliest since the 1970s (Cinco et al., 2016;Lander et al., 2014;Lagmay et al., 2015). TY Haiyan was a category-5 super typhoon that claimed the lives of at least 7,300 people, most of them from drowning due to the devastating 5 to 7-meter-high storm surge and coastal inundation (Soria et al., 2016). It also affected more than 16 million people (NDRRMC, 2014) and caused an 45 estimated USD 5-15 billion worth of damage particularly in agriculture and critical infrastructure (Brucal et al., 2020). Comiso et al. (2015) found that TY Haiyan coincided with the warmest sea surface temperature (SST) observed over the Pacific warm pool region which may have contributed to its intense nature. This relation between intense TCs and warmer tropical SSTs, has also been found in the Atlantic (Emanuel 2005), and suggests that continuous warming may lead to more intense TCs in the future. Consistently, an increasing trend in intense TC frequency affecting the Philippines since the 1970's has been 50 observed (Cinco et al. 2016;Comiso et al. 2015). TC rainfall is also expected to increase in the future as TCs intensify (Intergovernmental Panel on Climate Change, 2021; Patricola and Wehner 2018), potentially increasing the risk of flooding and landslides. Given TY Haiyan's intensity and impacts, it is regarded as a benchmark for an intense and destructive TC.
Hence, it is important to test how well it can be simulated in current models in the present climate and what are the sensitivities to model formulation. Improving the representation of intense TCs like Haiyan in regional climate models (RCMs) such as 55 the Weather Research and Forecasting (WRF) model (Skamarock et al. 2008) is essential for simulations of such TCs in different future climate change scenarios to provide credible impact assessments.
While global climate models (GCMs) are very useful for looking at the changes in TC activity under different climate change scenarios (e.g. frequency, intensity, genesis from a climatological and global/regional perspective) (Gallo et al., 2019;60 Villafuerte et al.;2021, Patricola andWehner 2018) and some advances have been made in the past few decades in the use of global convection-permitting models (Judt et al., 2021), previous studies still demonstrate the need for (convectionresolving/convection-permitting) RCM simulations to better simulate the processes relevant to the TC formation and development as well as their properties, particularly the most intense ones (e.g. Walsh et al., 2013). In consideration of the computational cost in resolving important TC processes, the use of RCMs is a valuable and complementary approach to using 65 GCMs in investigating the potential changes in TCs in the future. One such RCM is the WRF Model (Skamarock et al. 2008) which is currently used for operational forecasting by the country's meteorological office -Philippine Atmospheric, Geophysical and Astronomical Services Administration (PAGASA) (Flores 2019;Aragon and Pura 2016) and also used in studies projecting future changes in temperature and rainfall over the Philippines Daron et al., 2018;Tolentino and Bagtasa, 2021). 70 A consensus of past modelling studies of western North Pacific TCs, including TY Haiyan, show the cumulus convection scheme as having the most influence on its intensity over other model parameters such as the planetary boundary layer (PBL) and/or microphysics schemes (Islam et al. 2014;Di et al. 2019). In particular, the Kain-Fritsch (KF) (Kain, 2004) cumulus convection scheme has been found to produce the best TC tracks and wind intensity estimates (Shen et al., 2019;Zhang et al., 75 2011;Spencer and Shaw, 2012;Prater and Evans, 2002;Mohandas and Ashrit, 2012). Furthermore, the often-selected KF scheme was shown to be also sensitive to model resolution (Li et.al., 2018). However, the use of the KF scheme has also shown certain limitations. A study by Torn and Davis (2012) found that the KF scheme produces larger TC track biases than the Tiedtke (TK) cumulus convection scheme. Other than the said parameterization schemes, improvements in simulations of TC intensity are also found to be influenced by the surface flux options (Kueh et al., 2019). 80 In this paper, we revisit the simulation of TY Haiyan using downscaling and assess its sensitivity to model formulation and the driving initial and boundary conditions in preparation for pseudo-global warming and CMIP6 climate projection https://doi.org/10.5194/nhess-2021-400 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License. experiment studies. This study builds on the work of Islam et al. (2015) who assessed the effects of different combinations of Planetary Boundary Layer, microphysics and cumulus convection scheme using WRF but found substantial underestimation 85 of TY Haiyan's intensity regardless of the sensitivity to physics parameterization; Li et al. (2018) who used WRF to look at the effects of the cumulus parameterization at different resolutions (9-2 km) and found that the most effective resolution to simulate TY Haiyan with no cumulus parameterization or a revised KF scheme is at 2-km and 4-km resolution, respectively; and that of Kueh et al. (2019) who looked at the influence of the different surface flux options in simulating TY Haiyan's intensity using one cumulus convection scheme and found that a better representation of surface flux formulas improved the 90 simulated intensity in WRF. Here, we investigate the effects of the different combinations of model cumulus convection schemes, spectral nudging and surface flux options on the TY Haiyan track, intensity and rainfall simulations. From this study the best combination is determined which will then be used for investigating the effects of future climate change on TY Haiyan and other TC cases. The associated storm surge of TY Haiyan (Mori et al. (2014), Nakamura et al. (2017) and Takayabu et al. (2015) is not considered here. Model parameterization scheme sensitivity studies that assess the simulation of TCs provide 95 guidance to future TC modelling studies (Villafuerte et al., 2021). This study also seeks to contribute to such studies with a particular focus on the Philippines by assessing the skill and sensitivity of a TC case study using a mesoscale numerical weather prediction model. In particular, it aims to study the influence of the combination of cumulus convection scheme and the different surface flux options to the different TC characteristics. This study adds on existing literature by looking at the effects of cumulus convection schemes combined with different flux options and spectral nudging. Specifically, it aims to address the 100 following questions:  How sensitive are the TY Haiyan simulations to convective schemes, surface flux options and spectral nudging?
 How sensitive are the simulated track and intensity of TY Haiyan to the uncertainty in the initial and boundary conditions?

105
The results will provide valuable information for regional climate downscaling of intense TCs, which can be used in evaluating the sensitivity of future TCs in climate change simulations. Section 2 provides a description of the methodology. Then the paper continues with the results of the sensitivity experiments followed by the discussion, and finally, Section 4 provides a summary of the findings and recommendations for future work.  (Comiso et al 2015). It then rapidly intensified into a TY 115 on November 5 at 142.9°E, 6.9°N and was classified as a category-5 equivalent super typhoon by the Joint Tropical Warning Center (JTWC) and was classified as a Typhoon by PAGASA, its highest classification at the time. It further intensified before making landfall on November 7 2040 UTC. It traversed the central section of the Philippines and started to slowly weaken to a tropical depression on November 11 (JMA, 2013). Typhoon Haiyan claimed the lives of more than 7,300 people, mostly due to the associated storm surge and coastal inundation. It is estimated to have caused between USD 5-15 billion worth of direct 120 damages in agriculture and infrastructure (Brucal et al., 2020) and affected more than 16 million people (NDRRMC 2014). https://doi.org/10.5194/nhess-2021-400 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.

Model Description
Simulations were conducted using WRF version 3.8.1 (Skamarock et al. 2008), a non-hydrostatic numerical weather prediction RCM developed by the National Center for Atmospheric Research (NCAR). It is used for atmospheric research and operational 125 forecasting, and increasingly for regional climate research (Powers et al., 2016). The model includes a variety of physical parameterization schemes, including cumulus convection, microphysics, radiative transfer, planetary boundary layer, and land surface. The WRF Advanced Research WRF (WRF-ARW) solver uses the Arakawa-C grid as the computational grid and the Runge-Kutta 3 rd order time integration schemes (MMML-NCAR, 2019). Skamarock et al. (2008) provides a more detailed description of the model specifications. PAGASA uses WRF for its operational forecasting over the Philippine Area of 130 Responsibility (Flores 2019;Aragon and Pura 2016) and it is also used in studies projecting future changes in temperature and rainfall over the Philippines .
The land surface information comes from the 30 arc s (*1 km) resolution Moderate Resolution Imaging Spectroradiometer (MODIS) satellite dataset with 20 global land use categories. 135

Initial and boundary conditions
The European Centre for Medium-Range Weather Forecasts (ECMWF) Re-analysis 5 th Generation (ERA5) data is used for both the initial and boundary conditions. It is the latest generation of reanalysis products produced by ECMWF with horizontal resolution of 31 km, hourly temporal resolution and 137 vertical levels (Herbatsch et al., 2020). ERA5 uses observations 140 collected from satellites and in-situ stations, which are quality controlled and assimilated using 4D-VAR, a model based on the ECMWF's Integrated Forecast System (IFS) Cycle 41r2.
Alongside the release of the ERA5 single realisation deterministic data from 1979 to the present, data from the Ensemble of Data Assimilations (EDA) system was also made available. The EDA system is a 10-member ensemble at a lower resolution 145 than the deterministic data (60 km horizontal resolution and 3-hourly) (Hennermann, 2018). The EDA system provides estimates of analysis and short-range forecasts through one control and nine perturbed members, which provide background error estimates for the deterministic forecasts. This system allows for estimating uncertainty since it provides estimates of the analysis and short-range forecast. These are provided as an uncertainty measure, albeit with half the resolution of the reanalysis.
To test the sensitivity to varying boundary conditions, simulations were also performed using four randomly selected 150 representatives of the 10-member ERA5 EDA system. The selected ensemble members were used to test the sensitivity to different perturbed observations, sea surface temperature fields and model physics the (Isaksen et al., 2010).

Design of sensitivity experiments and analysis
In this study, the WRF-ARW model has been configured with two nested domains centered over the point of 18.3° latitude, 155 135° longitude. The outermost grid has 294 x 159 grid points with 25-km grid spacing, while the innermost domain has 745 x 550 grid points with 5-km grid spacing, with 28 vertical levels and the model top pressure level was set to 50 hPa. This model resolution was chosen in favour of using supercomputing resources for the systematic testing of different paramaterization schemes, and in consideration of additional simulations under future climate conditions. Different domain configurations were tested prior to selecting this particular configuration, with the current domain configuration having the track and intensity 160 closest to that observed (not shown). Additional simulations using convection-permitting resolution (single domain, 4.5km) were also performed and showed no significant change in simulated intensity from the configuration used here (not shown).
The model domain set-up is shown in Fig. 1.  Convection is mostly simulated in models with resolution coarser than 10-5km through the cumulus parameterization scheme.
WRF's cumulus parameterization scheme simulates the effects of cumulus convection on heat, moisture, and precipitation at the sub-grid scale (Skamarock et al. 2008). The choice of cumulus parameterization schemes has an impact on WRF's ability 180 to simulate the TC track, intensity, and structure (Zhang et al., 2011;Shepherd and Walsh 2016;Parker et al., 2017). Only two schemes were investigated in this study the KF scheme and TK scheme, the difference of which are summarized below in Table 1. These cumulus schemes are used because PAGASA uses KF for its operational forecasting configurations (Flores 2019). It has also often been used for TC simulation studies in the Philippines and has been found, in several studies, to be the best choice for simulating TC track and intensity (e.g. Sun et al., 2019;Li et al., 2018) and rainfall (e.g. . The 185 TK scheme, on the other hand, has been suggested to be the more appropriate cumulus scheme in tropical weather/climate applications of the WRF model (Parker et al., 2017).
Using a mass flux approach with downdraft removal, and utilizing Convective Available Potential Energy (CAPE), KF is a deep and shallow convection sub-grid scheme that includes clouds, rain, ice, and snow detrainment and cloud persistence (Kain 190 2004). Although KF can account for relatively small-scale processes that drive convection, it has inherent limitations in simulating shallow convection over tropical oceans (Parker et al., 2017). On the other hand, the TK scheme assumes that the moisture flux through the cloud base is equivalent to the surface moisture flux, as well as momentum transport, cloud detrainment, and ice detrainment (Tiedtke 1989;Zhang et al., 2011). According to Parker et al. (2017), the TK scheme is more appropriate for simulating intense TCs in tropical oceans. 195 https://doi.org/10.5194/nhess-2021-400 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.  (2016) 200 Experiments were also conducted to examine the sensitivity to the available parameterizations for surface flux options. For TC applications, WRF-ARW provides three different formulations of aerodynamic roughness lengths of the surface momentum and scalar fields as surface flux options (isftcflx = 0,1, and 2) (see Kueh et al. 2019 for a detailed description of the differences between these options). It has been shown that surface fluxes can influence the model's ability to simulate TC intensity and structure (Green and Zhang, 2013;Kueh et al., 2019). For the default flux option (referred to here as sf0), the 205 momentum roughness length is given as in Charnock (1955) the momentum roughness length is given as Charnock's (1955) expression plus a viscous term, following Smith (1988) -Eq. (1): where is the Charnock coefficient and v the kinematic viscosity of dry air, for which a constant value of 1.5×10−5 m2 s−1 is used. A constant value of α=0.0185 is used for sf0. 210 Since the roughness length formulas in sf0 are demonstrably inconsistent with a substantial amount of research (Kueh et al. 2019), two more options were developed (hereinafter referred to as sf1 and sf2) (Kueh et al. 2019). Based on the findings that the drag coefficient (CD) seemed to level off at hurricane force wind speed (e.g., Powell et al., 2003;Donelan et al., 2004), the surface flux option 1 (sf1) was developed and implemented in WRF as a blend of two roughness length formulas (Green 215 and Zhang 2013). The sf1 option was first implemented in version 3.0 of WRF (Kueh et al. 2019). The sf1 and surface flux option 2 (sf2) has the same momentum roughness length, but in sf2, the temperatures and moisture roughness lengths are expressed in accordance with Brutsaert (1975) (MMML-NCAR, 2019).
A set of experiments is conducted to explore the impacts of nudging on the ERA5 large-scale environment by applying spectral 220 nudging (snON). It has been shown that spectral nudging can improve TC track simulations (Guo 2017;Tang et al., 2017)    Surface flux option (isftcflx)* With spectral nudging With spectral nudging The control simulation is the experiment with KF as cumulus scheme, with spectral nudging turned off and surface flux option of sf0 (KFsnOFFsf0). This configuration was also used in the experiments using the different members of EDA, to test the 235 sensitivity to different initializations.
Other parameterization schemes (adapted from Li et al., 2018) in the model that remained the same in all the experiments, includes: the Rapid Radiative Transfer Model (RRTM) scheme (Mlawer et al., 1997) and the Dudhia scheme (Dudhia 1989) for the longwave and shortwave radiation, respectively; surface layer is from the MM5 Monin-Obukov scheme (Monin and 240 Obukhov 1954); WRF Single-moment 6-class Scheme for the cloud microphysics (Hong and Lim 2006); and Yonsei University (YSU) PBL scheme ; and the Unified Noah Land Surface Model (Chen and Dudhia 2001;Tewari et al. 2004) for the land surface processes and structure.

Verification Data 245
To determine the model's skill in simulating TY Haiyan, we used the International Best Track Archive for Climate Stewardship (IBTrACS) which compiles best-track information from various agencies worldwide (Knapp et al., 2010). We compared the

TC Tracking Method
The simulated track and intensity values were obtained every 6 hours using the TRACK algorithm (Hodges et al., 2017) as used in Hodges and Klingaman (2019). TRACK determines TCs as follows: first the vertical average of the relative vorticity 265 at 850-, 700-, and 600-hPa levels is obtained. The field is then spatially filtered using 2D discrete cosine transforms equivalent to T63 spectral resolution and the large-scale background is removed. The tracking is performed by first identifying the vorticity maxima > 5.0 × 10 -6 s -1 . Using a nearest neighbour method, the tracks are then initialized and refined by minimizing a cost function for track smoothness subject to adaptive constraints (Villafuerte et al., 2021). The feature points are determined by first finding the grid point maxima which are then used as starting points for a B-spline interpolation and steepest ascent 270 maximization method, to determine as the off-grid feature points (Hodges 1995 as cited by Hodges and Klingaman, 2019).
The tracking is done for the entire simulation period. Additional variables are added to the track data after the tracking is complete, such as the maximum 10-m winds within a 6 0 geodesic radius and the minimum sea level pressure (MSLP) within a 5 0 radius using the B-splines and minimization method (Hodges and Klingaman, 2019).    observed track compared to simulations utilizing the TK scheme which are closer to the observed track. The minimum DPE obtained from the simulations using the TK scheme is 8km after 18 hours of simulation for the simulation using TKsnOFFsf2. Fig. 3 shows the sensitivity of the tracks to the cumulus parameterization scheme, surface flux options, and to spectral nudging.

Results and Discussion
The results show that these three model settings individually lead to significant reductions in DPE values. The difference 295 between the mean DPE of simulations using KF and TK schemes (p-value: 0.010) were found to be statistically significant at 99% confidence levels using a Student's t test. The simulations using the TK scheme have mean DPE of 47 ± 5 km and KF scheme have mean DPE of 55 ± 7 km (Fig. 3a). Overall, we found the TK scheme to be best in simulating the track of TY Haiyan.

300
Our results show that the tracks are also slightly sensitive to the use of spectral nudging, especially in the latter half of the simulation (Fig. 3b). The evolution of DPE in Fig. 3 shows gradual increases in its value in the first half of the simulation, as the typhoon approaches land (between 48 h-54 h), the DPE then starts to abruptly increase until the end of the simulations.
This suggests that the spectral nudging configuration does not constrain the model strongly. Nevertheless, simulations run with spectral nudging consistently show lower DPE in the 2nd half of the simulation compared to the no nudging experiments. 305 Moreover, the TK runs appear to also have less dependence on nudging (not shown). The mean DPE of the TK simulations with nudging is 68 km while the simulation without nudging is 87 km. This is consistent with previous studies where spectral nudging improves TC tracks in the WNP (Guo 2017;Moon et al., 2018). Overall, the surface flux options did not have a significant effect on the tracks of the simulated TY (Fig. 3c).   pressure and lower maximum wind speeds. A Student's t test indicates that difference between the minimum sea level pressure simulations using KF and TK schemes (p-value: 0.008) is significant at the 99% confidence level. However, the simulations 325 were not able to capture TY Haiyan's rapid intensification phase as in previous studies (Islam et al., 2015;Kueh et al., 2018).  The KF and TK schemes represent shallow convection differently, resulting in different simulated TC intensities (Torn and Davis, 2017). The TK scheme allows both upward transport of moisture across the boundary layer, as well as vertical advection 355 of evaporation from the ocean surface (Parker et al., 2018). Consequently, this reduces the mass flux in deep convection, thereby lowering the rate of TC intensification and resulting in lower simulated intensities. The KF scheme, however, is less likely to reduce the deep convective mass flux that allows for intensification rates to increase. These results are consistent with the differences in the simulated intensities shown in Parker et al. (2017) and Shepherd and Walsh (2016). Parker et al. (2017) found that the KF scheme produces more intense TC systems (lower MSLP values) than the TK scheme for TY Yasi in 360 Australia. Shepherd and Walsh (2016)  The choice of surface flux option (sf0, sf1, sf2) also affects the ability to reproduce both minimum sea level pressure and maximum winds, as shown by the lower RMSE of sf1 (Fig. 6). Simulations with sf1 have generally shown to have the highest 370 correlation coefficients. While both wind speed and MSLP intensity are strongly dependent on the surface flux option, sf1 shows to simulate the highest intensity for TY Haiyan. As in Kueh et al. (2019), the default option (sf0), in which CD does not level off, the simulation of Haiyan has the weakest wind speeds. The sf1 option is expected to have the highest intensity since it has the largest enthalpy and momentum (Ck/CD) ratio at high wind speeds and lowest CD. This gives less friction at high winds, thereby favouring higher intensity (Kueh et al. 2019). 375 https://doi.org/10.5194/nhess-2021-400 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.
Comparing the simulations with the KF and TK schemes shows that the former produces better simulated intensities with lower biases, RMSE and higher correlation coefficients (Fig. 6), consistent with Zhang et al. (2011) andParker et al. (2017), and minimum sea level pressure as with Spencer and Shaw (2012). 380 Table 3: The resulting deviation from landfall location (km, rounded to nearest whole number); translation speed (ms -1 , rounded to two decimal place); and deviation from observed translation speed (ms -1 , rounded to two decimal places); and deviation of the simulated MSLP at landfall (hPa, rounded to the nearest whole number); compared to observations We also considered the wind-pressure relationship of the simulated intensities of all experiments, which according to Green 385 and Zhang (2013) is affected by surface flux options. The scatterplot in Fig. 7 indicates the relationship between the MSLP and maximum wind, based on the different simulations. The IBTraCS data (black square markers) are also included in this plot. Almost all simulations show a decreasing trend of the MSLP and maximum winds as the storm intensifies; however, the intensities are evidently underestimated (MSLP and maximum wind speeds). Based on Manganello et al. (2012), the maximum wind speed is usually underestimated in RCMs when the simulated MSLP is below approximately 980 hPa. It is worth pointing 390 out that of the different simulations, those utilizing the surface flux option 1 (sf1, blue) give the most intense storm by wind speed (Fig. 8). The simulated maximum wind speeds in the simulations using the default surface flux option (sf0, red) only ranges between 35 to 55 ms -1 while the simulations using the other options (sf1, blue and sf2, cyan) are well distributed from ~40 to 73 ms -1 , consistent with the result of Kueh et al. (2019). Most simulations have an underestimated maximum wind speeds for MSLP below 910 hPa, which is consistent with a study using WRF that produced lower wind speed compared to 395 IBTrACS for a given MSLP (Hashimoto et al., 2015). However, the simulations were able to generate considerable intensity for the maximum wind speed for TY Haiyan compared to that of Islam et al. (2015) who used different model physics options i.e. WRF single moment 6-class (WSM6), WRF single moment 3-class (WSM3), new Thompson (THOM), Milbrandt-Yau double moment (MY2) 7-class scheme, and the Goddard GCE (GGCE) schemes. Previous studies using lower resolution generated insufficient wind speeds in the regime higher than 45 ms -1 (Jin et al., 2015) which are primarily attributed to low 400 model resolutions and deficiencies in surface drag representations at high wind conditions (Jin et al., 2015;Shen et al., 2017).

405
In simulating TCs, it is important to get the timing and intensity at landfall right as it gives a good indication of the potential damage along coastal areas (Parker et al. 2017). TY Haiyan made landfall in the eastern-central Philippines (Guiuan, Eastern Samar) on 7 Nov 2013 at 2040 UTC. Table 3 shows that the simulation with the closest landfall time and location occurs for the KFsnoffsf2 simulation. The deviation from the observed landfall point, minimum deviation is 3 km for KfsnOFFsf2 and 410 TKsnONsf2 and maximum deviation is 76 km for KFsnONsf1, is within the average forecast error for tropical cyclones at 24hour lead time in Western North Pacific (Peng et al, 2017). Fig. 8 also shows that the simulated TY is slightly slower (farther from land on 7 Nov 2013 at 00UTC) than observed, with the timing of landfall delayed between approximately 2 to 6 hours in the simulations. Based on data from IBTrACS, Haiyan's 415 translation speed before landfall is approximately 9.48 ms -1 , while the mean translation speed of all the simulations is 9.43 ms -1 as shown in Table 3. Fig. 8 also shows that the extent of the wind field of the simulations using KF scheme are wider than the ones using TK scheme. The KF scheme simulations have a bigger radial extent, for winds speeds larger than 35 ms -1 or 80 miles per hour (mh -1 ), than the simulations using the TK scheme. The wind field extent is also bigger in simulations with sf1 and sf2 than the ones using the default surface flux option (sf0), with sf1 having a wider and more symmetric radial extent of 420 winds greater than 50 ms -1 or 110 mh -1 . TY Haiyan's radius of maximum wind was estimated to be between 25-29 km (Shimada et al, 2018). In addition, the radial extent of winds of approximately 15 ms -1 (30 mh -1 ) is bigger in simulations using KF than simulations using TK scheme, with radius of maximum wind extending up to ~52 km and ~42 km, respectively. The TKsnONsf0 and TKsnOFFsf0 both have radial extent of winds of 15 ms -1 (30 mh -1 ) that are closer to what is estimated using the OSCAT scatterometer data. 425

Simulated Track and Intensity from ERA5 EDA Ensemble Members
The simulated tracks of TY Haiyan, using the four ERA5 EDA members as initial and boundary conditions and configurations that are the same as used for the control simulation, are found to be within the variability of the simulations using the different 430 https://doi.org/10.5194/nhess-2021-400 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License. parameterizations (Fig. 9). The average DPE of the ensemble mean is 86 km compared to the average DPE of the simulations using different parameterizations which is 78 km with a range from 7 km to 250 km throughout the whole simulation period.
There is no significant difference between the mean DPE of the simulations using the different ensemble members with the simulations using the different parameterization schemes (p-value = 0.464).

445
The spread in the mean bias of the simulated intensities (MSLP and maximum winds) using the ensemble members as boundary conditions is similar to or within the spread of the correlation between the experiments with the different parameterization schemes and spectral nudging option (Fig. 10). Judging from the spread of the simulate intensities found in the boundary conditions experiments, the use of different ensemble members has relatively less effect on the simulated intensities as 450 compared to the sensitivity to cumulus and surface flux parameterizations.

460
The simulated rainfall in WRF is represented implicitly to demonstrate the effects of sub-grid scale processes through the cumulus scheme and explicitly through the microphysics scheme. In this study, we used the combination of both implicit and https://doi.org/10.5194/nhess-2021-400 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License. explicit precipitation as the total rainfall. The spatial distribution of rainfall (mm) from 00 UTC 7 November 2013 to 18 UTC 8 November 2013 from the different experiments without spectral nudging are presented in Fig. 11. These results show a discernible difference between the spatial distribution and magnitude of the simulated rainfall, which indicates high sensitivity 465 to the cumulus schemes. The accumulated 6-hourly rainfall was generally larger in magnitude and spatial extent for the simulations using the KF scheme ( Fig. 9 b-d) than those that used TK scheme (Fig. 9 e-g). There is not much difference in the magnitude and distribution of rain among the different surface flux options.
It is also important to note the delay in the rainfall at landfall, primarily due to the relatively slower movement of the simulated 470 TCs. The extent of the distribution of rainfall outside of Haiyan's inner rain bands were also not captured well by the simulations when compared with the satellite-derived GPM rainfall (Fig. 11a). In comparison with the GPM rainfall, the distribution of the simulated high rainfall using the KF scheme show more similar patterns unlike with the TK scheme. The areas of high rainfall appear to be similar in the simulations using different flux options but different in simulations using KF and TK scheme. The simulations using KF scheme also seem to capture the outer rainbands of TY Haiyan, but extending 475 further southeast compared to the GPM rainfall. Previous studies have also indicated the sensitivity of TC-associated rainfall to different physics parameterizations in WRF. Satuya et al. (2019) and Duc et al. (2019) found that KF better predicts rainfall than TK but both generally perform poorly in simulating rainfall, and WRF TC-associated rain are underestimated (Bagtasa 2021).

Environmental Factors
Based on previous similar studies (Parker et al. 2017;Torn and Davis 2012) and as shown in Fig. 12, the KF scheme results in a warm temperature bias (at 700hPa) and a steering flow bias (at 850 hPa). In particular, the TK scheme produces cooler temperatures and the KF scheme simulates up to approximately 1.5 to 2°C warmer temperatures relative to ERA5, while the ones using the TK scheme have a colder bias at 700hPa (Fig. 12), which is consistent with previous studies (Parker et al. 2017 490 andShepherd andWalsh 2016). On the other hand, the KF scheme is likely to simulate the deep convective mass flux which allows for an increase in intensification rates (Zhu and Smith 2002;Emanuel 1989 as cited by Torn and Davis 2012).
where u, v are the zonal and meridional wind components, respectively, at 200 and 850 hPa. The simulated vertical wind shear is weaker along the track of TY Haiyan for both simulations using KF and TK, but the simulation using KF has a bigger area with weaker shear. It is likely that the more homogeneous temperature field in KF resulted in less vertical wind shear, while 505 the simulation using the TK scheme led to a more heterogeneous temperature increasing the vertical shear. A previous study by Floors et al (2011) showed that the temperature differences through the atmospheric profile leads to geostrophic wind shear in WRF simulations. With the weaker vertical shear, the intensity is higher in the simulation using KF than the simulation using TK scheme. Weaker vertical shear has been found to be favourable in maintaining TC development and intensity (Shen et al., 2019). 510 To further investigate the difference in the track between KF and TK simulation runs, we analysed the 500mb geopotential height. The 5800 m geopotential height contour at 500 mb is used to depict the western North Pacific sub-tropical high (WNPSH) (Xue and Fan 2016). With the ridge location at 20°E, the WNPSH extends to the north of South China Sea (Shen et al., 2019). It has been found that the westward extent and location of the subtropical high ridge directly affect TC tracks in 515 the WNP basin (Ho et al., 2004). In the simulation using the KF scheme, the subtropical high is weaker and is substantially in a more northward position compared to the simulation using the TK scheme (Fig. 13), which likely causes the tracks of the https://doi.org/10.5194/nhess-2021-400 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License. simulations using the KF scheme to drift northward, while the simulations using the TK scheme are much closer to the observed. According to Sun et al. (2014), deep convection in mass flux schemes, such as KF, produces large amounts of anvil clouds that warm the upper troposphere and cause latent heating south of the WNPSH that leads to the weakening of the 520 WNPSH and the movement of the TCs northward. Villafuerte et al. (2021) further added that the use of cumulus schemes results in a weaker subtropical high resulting in shifts in northward re-curvature of TC tracks.

525
The TK scheme also produced relatively drier storm environments along the TC path compared to the simulation using KF scheme and as a result, less convection, which translates into weaker intensity (lower wind speeds), whereas simulations using 530 the KF scheme are ~15% higher relative to the simulation using TK. The TK scheme has relatively drier bias with respect to ERA5 along the TC track (Fig. 14). According to Villafuerte et al. (2021), the TK scheme underestimates mid-tropospheric relative humidity, providing a drier environment thereby constraining deep convection and inhibiting TC development.
Furthermore, Shen et al. (2019) demonstrated that drier lower troposphere enhances downdrafts and inhibit convection, resulting in weaker intensities and less rain. When comparing the distribution of mid-tropospheric relative humidity as shown 535 in Fig. 14, KF shows a lower relative humidity along the track of Haiyan, which indicates better vapor to precipitation conversion of the KF scheme. From this distribution of relative humidity, it can be easily assumed that the KF scheme produces more convection and generated significant rainfall associated with the system, as compared to the weaker convective organization (hence less rainfall) of the simulations using the TK scheme.

545
Typhoon Haiyan (2013) was one of the most intense and destructive tropical cyclones ever to hit the Philippines. As climate models project more intense storms, such as Typhoon Haiyan, will occur more frequently in the future due to climate change, it is important to improve their representation in high-resolution models, in order to improve understanding of TCs under climate change and improve confidence in model projections, and more importantly, for risk and impact assessments. The intensity of TY Haiyan proved difficult to simulate using the Weather Research and Forecasting model at 25-km and 5-km 550 domain configurations as with other previous studies. This study was able to assess the sensitivities to different parameterizations in WRF that can be useful in future simulations of TC cases under future climate conditions. Despite the failure to simulate Haiyan's rapid intensification phase, the simulations were still able to capture the tracks and intensity reasonably well. Based on the results, there seems to be a trade-off between utilizing KF and TK cumulus schemes that has not been previously discussed in previous studies of tropical cyclones in the Philippines. 555 The simulated intensity of TY Haiyan is most sensitive to changes in the cumulus scheme and surface flux options; on the other hand, simulated track is most sensitive to cumulus scheme and spectral nudging. However, the TK cumulus scheme produces better track and the KF scheme produces better intensity. There is a statistically significant difference in the simulated tracks and intensities between the use of the two cumulus schemes. The TK scheme simulates the track better, while the KF 560 scheme produces higher intensities, with the KF scheme simulating a mean bias of 16 hPa and 2 ms -1 and the TK scheme with a mean bias of 31 hPa and -6 ms -1 , respectively. The KF scheme has larger DPEs (mean DPE of 94 km compared to mean DPE of 61 km for TK scheme) due to a more northward steering flow. On the other hand, simulations using the TK scheme had weaker wind and higher MSLP due to the suppression of deep convection by active shallow convection. Simulated rainfall is also sensitive to the cumulus schemes, with simulations using TK having less and smaller rainfall extent than simulations 565 using KF cumulus convection scheme.
The results also show the simulated tracks are sensitive to spectral nudging which results in a reduction in the mean DPE by 20 km. The intensity varies as well with different surface flux options, with surface flux option 1, where the momentum roughness length is expressed using a combination of two roughness length formulas (Green and Zhang, 2013), in which the 570 first is Charnock (1955) plus a constant viscous term and the second is the exponential expression from Davis et al. (2008) https://doi.org/10.5194/nhess-2021-400 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License. with a viscous term (as cited by Kueh et al., 2019) simulating better intensities than the other two options (default surface flux option and surface flux option 2). The use of boundary conditions from different ensemble members also resulted in variations in the simulated tracks and intensities, but still within the range of variability of the different parameterization experiments.
The use of the KF convective scheme and a more reasonable surface flux option (sf1) can help improve the simulated intensity. 575 While the use of the TK convective scheme and application of spectral nudging can improve the track simulation.
This study is part of an on-going effort to investigate the effect of future climate on the intensity and track of selected destructive TC case studies in the Philippines such as Haiyan using a regional climate model. The resulting sensitivities to the cumulus schemes will be an important consideration in simulating the TC case studies with climate change forcing. Our 580 findings further stress the need for choosing the appropriate cumulus schemes and surface flux parameterization given its impacts on different TC characteristics, e.g. the KF scheme and surface flux option 1 for simulating better intensities of extreme TCs such as Haiyan, besides higher grid resolutions as noted in previous studies (Kueh et al., 2019, Li et al., 2018. The results presented here can also be used in further improving the value of downscaling for simulating intense TCs like Haiyan. These and future results will be useful in addressing the growing need to plan and prepare for, and reduce the impacts of future TCs 585 in the Philippines. As shown in this study, there are uncertainties associated with the use of cumulus parameterizations schemes, spectral nudging and surface flux parameterizations. To cover these uncertainties, the use of ensemble simulations can be applied. For operational applications, an ensemble of cumulus parameterizations can be used to take into account the uncertainty in the track and intensity of simulation intense TCs. This study can facilitate research on regional climate modeling to improve simulations of intense TCs like Haiyan. Furthermore, it is important to study RCMs with a model resolution less 590 than 5 km that can be extremely useful in simulating TCs and associated rain. Li et al. (2018) suggested that a 2-km convectionpermitting resolution is needed to reproduce intense TCs such as Haiyan. Other model parameterizations such as cloud microphysics and planetary boundary layer, and ocean coupling may help further improve the intensity simulations of extreme TC such as Haiyan but are beyond the scope of this paper. Simulations using a higher resolution convection-permitting model is needed. Additional simulations and further investigations on these aspects, as well as for other similar TCs will be useful. 595 Appendices none

Data availability
Simulation data are stored at the JASMIN data storage facility and is available upon request from the corresponding author.
Author contribution: 605 RJD designed the experiments with guidance from GB, PLV and KH. RJD performed the simulations and analysis with inputs from all co-authors particularly the interpretation of the results. RJD wrote the article with contributions from all co-authors. https://doi.org/10.5194/nhess-2021-400 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.