Response of damaging Philippines tropical cyclones to a warming climate using the pseudo global warming approach

The potential changes in the characteristics and damage potential of three of the most damaging tropical cyclone (TC) events (Haiyan 2013, Bopha 2012, Mangkhut 2018) in the Philippines have been simulated using the pseudo global warming (PGW) technique. Simulations were performed using the Weather Research and Forecasting model at 5 km resolution with cumulus parameterization (5 kmCU) and 3 km without cumulus parameterization (3 kmNoCU), with PGW deltas derived from a selection of the CMIP6 models. We found that re-forecasting the three TCs under future warming leads to more intense TCs, with changes in maximum wind of 4%, 3%, and 14% for the 5 kmCU runs, and 14%, 4%, and 12% for the 3 kmNoCU runs of Typhoon Haiyan, Bopha, and Mangkhut, respectively. The changes in track, translation speed, and size are relatively small. The TC cases have a higher impact potential in the future, as expressed by the cyclone damage potential index, ranging from ~ 1% to up to 37% under the SSP5-8.5 scenario. Based on the pre-industrial runs, climate change has had, so far, only a weak influence on TC intensity and not much influence on track, translation speed, and size. Simulations without convective parameterization show similar changes in the sign of the projected TC intensity response, but different signals of change in translation speed and size.


Introduction
The Philippines has high exposure to tropical cyclones (TCs) which often result in casualties and significant damage to property due to strong winds and flooding from storm surges and associated rainfall. The effect of climate change on TC activity and its characteristics is of major interest and practical importance in the Philippines, given the socioeconomic consequences produced by TCs. Most studies predict a decrease in the frequency of TCs in the future, but an increase in the intensity (high confidence) and number of intense TCs (low confidence), and TC-associated rainfall globally (Walsh et al. 2019;Knutson et al. 2020;Knutson et al. 2021) as well as in the Western North Pacific (WNP) Basin (Christensen et al. 2013). There is still no clear consensus on TC translation speed projections , which could lead to damage from accumulated and prolonged rainfall, as well as prolonged exposure to strong winds. In the Philippines, estimates are broadly consistent with global research results on TC activity, which predict a decrease in TC frequency in the future (Gallo et al. 2019), but an increase in the frequency of strong TCs (Walsh et al. 2016 andSugi et al. 2017) and an increase in TC-related rainfall (Daron et al. 2016).
Global climate models (GCMs) are critical tools for projecting future changes and trends in extreme weather and climate events. However, owing to their typically coarse resolutions and model biases, GCMs have problems to replicate observed TC properties. This is despite considerable breakthroughs in the use of high-resolution GCMs (Roberts et al. 2020;Manganello et al. 2014;Murakami et al. 2012) and global convection-permitting models in modeling TCs in recent decades (Gutmann et al. 2018;Judt et al. 2021;Kendon et al. 2021;Yamada et al. 2017) However, GCMs at high resolutions require very high to extensive computing resources. A useful alternative approach to increase the understanding of the changes in TC activity and their characteristics, for the purpose of simulating event-based, or seasonal to inter-annual variability and long-term climatological changes, is the use of downscaling. It encompasses both dynamical downscaling using Limited Area Models (LAMs) and statistical downscaling methods (Knutson et al. 2013;Xu et al. 2019;Chen et al. 2020;Emanuel et al. 2008Emanuel et al. , 2020. Statistical downscaling methods have been found to be good over the historical period since they are tuned with observations. However, the empirical relationships between largescale predictors (e.g., GCM-derived atmospheric parameters) and local predictands (e.g., winds, rainfall) are assumed to also hold in the future (Camargo and Wing 2016) which may not be correct. Dynamical downscaling with LAMs, used as regional climate models for long-term seasonal or climate simulations or as numerical weather prediction for short-term weather forecasting, can realistically capture the processes relevant to the formation and evolution of TCs but are still dependent on parameterizations (Seneviratne et al. 2021;Gallo et al. 2019;Daron et al. 2018). However, it is important to note, that the performance of dynamical downscaling using LAMs is strongly affected by the quality of the boundary forcing data i.e., output from GCMs (Holland et al. 2010;Liu et al. 2019).
There are different approaches to understanding the climate change effects on TC characteristics at higher resolution using dynamic downscaling (e.g., Lynn et al. 2009;Lackmann 2015). Adachi and Tomita (2020) provide a summary of these different approaches including the Pseudo Global Warming (PGW) approach. Lynn et al. (2009) performed current and future simulations of Hurricane Katrina (2005) using the PGW method and found decreases in the minimum central pressure but also decreases in the mean and maximum wind speeds, maybe due to the shift in Katrina's simulated track. The PGW technique allows us to examine the change of characteristics of historical TC cases under future and even past climate conditions (Takayabu et al. 2015;Parker et al. 2018;Patricola and Wehner 2018;Chen et al. 2020). Lackmann (2015) used convection-permitting simulations to look at changes in Hurricane Sandy (2012) due to the mean changes in temperature and humidity in the future and found decreases in the minimum central pressure. Nakamura et al. (2016) used the same technique for Typhoon Haiyan and showed that Typhoon Haiyan became more intense if only the sea surface temperature (SST) is changed and less intense if SST, atmospheric temperature (ATM) and relative humidity (RH) are changed in the future conditions. Parker et al. (2018) conducted PGW simulations for three TCs that made landfall in Australia and found that the TC cases are more intense (reduced sea level pressure and increased wind speeds) and have more rainfall. Patricola and Wehner (2018) also applied the same technique, but with higher resolution convection-permitting regional climate simulations for several TC case studies across the globe (including Typhoon Haiyan) and found that these TCs intensify in the future simulations and have greater rainfall amounts. Chen et al. (2020) looked at the potential future changes of three landfalling TCs in the Pearl River delta using the PGW technique and simulated the potential changes in storm surge activity. Their results showed increases in the peak intensities of the TCs and a corresponding increase in storm surge. Mittal et al. (2019) and Reddy et al. (2021) also conducted similar experiments on different TC cases in the Bay of Bengal. The former simulated Phailin (2013) and found a stronger and bigger storm in the future, with no change in translation speed. The latter looked at three TC cases and found increased winds and rainfall, reduction in translation speeds, and deepening of the TC core and higher damage potentials. Gutmann et al. (2018), on the other hand, performed a 13-year convectionpermitting simulation of 22 TC cases in the Atlantic Ocean and found increased maximum winds, lower central pressure, slower translation speeds, and higher precipitation rates, with varying degrees and magnitude of responses across the 22 TC cases. With the differences in changes in TC characteristics shown in previous studies in different basins, it is important to look at the potential response of TC characteristics to climate change in the WNP basin, particularly in the Philippines.
In addition, the PGW technique has also been used to investigate whether present-day climate change has already affected TCs in the Atlantic Basin (Patricola and Wehner 2018) and in Japan (Kawase et al. 2021). Using simulations from a 4.5 km convection permitting Weather Research and Forecasting (WRF) model, Patricola and Wehner (2018) found that climate change has so far increased rainfall of Hurricanes Irma, Katrina, and Maria but did not influence their intensities. On the other hand, Kawase et al. (2021) found that historical warming has already increased the intensity of Typhoon Hagibis (2019), along with its associated precipitation over Japan. So far, there is no consensus regarding whether climate change has already influenced tropical cyclone frequency and intensity due to the large natural variability and long-term data heterogeneity. Studies such as that of Patricola and Wehner (2018) and Kawase et al. (2021) may help provide evidence as to whether climate change has already influenced TC activity. We provide additional evidence towards understanding the link between climate change and TCs, primarily looking at how climate change has so far affected TCs in the Philippines.
Previous studies that used the PGW technique highlighted the importance of looking at multiple TC cases (Nakamura et al. 2016;Parker et al. 2018;Gutmann et al. 2018;Mittal et al. 2019) with varying TC intensities (Nakamura et al. 2016) in different seasons to identify outliers in TC response to future climate conditions (Parker et al. 2018). These past studies did not investigate the uncertainties associated with the use of different GCMs (Parker et al. 2018;Patricola and Wehner 2018;Mittal et al. 2019;Reddy et al. 2021) which is relevant in accounting for the uncertainty due to the potential range of sensitivities brought about by different models. The innovations in our study are: (1) investigating the different TC cases with observed varying intensity, landfall area, and months of occurrence; (2) the simulations were forced with initial and boundary conditions from several GCMs contributing to the Coupled Model Inter-comparison Project Phase 6 (CMIP6) (Eyring et al. 2016) and looking at different factors such as sea surface temperature, atmospheric temperature and relative humidity; (3) the WRF's simple mixed-ocean layer model was used in order to incorporate the atmosphere-ocean feedback during the passage of TCs; and (4) we tested the sensitivity to different cumulus schemes and initialization times. We have also attempted to identify the physical mechanisms driving the simulated TC responses and discuss relevant underlying uncertainties.
This study aims to analyze how the characteristics and potential impacts of the most damaging TC events in the Philippines might change under future climate conditions using dynamical downscaling with the high-resolution WRF, configured for the Southeast Asia region centered on the Philippines, as in Delfino et al. (2022a, b), and using the PGW technique. We aim to answer the following questions: • How might the most damaging TC events in the Philippines i.e., Typhoons Haiyan (2013), Bopha (2012) and Mangkhut's (2018) respond to past (pre-industrial) and future climate change perturbations derived from the latest CMIP6 GCMs? • How do the added PGW deltas (warming signals i.e., surface and atmospheric warming; with or without relative humidity) affect TC characteristics? • How much of the uncertainty in the TC responses is caused by the convective parameterization?
The rest of this paper is organized as follows: Sect. 2 presents a brief description of the TC cases, experimental design and methodology; Sect. 3 describes the simulated changes in the large-scale atmospheric variables (Sect. 3.1), and changes in the TC characteristics (3.2), discussion of the physical mechanisms behind the changes in TC characteristics (3.3) and cyclone damage potential (3.3) Sect. 3.2 is further separated in two sub-sections presenting the results of the different forcings under future climate (3.2.1) and the changes due to the different resolutions under past and future climate (3.2.2) The summary and conclusions are presented in Sect. 4.

Data
The European Centre for Medium-Range Weather Forecasts (ECMWF) Re-analysis 5th Generation (ERA5) is used for both the initial and boundary conditions for the simulations under the current climate. The data is available from https:// cds. clima te. coper nicus. eu/. It is the latest generation of reanalysis products produced by ECMWF with horizontal resolution of 31 km, hourly temporal resolution and 137 atmospheric levels (Hersbach et al. 2020). ERA5 utilizes the best available observational data from satellites and in-situ stations, which are quality controlled and assimilated using a state-of-the-art 4-dimensional data assimilation system (4D-VAR) (Isaksen et al. 2010) in and the ECMWF's Integrated Forecast System (IFS) Cycle 41r2.
The TC best-track information used here for verifying the simulations is taken from the World Meteorological Organization (WMO) subset of the IBTrACS (IBTrACS-WMO, v03r09) which was taken from the best-track data provided by the Japan Meteorological Agency (JMA). The CMIP6 GCM data were downloaded from https:// www. esgf-index1. ceda. ac. uk/ search/ cmip6-ceda/.
We used data from CMIP6 models for the far-future period of 2070-2099 and the worst-case/high-emission scenario from the suite of Shared Socio-economic Pathways (SSP) 5-8.5. The pre-industrial run covers the period between 1850 and 1899 and the current climate covers the period between 1970 and 2000. Under the pre-industrial climate, the PGW delta was calculated by adjusting the same variables from the historical 1970-2000 monthly mean minus the pre-industrial 1850-1899 monthly mean, then added to the current ERA5 initial and boundary condition. Four CMIP6 models were selected based on a subset of the CMIP6 GCMs (see Supplementary Fig. 1) analysed in Emanuel (2020) and that correspond to a range of Transient Climate Response (TCR) and Equilibrium Climate Sensitivity (ECS) (Meehl et al. 2020) which are expected to provide a good range of potential storylines for the TC cases . The TCR and ECS were used as proxies, to make sure that our simulations would experience a wide range of warming levels and of changes in environmental conditions. We have used these proxies, therefore, to provide a representative and relevant set of storylines considering the following: (a) small to large changes in magnitude of global warming; (b) magnitude of SST change in the region; and (c) patterns of SST change 1 3 in the region. The description of the selected models is presented in Table 1.

TC case studies
For this study, we have chosen the three of the most damaging TCs to have affected the Philippines between 1970 to 2020 (Lara 2020) -Typhoon Haiyan (2013), Typhoon Bopha (2012), and Typhoon Mangkhut (2018). Typhoons Haiyan and Bopha are also in the top five deadliest TCs in recorded history in the Philippines. These TC cases were also selected to represent the three main regions of landfall-Luzon for Mangkhut, Visayas for Haiyan, and Mindanao for Bopha.

Experimental design and model setup
The three TCs were simulated using the Advanced Research Weather Research and Forecasting (WRF-ARW) model version 3.8.1 (Skamarock et al. 2008). The WRF model has been used extensively in simulating TCs in the Philippines (Flores 2019;Aragon and Pura 2016;Spencer et al. 2012;Islam et al. 2015;Lee and Wu 2018;Cruz and Narisma 2016). The model setup and configurations, based on previous sensitivity experiments (see Delfino et al. 2022a, b;Bagtasa 2021), are detailed in Table 3. Data from ERA5 was used as initial and boundary conditions for the control simulations and the augmented initial and boundary conditions obtained from adding the PGW deltas were used for the future simulations. By comparing these future simulations with the control runs, we can therefore infer the responses of the TC cases to future climate conditions. Simulations under the current and future climate (with surface and air temperature, and relative humidity deltas) were performed with four different initialization times (with 6-hourly intervals) to construct ensemble members, performed using a two-way nest of 25 km and 5 km outer and inner domain, respectively, where the cumulus parameterization is turned on in both outer and inner domains. This experimental configuration is referred to as 5 kmCU (Table 4). Two-way nesting configuration has been used in previous TC case studies in the Philippines (e.g., Delfino et al. (2022a); Mori et al. (2014), Takayabu et al. (2015), Nakamura et al. (2016) and other PGW studies in other TC basins (Parker et al. 2018;Davis et al. 2008;Mittal et al. 2019;Reddy et al. 2021). We also performed simulations at a convection-permitting horizontal resolution of 3 km where the cumulus parameterization is turned off in the 3 km inner domain, referred to as 3 kmNoCU, to account for the uncertainty in the use of cumulus parameterization

Pseudo-global warming (PGW) technique
This study used the PGW method introduced by Schär et al. (1996) as a surrogate climate change method and later used in various studies Sato et al. 2007;Knutson et al. 2008;Nayak and Takemi 2020). The PGW technique adds a climate perturbation signal to the presentday conditions for the period of interest. The initial and boundary conditions are generated by combining 6-hourly reanalysis data (ERA-5) and climatological monthly mean perturbations that are extracted from the GCMs. The result is the forcing data for the alternative climate, called the PGW condition or climate change delta. It can then be used as the WRF initial and boundary conditions for climate simulations of the TC case studies. The PGW technique is one method of evaluating the impacts of different climate change scenarios such as global warming on past TCs. Other methods are presented in Adachi and Tomita (2020), a summary of which is presented in Supplementary Table 1, and a new approach is being proposed by Dai et al. (2020) to combine the transient weather signal from one GCM simulation with the monthly mean climate states from multi-GCM ensemble mean for the periods of interest. In most studies, the downscaled PGW condition is compared to the historical or present/currentday condition developed from reanalysis data to determine the projected change in the future climate for example. It is also well-suited for a storyline approach, which investigates the impact on an event intensity, as compared to the riskbased approach which investigates the impact on an event frequency (Shepherd et al. 2018). The PGW delta was calculated by subtracting the monthly means (i.e., November for Haiyan, December for Bopha, and September for Mangkhut) of the current simulations from the future climate projected by each of the GCMs. The calculated PGW deltas for each variable were then added to the 6-hourly ERA-5 initial and boundary conditions for the three TC cases to build the PGW conditions. The PGW delta was calculated in the initial and lateral boundary conditions for the following variables: surface (land and sea) and air temperature, pressure (surface and sea level), geopotential height and relative humidity (RH), as used in previous PGW studies (Patricola andWehner 2018, Parker et al. 2018). There is disagreement in the literature regarding whether the RH should be allowed to change in future simulations. Some PGW studies (e.g., Parker et al. 2018;Mittal et al.2019;Chen et al. 2020;Reddy et al. 2020;Lackmann et al. 2015) have applied a delta to the RH, but there are also examples where RH is not changed (e.g., Kawase et al. 2009;Patricola and Wehner 2018) where it is argued that: (1) an inconsistency is created when adding a broadly uniform, GCM-derived specific humidity delta onto a specific synoptic pattern and (2) RH is generally expected to remain constant in the future. Based on the results of Parker et al. (2018), the inclusion or exclusion of the RH delta does not affect the overall results. For this study, we have decided to experiment with changes in RH, in order to further understand its effects on TC response (FULL experiments).

TC tracking method
The simulated track and intensity values were obtained every 6 h using the TRACK algorithm (Hodges et al. 2017) as used in Hodges and Klingaman (2019). TRACK determines TCs as follows: first the average of the relative vorticity at 850-, 700-, and 600-hPa levels is obtained and averaged. The field is then spatially filtered using 2D discrete cosine transforms equivalent to T63 spectral resolution and the large-scale background is removed. The feature points are determined by first finding the grid point relative vorticity maxima > 5.0 × 10 -6 s -1 which are then used as starting points for a B-spline interpolation and steepest ascent maximization method to determine the off-grid feature points (Hodges 1995 as cited by Hodges and Klingaman 2019) that results in smoother tracks. The tracking is performed by first initializing the feature points into tracks using a nearest neighbour method and then refining them by minimizing a cost function for track smoothness subject to adaptive constraints (Villafuerte et al. 2021). The tracking is done for the entire simulation period. After completing the tracking, other variables are added to the tracks, including the maximum 10-m winds and minimum central pressure at full resolution. This is done by searching for the maximum 10-m winds within a 6° geodesic radius, and for the true minimum mean sea level pressure within a 5°-radius using the B-splines and minimization method (Hodges and Klingaman 2019). In addition, the TC vertical structure or profile was also taken using the TRACK algorithm by calculating the tangential, and radial winds at peak intensity at height levels between 0 and 18 km on a radial grid with geodesic radius of 10° (Bengtsson et al. 2007) for each simulation.

Cyclone damage potential
The potential damage that can be caused by the TCs is calculated using the cyclone damage potential (CDP) index Holland et al. 2016), formulated using critical cyclone parameters that involve the maximum wind (or energy dissipated at the surface), size (radial extent) and translation speed, and is given by: where Vmax and TS are the maximum surface wind speed and translation speed, respectively, both expressed in ms −1 . Vmax here is the maximum of 10 m wind speed over the inner domain. The translation speed of the TC for the current time step is calculated as a centered difference. The radius of the hurricane winds (HFWrad) is calculated by determining the radial distance between the storm center and the outermost point from the storm center where the wind speed is equal to 33 ms −1 . The CDP index is computed only if Vmax is greater than 33 ms −1 . If TS < 2.6 ms −1 , then TS is set to be as 2.6 ms −1 . The CDP index value ranges from 0 to 10. A higher CDP indicates an increased risk of TC-related damage costs, particularly damage brought about by wind ).

Large-scale conditions under future and pre-industrial climate
Based on the future climate change signals calculated from these four GCMs, the mean SST change in the experiment domain for the far future is projected to be between 1.  (Fig. 1g-i) and MPI ( Fig. 1j-l) models. The vertical distribution of air temperature is also shown on Fig. 2a-c for November, December, and September, respectively, which shows that the air temperature will be ~ 3-6 °C warmer on the average, with stronger warming in the upper level (near 250 hPa), leading to the atmosphere being more stable in the future under SSP5-8.5 (Fig. 2). This bulge of warming in the upper troposphere is consistent with previous studies (Hill and Lackmann 2011;Knutson and Tuleya 2004). The RH in the troposphere, particularly in the mid-troposphere, shows minimal change in the future climate under the SSP5-8.5 scenario ( Fig. 2d-f).

TS
In the pre-industrial climate, the change in SST ranges between − 0.03 to more than − 0.2 °C (see Supplementary Fig. 2), the changes in SST and the change in air temperature and relative humidity are also shown in Supplementary Fig. 3. Nakamura et al. (2016) showed that atmospheric and surface temperature warming has opposing influence on Haiyan's intensity under future climate conditions. The results of the experiments with different forcing and GCMs are Fig. 1 The November (left panels), December (middle panels) and September (right panels) monthly mean sea surface temperature delta (°C) added to the boundary and initial conditions of the control runs for Typhoons Haiyan, Bopha, and Mangkhut, respectively, to create the future climate scenario change in far future (2070-2099) from CMIP6 models (CESM2-a-c; HadGEM3-d-f; MIROC6-g-i; and MPI-j-l) according to the SSP5-8.5 scenario relative to the historical period . The domain-averaged delta SSTs are indicated in the lower right hand corner Fig. 2 The November (left panels), December (middle panels) and September (right panels) monthly mean air temperature delta (°C) (top panels) and relative humidity delta (%) (bottom panels) added to the boundary and initial conditions of the control runs for Typhoons Haiyan, Bopha, and Mangkhut, respectively, to create the future climate scenario change in far future (2070-2099) from CMIP6 models (CESM2, HadGEM3, MIROC6 and MPI) according to the SSP5-8.5 scenario relative to the historical period        Knutson et al. (2021).

Changes due to different forcings and GCMs
Overall, there is a range in the magnitude of changes due to the different GCM forcings which highlights the uncertainties associated with the use of different GCMs. However, the changes in intensity in the simulations forced by higher SST deltas (CESM2 and HadGEM3) and lower SST deltas (MPI and MIROC6) do not show systematic differences e.g., the simulations forced by HadGEM3 which has the highest SST delta (mean of 3.27, 3.12 and 3.65˚C for Haiyan, Bopha and Mangkhut, respectively) resulted in differences in maximum wind speed of 7 ms −1 , 7 ms −1 , and 10 ms −1 for Haiyan, Bopha and Mangkut, respectively under the FULL experiments. The simulations forced with MIROC6 which has a relatively lower SST delta (mean of 2.16, 2.13 and 2.12˚C for Haiyan, Bopha and Mangkhut, respectively) resulted in differences in maximum wind speed of 4 ms −1 , 9 ms −1 , and 14 ms −1 for Haiyan, Bopha and Mangkut, respectively, also under the FULL experiments.
The boxplots in Fig. 3 show that there is a systematic shift towards higher intensities in the lifetime maximum intensity over the simulated lifecycle of the storms under future climate conditions for all experiments, with a reduction in central pressure and increase in the maximum winds, which are most apparent in the more intense Mangkhut and Haiyan, and not so much for Bopha. The minimum central pressure (maximum surface winds) reached 890 hPa (73 ms −1 ) for Haiyan under the future climate (FULL) compared to 918 hPa (65 ms −1 ) in the current climate. Mangkhut reached up to more than 70 ms −1 in the future climate (FULL) runs compared to 917 hPa (61 ms −1 ) in the control.
To evaluate the response of the TC cases to climate change, the simulated TC tracks under current climate were first verified with observed TC tracks from the IBTrACS (Knapp et al. 2010) dataset. Indeed, the simulated tracks do not deviate substantially from the observed tracks under the  Supplementary Fig. 4), with the future climate simulations having a general tendency for the TCs to track more northwards than in the current climate (Bopha) and to recurve north and not make landfall in any part of the Philippines (Haiyan and Mangkhut). Despite some small deviations, the TC tracks under the future climate simulations are still within 3° away from the current climate simulation, a threshold that was found acceptable for this kind of analysis by Patricola and Wehner (2018).
We also looked at different parameters of TC size, which is defined here as the radius of different wind thresholds as used in Radu et al. (2014) i.e., gale-force winds (GFW, 17.5 ms −1 ), damaging-force winds (DFW, 25.7 ms −1 ), and typhoon-(or hurricane-) force winds (TFW, 33 ms −1 ). The SFC-only experiments present an increase of up to 112% in the radius of maximum winds (RMW) and 60% in TFW. While the changes in size under the SFC + PLEV and FULL experiments are relatively small, with changes ranging from -7% up to 12% for all three TC cases.

Changes due to different resolutions in different periods
To investigate the uncertainty in the response of the TCs to climate forcings caused by convective parameterization, which was found to be the source of uncertainty in the current simulations of the TC cases (Delfino et al. 2022a), additional simulations of the three TC cases at a convectionpermitting horizontal resolution of 3 km (3 kmNoCU) were performed. In addition, the three TC cases were also simulated under the pre-industrial climate to provide evidence as to whether CC has already influenced the TC characteristics. We have found that the simulations using the 3 kmNoCU have almost the same response to that of the 5 kmCU simulations i.e., the TCs have relatively lower minimum central pressure and higher maximum surface winds in the future under the SSP5-8.5 scenario over the lifetime of the TCs (Fig. 4), and all throughout the simulation period (Figs. 5 and 6). The future response in minimum central pressure ranges from -0.46%, -0.07%, and -1.20% for the 5 kmCU, and -0.07%, -0.43%, and -1.11% for the 3 kmNoCU for Typhoon Haiyan, Bopha and Mangkhut, respectively. The mean (spread) percent change in maximum surface wind ranges between 4.09% (0 to 7%), -2.68% (-7 to 2%), and 15.57% (13 to 20%) for the 5 kmCU runs, and 13.75%, 4.28%, and 11.62% for the 3 kmNoCU runs of Typhoon Haiyan, Bopha and Mangkhut, respectively.
There are also relatively bigger spreads in the lifetime intensity of the TCs in the 3 kmNoCU compared to the 5 kmCU runs. This might be due to the relatively lower intensity at the start of the simulations (first 24 h) in the 3 kmNoCU runs compared to the 5 kmCU run (Fig. 5). During the TC's lifetime, the 3 kmNoCU runs achieve a deeper and more intense state. Figure 4 shows that the minimum of the central pressure is relatively lower in the 3 kmNoCU runs, and the maximum of the wind speeds is higher than in the 5 kmCU runs. The spread in the boxplots is largely due to the lower intensities in the 3 kmNoCU runs at the start and end of the simulations.
Of all the simulations, the 3 kmNoCU runs under future climate shows significant increases in intensity (at 1% and 5% level) for all three cases with p-value of 0.0003 and 3.712e-05 for Haiyan's minimum central pressure and maximum winds, p-value of 0.0317 and 0.0003 for Bopha's minimum central and maximum surface winds, and p-value of 0.0290 and 0.0738 for Mangkhut's minimum central and maximum surface winds. Haiyan (at 10% level) and Mangkhut (at 1% level) also showed significant intensity increases under the future climate in the 5 kmCU runs (Table 7).
For the three TC cases, the ensemble of the minimum central pressure and maximum surface winds are almost similar (Figs. 4,5) under the current and pre-industrial climate conditions in both 5 kmCU and 3 kmNoCU runs, particularly for Typhoon Haiyan and Mangkhut, with some deviations after landfall for Bopha. However, in terms of peak intensity, in the more intense TCs Haiyan and Mangkhut, CC may have already had a weak influence on intensity in terms of relatively weaker winds than in the current climate (Fig. 3). In particular, the difference between the current and pre-industrial climate for minimum SLP (maximum wind) are 10 hPa (-4 ms −1 ) and -9 hPa (4 ms −1 ) for Haiyan and Bopha, respectively, while no mean difference for Mangkhut. For the 3 kmNoCU runs, Haiyan has a -1 hPa (2 ms −1 ), Bopha has -19 hPa (10 ms −1 )-significant at 1% level, and Mangkhut -3 hPa (2 ms −1 ) -significant at 5% and 10% respectively, current-pre-industrial difference (Table 7), which indicates that the TCs are relatively weaker under the pre-industrial climate than under the current climate simulations. Figures 7 and 8 show the vertical profiles of the acrosstrack tangential winds at peak intensity over a 10° radius for all the simulations for the 5 kmCU and 3 kmNoCU simulations, respectively. The profiling is done at peak intensity at height levels between 0 and 18 km and the radial grid extends to 10° for each simulation. The intensification of the TC cases under the future climate scenarios is also shown by the vertical cross-section of the composite of azimuthally averaged winds. Under the SSP5-8.5 scenario, Figs. 7 and 8 show the maximum tangential wind speed increases from 55 to 65 ms −1 to up to more than 70 ms −1 . This is similar to the result in Manganello et al. (2014) which projects an increase of around 5 ms −1 in tangential wind speed on the average for TCs over the WNP. The experiments under future climate scenario also project stronger outward and upward expansion of the region of stronger winds. The radial winds also tend to be stronger (inflow) near the surface. Under the pre-industrial climate, the vertical distribution of tangential winds is similar to that of the simulations under the current climate. It is also worthwhile to note here that the vertical distribution of strong tangential winds (> 50 ms −1 ) increases more in the 3 kmNoCU (Fig. 8) than in the 5 kmCU runs (Fig. 7), for all time periods and three TC cases. This is consistent with the results of an earlier study by Kueh et al. (2019) that showed an increase in the vertical extent of stronger tangential wind speeds, a decrease in the outward slope of the maximum wind axis, and a contraction of the eyewall when model grid resolution is increased.
To have a fair comparison between the TC cases under different climate conditions the simulated tracks must be robust and not deviate very much from the observed tracks. Indeed, for both the 5 kmCU and 3 kmNoCU simulations, there are no significant deviations in the track for both pre-industrial and future runs (Supplementary Figs. 4,5 and 6). In addition, there are no substantial changes in translation speed of all pre-industrial runs in both 5 kmCU and 3 kmNoCU experiments, with percent changes in translation speed ranging from ± 2 m/s. In terms of future projected translation speed, however, there is a small but opposing response between the 5 kmCU and 3 kmNoCU runs, with slower TCs in the 3 kmNoCU (-4 to -24%) Fig. 7 Radius (in °s) -height (in km) cross sections of tangential (across-track) for Typhoons Haiyan, Bopha and Mangkhut under preindustrial (left panels), current (middle panels) and future (right pan-els) climate conditions for the 5 kmCU simulations. The number in the right-hand corner indicates the maximum 10-m peak wind intensity in ms −1 and faster TCs in the 5 kmCU runs (+ 2 to 9%) under the future climate compared to the current climate simulations. We also found small changes in size, ranging between -7 to 12% in all future runs for all the TC cases in both the 5 kmCU and 3 kmNoCU runs.

Discussion
In this section we investigated some of the thermodynamic and dynamic aspects from the current, future, and past simulations associated with the changes in intensity, size, and speed of the different TC cases.

Intensity
The projected increase in intensity is consistent with the Potential Intensity Theory introduced by Emanuel (1987) based on the Carnot cycle heat engine. According to Emanuel (1987), there are two primary factors that control how strong a TC's winds can get -(1) ocean surface temperature and heat content and (2) the temperature and moist state of the atmosphere, thus a TC's potential intensity depends on the moist thermodynamic state of the atmosphere and ocean. Other factors can keep a TC from reaching its potential intensity like vertical wind shear, dry air entrainment and SST cooling from air-sea interaction. If CC increases TC Fig. 8 Radius (°) -height (in km) cross sections of tangential (acrosstrack) for Typhoon Haiyan, Bopha and Mangkhut under pre-industrial (left panels), current (middle panels) and future (right panels) cli-mate conditions for the 3 kmNoCU simulations. The number in the right-hand corner indicates the maximum 10-m peak wind intensity in ms −1 potential intensity, and other factors counteract this, then it is expected (in theory) that TCs should shift to stronger intensities. Studies like Patricola and Wehner (2018), Bhatia et al. (2018) and Kim et al. (2014) are consistent with this theory.
Based on the results presented in Sect. 3.2.1, there are substantial changes in the intensities in the SFC-only experiments, and more modest increases in the SFC + PLEV and FULL experiments. Such a large intensification is driven by more heat flux supplied from a warmer sea surface in the future. This result agrees with previous studies that demonstrate the important role of warm SST on TC intensification (Delfino et al. 2022b;Chen et al. 2020;Nakamura et al. 2016). On the other hand, as can be seen in the SFC + PLEV and FULL experiments, atmospheric warming has a large negative impact on TC intensity, but the TC cases still increased intensity under the future climate. Some previous studies have also pointed out that atmospheric warming itself can weaken TCs, due to the associated increase in atmospheric stability (Chen et al. 2020;Nakamura et al. 2016). Additionally, although the relative humidity is slightly increased under the future climate, the FULL experiments (with relative humidity changes) show no significant difference from the SFC + PLEV experiments, consistent with Parker et al. (2018) and Chen et al. (2020).
From the results presented in Sect. 3.2.2, we have seen increases in intensity under the future simulations. Based on maximum potential intensity theory (Emanuel 1987) and recent studies (Patricola and Wehner 2018, Nakamura et al. 2016, Chen et al. 2020, Mittal et al. 2019, Parker et al. 2018, TC intensity is expected to increase due to favorable conditions such as warmer ocean temperature, unstable atmosphere with a moist mid-troposphere, and weak vertical wind shear (Wu et al. 2022). Figures 9 and 10 shows increases in both the water vapor mixing ratio and latent heat flux under the SSP5-8.5 scenario in the future climate. The water vapor mixing ratio increase is 32% (66%), 27% (29%), and 26% (29%) under the future climate condition in the 5 kmCU (3 kmNoCU) experiments (Table 8). The increase in latent heat flux is 10% (30%), 16% (30%), and 5% (1%) under the 5 kmCU (3 kmNoCU) for Typhoons Haiyan, Bopha and Mangkhut, respectively (Table 8). This is also similar to the results of Gutmann et al. (2018) which showed that the increases in water vapor and thus latent heat feedbacks lead to increases in TC intensity. In the future climate, TCs may have a deeper TC core due to enhanced latent heating from increased precipitation and a significant decrease in the vertical wind shear (Mittal et al. 2019). The intensification may be driven by more latent heat flux supplied from a warmer ocean in the future (Chen et al. 2020;Nakamura et al. 2016) and due to the change to the SST gradient (Parker et al. 2018).
The observed changes in the vertical profile of temperature in a moist environment with enhanced SST will increase the convective available potential energy (CAPE) (Trenberth 2005). The CAPE values are also calculated, where a notable increase is seen in the far future scenarios for all Fig. 9 Simulated water vapor mixing ratio (g/kg), averaged over the domain, from the 5 kmCU (top panels) and 3 kmNoCU (bottom panels) experiments under pre-industrial (black), current (blue) and future (red) climate for the different TC case studies-(a, d) Haiyan (left), (b, e) Bopha (center) and (c, f) Mangkhut (right). Solid lines denote the mean, and the spread denotes the simulated intensity from the four ensemble members initialized at different times the TCs considered. There is a 41% (45%), 37% (51%), and 44% (36%) average increase in CAPE in the far future SSP5-8.5 scenario for the 5 kmCU (3 kmNoCU) simulations for Typhoons Haiyan, Bopha and Mangkhut, respectively. This increase in CAPE values in future climate is consistent with previous studies (Knutson and Tuleya 2004;Parker et al. 2018). The higher CAPE enhances the potential for more vigorous atmospheric convection and further intensifies TCs (Knutson and Tuleya 2004;Reddy et al. 2019). A warmer and moist environment with high SSTs is observed in the November (2.73) and September (2.75) climate change signal than in December (2.63) under the SSP5-8.5 (see in Fig. 2). Hence, TCs Haiyan and Mangkhut are seen to intensify more with substantial relative changes in CAPE compared to Bopha in the future climate (Fig. 11).

Track and translation speed
Based on the results presented in Sect. 3.2.1 above, there are substantial changes in the TC track under the SFC-only experiments for all three TC cases ( Supplementary Fig. 4).
We have also found that the SFC-only experiments, where the tracks recurved, had an enhanced westerly steering flow as well as a retracted Western North Pacific Subtropical High (WNPSH). This is consistent with previous studies, which showed that the TC track can be altered due to the changes in TC size as well as in the retraction of the WNPSH (Sun et al. 2015(Sun et al. , 2017 The observed and future changes in TC translation speed remain uncertain Yamaguchi et al. 2020), primarily due to lack of studies and clear evidence (Kossin 2021). In this study, we found minimal changes (± 2 ms −1 ) in translation speed with opposing signals from the 5 kmCU (faster) and 3 kmNoCU (slower). However, the changes in the large-scale atmospheric circulation patterns as well as the recurrent modes of climate variability affecting the WNP are not considered here, thus, the future changes in tracks and translation speeds are not fully accounted for.

Size
Larger TC size is found to be due to an increase in environmental humidity and temperature in the lower-troposphere leading to an increase in environmental CAPE (Mittal et al. 2019). Although the changes in mid-tropospheric RH are  minimal (between 1 to 13%) in all three cases, the RH are higher under future climate conditions relative to the current climate. This might explain the relatively small changes in size. There is also a relatively similar response between 5 kmCU and 3 kmNoCU experiments. In addition, there is also an increase in the mid-tropospheric RH in all three cases Haiyan, Bopha, and Mangkhut. The high mid-tropospheric RH supports TC development and intensification by restricting the negative influence on convection (Gray 1998). Figure 12 and Table 8 show the summary of the percent changes in the TC environment under the pre-industrial and future climate, relative to the current climate simulations. The changes in latent heat flux for the 5 kmCU (3 kmNoCU) runs ranges between 5% (5%), -7% (-5%), and 1% (-1%) for Typhoons Haiyan, Bopha, and Mangkhut under the pre-industrial climate, while under SSP5-8.5 scenario in the future climate, the latent heat flux changes range from 10% (30%), 16% (30%), and 5% (1%) for the 5 kmCU (3 kmNoCU) runs of Typhoons Haiyan, Bopha, and Mangkhut, respectively. The water vapor mixing ratio also increased in the future with changes ranging from 26 to 35% for the 5 kmCU runs and 29% to 66% for the 3 kmNoCU runs. There is a general decrease in the vertical wind shear which ranges from -2% (-4%), -2% (-3%), and -4% (-6%) for the 5 kmCU (3 kmNoCU) runs for Typhoons Haiyan, Bopha, and Mangkhut, respectively. These changes in the TC environment, particularly the increase in latent heat flux and water vapor mixing ratio; and the reduction in vertical wind shear, appear to be the main contributors to the further intensification of the TC cases in the future climate simulations. The increase in atmospheric convective instability in the TC outer region below the middle troposphere, which facilitates the local development of grid-scale ascending motion, low-level convergence, and the acceleration of tangential winds, may be responsible for the TC size response to ocean warming. A rise in environmental humidity and temperature in the lower troposphere leads to an increase in environmental CAPE, resulting in larger TC size (Mittal et al. 2019).

Cyclone damage potential
Overall, there is higher damage potential in the future for all three TC cases and in both the 5 kmCU and 3 kmNoCU simulations (Fig. 13). The increase in CDP ranges from 3.56, 2.43, and 5.50 for the 5 kmCU runs; and 7.55, 2.99, and 7.21 for the 3 kmNoCU runs under the far future SSP5-8.5 scenario. Under the pre-industrial (current) climate, the TC cases' CDP ranges from 3.43 (3.53), 1.97 (2.43), and 4.18 (4.25) for the 5 kmCU runs; and 4.65 (5.51), 1.67 (2.40), and 5.70 (5.46) for the 3 kmNoCU runs for Typhoons Haiyan, Bopha, and Mangkhut, respectively. The increase in CDP is primarily due to the increase in maximum winds, with small changes in both size and speed. The changes in CDP are generally larger in the 3 kmNoCU runs than in the 5 kmCU runs with changes ranging from + 1%, + 0.17%, and + 30% for the 5 kmCU runs and + 37%, + 25%, and + 32% in the 3 kmNoCU runs for Typhoons Haiyan, Bopha, and Mangkhut, respectively. Under the pre-industrial climate, however, there are minimal changes in the CDP for Typhoons Haiyan and Mangkhut, ranging between -3 to -15% (preindustrial minus current), but for Typhoon Bopha the change in CDP is -19% (-30%) in the 5 kmCU (3 kmNoCU) runs showing that Typhoon Bopha had more damage potential under the current climate as compared to if it happened in the past.
The percent changes in the damage potential of the different TC cases relative to the current are presented in Fig. 14. Typhoons Haiyan, Bopha, and Mangkhut's maximum winds increased under future climate conditions in both 5 kmCU and 3 kmNoCU simulations. Changes in size are relatively small and have mixed signal of responses, with increases/decreases for different TC case and in 5 kmCU and 3 kmNoCU simulations. The 5 kmCU experiments led to a change of + 30%, -11%, and -2% for Haiyan, Bopha and Mangkhut, respectively and the 3 kmNoCU experiments led to 0%, 6.67%, and 0% change for Haiyan, Bopha and Mangkhut. Under the pre-industrial climate, Haiyan did not show any change in the TFW and DFW under the 3 kmNoCU and 5 kmCU simulations, Bopha showed minimal changes, but Mangkhut showed up to -13% and -2% change in TFW under the 3 kmNoCU and 5 kmCU simulations (Fig. 14).
We also computed the changes in the translation speed of the TC cases since the future changes in speed remain uncertain. For the translation speed, we find an opposite signal of responses i.e., the TC cases move slower in 3 kmNoCU with -24%, -10%, and -4% and faster in the 5 kmCU simulations with 6%, 10%, and 4% change in speed for Typhoons Haiyan, Bopha, and Mangkhut, respectively (Fig. 14). Under the pre-industrial climate, the TCs response in terms of speed have changes of 18% for Haiyan, -1% for Bopha, and -2% for Mangkhut in the 3 kmNoCU runs. Additional work is needed to achieve more generalized conclusions over the TC response in terms of size and speed.

Summary and conclusions
The Weather Research and Forecasting (WRF) Model was used to simulate the potential impacts of large-scale environmental changes due to global warming on the simulation of three of the most damaging TC events in the Philippines. A set of simulations using 5 km resolution inner domain with cumulus scheme (5 kmCU), and a 3 km resolution domain without cumulus parameterization (3 kmNoCU), for four different initial times in each TC case, were able to capture the observed track and intensity of Typhoons Haiyan, Bopha, and Mangkhut.
The PGW technique was used to explore the response of the TC cases to historical and future warming climate. By using forcings from four different CMIP6 GCMs, the current weather patterns during the period of the TC events were re-simulated under different future (past) climates. The current, pre-industrial and future were computed from the monthly mean values obtained from the different GCMs, with the SSP5-8.5 scenario in the future. Based on the future climate change signals calculated from these four GCMs, the mean SST change is projected to be between 1.89 °C to 3.65 °C warmer and the air temperature will be ~ 3 °C to 6 °C warmer on the average, with stronger warming in the upper level (near 250 hPa), leading to the atmosphere being more stable in the future under SSP5-8.5 scenario. The relative humidity in the troposphere, particularly in the midtroposphere, shows a minimal change in the future climate under the SSP5-8.5 scenario.
Based on the results for the far future SSP5-8.5 scenarios, we have found that the intense TC events would be even more intense in the future, with significant increases in maximum peak winds and reduction minimum central pressure in both 5 kmCU and 3 kmNoCU runs. Previous studies that have done similar kinds of work for TCs from different basins (Patricola and Wehner 2018), for three TC cases in the Pearl River Delta (Chen et al. 2020), Bay of Bengal (Mittal et al. 2019), and Australia (Parker et al. 2018) has also reported mostly increases in peak intensity, particularly maximum winds. There are also higher rates of intensity change if we consider surface warming only, while atmospheric warming offsets the intensification caused by surface warming only, consistent with Nakamura et al. (2016) and Chen et al. (2020). Our study further found that in the future, the increase in intensity is mainly due to warmer temperatures, higher latent heat fluxes, increased water vapor mixing ratio, and less vertical wind shear. In addition, under the future climate, the more intense TCs like Haiyan and Mangkhut would have a deeper TC core or greater vertical extent of stronger winds primarily due to enhanced latent heating.
The results also show relatively small changes in track and translation speed in the future simulations compared to the current simulations, due to small changes in the TC steering flow. In fact, the TCs are slightly faster (slower) with an average of + 6% (-24%), + 10% (-10%), and + 4% (-4%) for Typhoons Haiyan, Bopha and Mangkhut in the 5 kmCU (3 kmNoCU) simulations compared to the current runs. Our study also found small changes in size, ranging between -7 to 12% in all future runs for all the TC cases, which might potentially add insights to the general lack of studies on the response of TC size to future warming (Kossin 2021).
By comparing with the current simulations, the simulations under the pre-industrial climate showed that climate change has so far weakly influenced the intensity of the TC cases but did not have much influence on the TC cases' track, translation speed, and size. With the increase in maximum winds, the damage potential of Haiyan and Mangkhut are slightly increased in the current climate simulations than the pre-industrial climate simulations. More importantly, the TC cases are expected to further intensify with continued warming in the future under the SSP5-8.5 scenario simulations compared to the current climate. The 3 kmNoCU simulations have the same sign of projected changes in TC intensity with the 5 kmCU simulations, with different responses in size and speed among the simulations. These results suggest that convective parameterization introduces minimal uncertainty in terms of TC intensity, however, TC speed and size need further investigation.
It is very important to note here that there is a range of responses i.e., signal and magnitude of change from different TC cases and GCM forcings. This highlights the uncertainty due to the potential range of sensitivities brought about by different models. However, the changes in intensity in the simulations forced by GCMs with higher mean SST deltas (CESM2 and HadGEM3) and lower mean SST deltas (MPI and MIROC6) do not show systematic differences. This highlights the uncertainties associated with the use of different GCMs and TC cases, which have been lacking in similar studies that used the PGW technique to investigate the response of TCs to climate change. These past studies did not investigate the uncertainties associated with the use of different GCMs (Parker et al. 2018;Patricola and Wehner 2018;Mittal et al. 2019;Reddy et al. 2021). Past studies have highlighted the need for and importance of studying multiple TC cases (Nakimura et al. 2016;Parker et al. 2018;et al. Mittal et al. 2019) with varying TC intensities (Nakamura et al. 2016) and in different seasons, to identify outliers in TC response to future climate conditions (Parker et al. 2018), which we have attempted to address in this study. Therefore, the novelty in our study is in investigating (1) simulations forced with initial and boundary conditions from several GCMs; (2) the sensitivity to the use of different cumulus schemes and initialization times (Delfino et al. 2022a); (3) using WRF's simple mixed-ocean layer model in order to capture the atmosphere-ocean feedback during the passage of TCs; and (4) different TC cases with observed varying intensity, landfall area, & months of occurrence. In addition, this study also tries to address the lack of attention of past studies to the physical mechanisms involved (Wu et al. 2022). We attempted to identify the physical mechanisms driving the simulated TC responses and discuss relevant underlying uncertainties.
However, there are still some aspects of the results of this study that need to be interpreted with caution due to (1) the limitations of the PGW technique i.e. it does not address the uncertainties surrounding the different modes of climate variability that are relevant in TC activity; (2) the small samplesize of the TC cases; and (3) it does not account for the range of uncertainty associated with the use of different limited area models and reanalysis datasets as initial and boundary conditions. The use of a single limited area model (WRF) for three TC cases and different GCM forcings allowed us to assess the robustness of the TC responses to future and past climate more directly among the TC cases and the different GCM forcings. However, the uncertainty in relation to the use of different limited area models is not within the scope of this study and can be an important aspect to look at in future studies. Global and regional modeling initiatives such as Coordinated Regional Climate Downscaling Experiment (CORDEX) could help address this uncertainty. Whilst the PGW method is useful to provide more detail of TC properties and their possible changes than global simulations, it cannot account for the remote influence of the large-scale atmospheric circulation patterns or the recurrent modes of climate variability affecting the WNP that can affect the steering flow of the TCs. This means that any changes seen in the tracks and translation speeds must be considered with caution. In addition, any large changes in track in the SFConly experiments mean the TCs will no longer experience the original environment modified by the climate change delta/SST anomalies.
In conclusion, based on our simulations of the different TC cases using the PGW approach, it is found that the most damaging TCs like Haiyan, Bopha, and Mangkhut will have higher damage potential in the future. The increase in the CDP ranges from ~ 1% to up to 37% in the future under the SSP5-8.5 scenario, primarily due to the increase in maximum winds. TCs of such intensity and damage potential in the future will have serious implications with the increasing exposure and vulnerability in the Philippines.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.