Climate change effects on the worst-case storm surge: a case study of Typhoon Haiyan

Human activities have changed the global climate and have affected some extreme weather events (Bindoff et al 2013). A typhoon (i.e., a tropical cyclone or TC) is one of the most potentially destructive extreme weather events. Estimation of the power dissipation index (cube of the maximum surface wind speed integrated for the whole lifetime of the event), a measure of the potential destructiveness of a TC over its lifetime, has already shown that TCs have become stronger in the past thirty years (Emanuel 2005). Moreover, the increase in this index is highly correlated with the observed increasing trend in sea surface temperature (SST) during the same period. However, studies disagree about whether changes in TC activity can be attributed to human influence, owing to insufficient observations and physical understanding (Bindoff et al 2013).

There is increasing interest in whether not only long-term trends of extreme events but also the characteristics of specific recently observed extreme events can be attributed to external drivers and natural climate variability. Simulations to estimate how predicted sea level increases up to 2050 or 2100 will change the probability of future inundations in the eastern coast of the United States on the scale of that caused by Hurricane Sandy (Sweet et al 2013) have indicated that the return period of a similar inundation is likely to decrease.

Pall et al (2011) proposed a new approach for quantifying the roles of external drivers and natural variability on specific extreme events by performing two types of ensemble simulations with an atmospheric general circulation model (AGCM). One took into account anthropogenic changes in the well-mixed greenhouse gas (GHG) concentration and the other did not. The first was forced by historical anthropogenic and natural forcing factors and observed SSTs, whereas in the second; historical changes in GHG concentrations and estimates of their effect on SST were omitted. Then, by comparing the results, they inferred how human activity influenced the likelihood and intensities of events. This approach of Pall et al (2011) is called Probabilistic Event Attribution (PEA), and has been applied in some research of recent extreme weather and climate events, such as heat waves and river floods (Otto et al 2012, Christidis and Stott 2014, Shiogama et al 2014, Wolski et al 2014).

The aim of this study was to estimate a robust signal of climate change in increasing the severity of a coastal hazard by Typhoon Haiyan, as an example of a worst case scenario in the present climate (e.g. Mori et al 2014, Lin et al 2014). Here we focus on the 'worst case scenario' in order to assess a physical upper bound that spawns disasters; we do not take into account the frequency of its occurrence. For this purpose, we performed two types of ensemble simulations of Typhoon Haiyan. The first used observed SSTs, atmospheric properties and GHG concentrations (ALL simulations), and the second ensemble used counterfactual natural external conditions (NAT simulations). Although the idea of our approach is partly based on that of PEA, it should be noted that we do not intend to investigate changes in the frequency of typhoons. To evaluate possible amplification effects due to climate change, we adopted ensemble prediction (Saito et al 2010) and applied dynamical downscaling from a global model (WEP: Weekly Ensemble Prediction system of Japan Meteorological Agency with the equivalent grid point resolution, 60 km) (Sakai 2009, Saito 2011) to a high-resolution Weather Research and Forecasting (WRF) model (Skamarock et al 2008) at the horizontal grid spacing of 1 km that was able to represent the strength of a category 5 TC (Gentry and Lackmann 2010, Kanada et al 2012). With the dynamical downscaled WRF 1 km models, we selected 16 experiments out of all 51 ensemble members.

When we discuss the activity of typhoons, frequency of genesis and tracks are important, as well as maximum intensity. To discuss the disaster prevention, we should handle these three metrics jointly. However, for the purpose of disaster mitigation, it is worth discussing at least the change in maximum intensity of the worst typhoon. Because Typhoon Haiyan was the most powerful typhoon to make landfall to date (Mori et al 2014), it is very important to discuss how the worst case event would change under the warming climate environment.

The paper has the following structure. The dynamical and ensemble downscaling procedure has been introduced in section 2. Section 3 describes the representativeness of the model typhoon, and also a signal of climate change appeared in the downscaling integrations. In section 4 we present some discussions, and section 5 is the conclusion.

The dynamical downscaling procedure adopted in this study is shown schematically in figure 1. The parent model was the Japan Meteorological Agency (JMA) weekly ensemble prediction model (WEP), which is the operational model used for weekly forecasts (Sakai, 2009). WEP, which is a global spectral model (GSM) with an equivalent grid point resolution of 60 km, was used to obtain the initial and boundary conditions, but the initial perturbations were calculated with a lower resolution GSM (equivalent horizontal grid resolution is 200 km). The ensemble had 51 members, prepared using the singular vector method (Buizza and Palmer 1995). Time integration was from 12 UTC on 4 November 2013 until 12 UTC on 11 November 2013.

Zoom In Zoom Out Reset image size

To avoid the occurrence of a shock at the lateral boundary, we used a cascade of nested models, downscaling finally to a WRF model with a 1 km grid resolution. The first downscaling step was to a regional model with a 20 km grid resolution (NHRCM20) (Sasaki et al 2011). To reproduce the TC track of the parent model in the downscaled model, we adopted spectral nudging above 15 000 m. The next downscaling step was to the regional model with a 5 km grid resolution (NHRCM05), and then two-way nesting from 3 km WRF to 1 km WRF models. In NHRCM05, the Kain-Fritsch convective parameterization scheme (Kain and Fritsch 1990, 1993) is used, and in the WRF models, cloud microphysical processes are explicitly handled.

Dynamical downscaling using 3 km/1 km WRF models has been done, but not for all ensemble members of WEP calculations because of the limitation of the computer resources. To cover the wide range of ensemble members, the following three cases are selected. They are, Case m02: where the track is nearest to the observed track (best track) around Leyte and Samar Island (124.8E). Case m11: where the minimum central pressure (MCP) of the typhoon calculated in NHRCM05 is lowest among all. Case p18: where the MCP of NHRCM05 is the highest among all. Because the timing of starting of the integration also affects MCP, the following three times are selected as the initial conditions. They are, #1000: 12UTC, 5, November, #1001: 18UTC 5, November, and #1002: 00UTC 6, November. Thus we have 9 experiments following the above conditions. Furthermore, to estimate the storm surge around Tacloban caused by Typhoon Haiyan, we select all cases where the typhoon passed around Leyte Island (125E), within the range of 50 km to the best track data. Seven cases (05 m, 12 m, 15 m, 21 m, 25 m, 06 p, 09 p) are found other than Case m02. For these seven cases, only #1001 (initial condition at 18UTC 5, November) has been calculated. As a result, we calculate 16 cases by using 3 km/1 km WRF models.

For the NAT ensemble simulations, we applied Pseudo Global Warming Downscaling (PGWD), in which the boundary conditions were assumed to be a linear coupling of the WEP data and the difference component of the climate change of air temperature between the middle of the nineteenth century and 2013.

To estimate the hypothetical natural SST, we removed linear trends in the monthly data from the HadISST (Hadley Centre Sea Ice and Sea Surface Temperature) data set (Rayner et al 2003) for 1870–2012 from the original SST data used in the ALL runs (Christidis and Stott 2014, Shiogama et al 2014). The difference between ALL-SST and NAT-SST in the target area (100°–180°E, 5°S–25°N) was between 0.2 K and 0.8 K (as shown in figure S2). To remove anthropogenic atmospheric warming, we subtracted the differences of atmospheric temperature between two 100-member ensemble simulations of the MIROC5 AGCM (Shiogama et al 2014, Watanabe et al 2010), performed with and without human influence, from the WEP data. The hypothetical anthropogenic SST signals of the MIROC5 ensemble are the same as in our regional climate model runs. We averaged the atmospheric temperature differences over the whole target area during November of 2010–2012 across 100 ensemble members. The averaged tropospheric and stratospheric temperature differences were +0.5 K and –1.5 K, respectively. We call this effect 'ATM'. For the downscaling calculations, we used GHG (CO₂, CH₄ and N₂O) levels in the 1850s reported by the Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC 2013).

The Surge-Wave-Tide coupled model (SuWAT) (Kim et al 2008) was forced by the 1 km WRF model calculation result. SuWAT is a fully coupled model of storm surge and ocean waves based on the non-linear shallow-water equation and spectral wave model, and takes into account atmospheric pressure, wind stress, and wave radiation stress (wave effects on current). Input data are surface wind and sea level pressure. We applied SuWAT to three domains with spatial resolutions ranging from 0.1° (D1) to 740 m (D3), with two-way nesting. The astronomical tide was excluded for computation and survey for simplicity, although the maximum astronomical tidal range is 0.7 m and was close to mean water level at the landfall of Haiyan.

3.1. Ensemble simulation of Typhoon Haiyan

The evolution and intensity of Haiyan were well represented in all of the experiments. In WEP, the simulated typhoon tracks around Leyte and Samar Island (125°E), on which Tacloban is situated, were distributed in a band about 500 km wide in the north–south direction with its centre north of Tacloban, and the tracks in the downscaled results followed the WEP result (figure 2). The westward-propagation speed of the typhoon was slower in the model simulations, compared with the observed track (best-track data), and as a result, the TC landfall on Leyte and Samar Island was delayed about a half day in the model simulations. This phase lag of the simulated typhoon is not attributable to the WRF model itself, because when Typhoon Haiyan was simulated with the same WRF 3 km/1 km model system with reanalysis data as the parent model, the result showed no phase lag (Mori et al 2014). This result suggests that the phase lag in this ensemble downscaling experiment is attributable to the performance of WEP. Because of this phase lag, however, we compared the evolution of the typhoon in relation to its longitudinal position (figure 3). The comparison indicates that in the ALL simulations both central surface pressure and maximum wind speed agree well with the best-track data, though the intensification of the TC that occurred before its landfall on Leyte and Samar Island was slightly delayed. The minimum central pressure (MCP) of the simulated typhoon was as high as 906 hPa (in experiment m02 (#1001)), which corresponds well to the intensity of the actual TC of 895 hPa as estimated by the Dvorak method (Dvorak 1975). The rapid decline of the typhoon strength after its landfall on Leyte and Samar Island was also reproduced by all of the ensemble simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To evaluate how well the experiments could reproduce Typhoon Haiyan, we used the results of experiment m02 (#1001), in which the simulated track was nearest to the observed track (best track). The structure of the eye-wall clouds was represented by cascade downscaling to a 1 km grid scale. Because the inner core convective band just around the eye is the engine of a typhoon, it is necessary to represent the structure of the eye-wall clouds for the sufficient intensification of a typhoon. Compared with weaker TCs, however, the radius of the eye of category 4 or 5 typhoons becomes smaller (about 20 km) (Weatherford and Gray 1988). Therefore, to represent the structure of the eye-wall clouds, the model must have a horizontal resolution of less than 10 km, and use of a model with a 2 km or 1 km grid is better because the simulation can be performed without a cumulus convective scheme (Gentry and Lackmann 2010, Kanada et al 2012). Supplementary figure S1 depicts the convective rain bands within and around the TC core in the 20, 5 and 1 km models. In the coarser resolution models, the eye-wall structure appears smoothed. To simulate the intensity of the disturbance better, the structure within the TC's core must be sufficiently resolved. Comparison of our simulation results (experiment m02 (#1001)) with the radar reflectivity data for Typhoon Haiyan obtained by Guiuan station (figure 4) showed that WRF 1 km was able to reproduce the structure around the eye-wall clouds well.

Zoom In Zoom Out Reset image size

3.2. Difference in the development of the Typhoon between NAT and ALL conditions

Next, we compared the ALL results with the NAT results obtained by the cascade downscaling method and found that the influence of the present climate state in comparison to the hypothetical natural condition on the intensity of the worst case storm was robustly indicated. In 15 of the 16 WRF experiments, the MCP was lower in the ALL simulations (mean difference and standard deviation; −6.44 and 4.98 hPa) than in the NAT simulations (figure 5(a)). In addition, the maximum surface wind speed was stronger in 15 of 16 ALL simulations (mean difference and standard deviation; 2.89 m s⁻¹ and 2.06 m s⁻¹) (figure 5(b)). This intensification has a 1% significant level.

Zoom In Zoom Out Reset image size

To investigate whether the intensity of the worst case storm was accounted for by the differences of climate states between the real climate in November 2013 and the hypothetical natural climate, we used the maximum potential intensity (MPI) following Emanuel (1986) to examine the differences in the potential development of the typhoon between NAT and ALL. We estimated the effect of changes in SST and atmospheric temperatures (ATM) on the MPI separately (supplementary table S1). The total contribution of SST and ATM to the MPI difference was −12 hPa (supplementary table S2). The contribution from SST changes alone ranged from −17 to −16 hPa, and that from ATM alone ranged from +3 to +4 hPa. The difference in the sign of the contributions from SST and ATM is consistent with the results of a previous study (Knutson and Tuleya 2004). Our result suggests that the SST change is the main contributor to the MPI difference between NAT and ALL. In the oceanic region around the Philippines, the SST was from +0.2 to +0.8 K higher in the ALL simulations compared to the NAT simulations (supplementary figure S2), and these higher SSTs can account for the MPI difference.

3.3. Storm surge estimation

Finally, we used the nonlinear shallow-water long-wave model SuWAT (Kim et al 2008) to estimate increases in the water surface elevation due to the worst case storm. We focused on the maximum water elevation, the storm surge height at Tacloban, where Typhoon Haiyan caused catastrophic damage. To estimate a disaster at such a specific point, the largest permissible bias in the track compared to the best track is 50 km (empirically estimated from the maximum wind speed radius of the typhoon). Because only ten ensemble members satisfied this condition for each of ALL and NAT experiments, we applied SuWAT to estimate the surge at Tacloban to only 2 × 10 experiments. The average maximum surge height in the ten D3 experiments was 2.60 m (standard deviation 1.36 m) and 2.19 m (standard deviation 1.00 m) in ALL and NAT, respectively, an amplification of around 20% in ALL (supplementary table S3). The maximum surge height along the coast, simulated in experiment m02 (#1001), was 4.27 m in ALL, whereas in NAT it decreased to 3.80 m (figure 6, shown with land inundation survey data (Mori et al 2014, Tajima et al 2014) which influenced local bathymetry and wind waves). Although there is not real surge data, the reanalysis of maximum surge height shows 5.15 m (Mori et al 2014). While the surge height in ALL is underestimated in comparison with the reanalysis, the spatial pattern of maximum surge height is similar to reanalysis by Mori et al (2014). As the first seiche mode in Leyte Gulf increased the amplitude of the storm surge at Tacloban, the surge height at any specific point is very sensitive to the position of the TC track. We estimated the mean maximum surge height using the ALL surface wind amplification data (5.82%; the difference between NAT and ALL) with the NAT track data (ALL_Pseudo). The results showed that the difference in the track between the ALL and NAT simulations had limited influence on the mean maximum surge height. Thus, the difference in maximum surge height between them is explained mainly by the difference in the cyclone intensity. The difference of maximum surge height between NAT and ALL is 0.47 m and amplification of surge height is much larger than surface wind amplification due to nonlinear characteristics of momentum transfer between atmospheric and ocean interface.

Zoom In Zoom Out Reset image size

There are many studies of historical sea level rise (IPCC-AR5, WGI, chapter 13). For example, Church and White (2011) investigate sea level change from the late 19th to the early 21st century. Although there were only a few tide gauges that operated in the middle 19th century, their papers figure 7 shows some sea level change records estimation from several data sources. From their figure, c.a. 0.25–0.30 m sea level rise is found and other studies also show similar value around this area if we exclude local land subsidence. Wind driven circulation or ENSO signal may also produce sea level change locally, but as shown in Woodworth et al (2008), the time scale caused by such ocean circulation is much shorter than 150 years. The IPCC Fifth Assessment Report also discusses the vulnerability of coastal regions to other physical processes, such as storm surges at extreme sea level. Changes in the severity of storm surges are some of the examples of how climate change affects coastal regions, but it is difficult to make an impact assessment at particular regions quantitatively. In our paper, the storm surge height increases at most 0.50 m between NAT and ALL experiments, caused by the intensification of the typhoon which is larger than historical sea level rise (estimated as 0.25–0.30 m).

We also examined the influence of SLR of 0.30 m, on the storm surge level with higher mean sea level condition. However, the maximum sea surface level considering a higher mean sea level condition didn't give significant changes in the storm surge height. Theoretically, the storm surge height is inversely proportional to the depth of the bay. Since the water depth near Tacloban in Leyte Gulf has depth of approximately 10 m, 0.30 m mean sea level rise almost cancels out the change of total surge height, although there remain some dependencies on characteristics of the bay and cyclones. The combination of a storm surge and sea level rise is nonlinear, thus it is important to know how sea level rise and dynamic effects, wind and pressure surges, will change in the future.

Typhoon Haiyan (local name Yolanda), the most catastrophic tropical cyclone ever to land in the western North Pacific Ocean, struck the Philippines on 8 November 2013. The typhoon and especially the storm surge in the Leyte Gulf that accompanied it killed more than 6000 people in Tacloban (Schiermeier 2013). We conducted ensemble simulations with very high resolution regional climate models and a surge model, and reproduced well the pressure depression, wind speed and surge level of Typhoon Haiyan, as an example of a worst case scenario. Furthermore, we compared these results with the results of ensemble simulations of a hypothetical natural event, one without human influences, and found that the simulated worst case typhoon and the accompanying storm surge in the real condition became worse than those in the hypothetical natural climate without anthropogenic forcing. In 15 of 16 ensemble simulations, the typhoon became stronger than it did in the hypothetical natural cases, and the height of the storm surge around Tacloban increased by around 20%.

When we made the NAT external condition, we omitted anthropogenic warming in air temperature difference as well as that in SST. It is well known that vertical wind shear and mid-level entropy deficit are important in determining the intensity of tropical cyclone (Gray (1975), Tang and Emanuel (2010, 2012)). Although previous studies have suggested the weakening of Pacific Walker circulation (reductions of wind shear in the western Pacific region), and anthropogenic contributions to this observed trend, there are substantial uncertainties in observations and biases in model simulations (Tokinaga et al 2012, DiNezio et al 2013). Supplementary figure S3 shows the histograms of vertical wind shear in the 100-member ensembles of ALL-gcm and NAT-gcm runs of MIROC5-AGCM (Shiogama et al 2014). Although anthropogenic factors weaken the vertical wind shear in the ensemble mean (which would increase the Typhoon intensity (Tang and Emanuel, 2010, 2012)), the variances due to internal variability are large. Because of the poor signal-to-noise ratio, we have not subtracted these differences in wind shear from the NAT conditions of RCM runs. It might cause our estimates of anthropogenic influence on the Typhoon intensity to become more conservative.

We should also note here, that upper-ocean mixing is not accounted for because we do not have a sound assumption for prescribing the upper-ocean structure in natural climate conditions. Vincent et al (2014) insisted that the mixing has a negative feedback on TC intensity for very strong storms, but the influence of such a mechanism is not accounted in our study.

As we introduced in section 1, to assess the disaster risk reduction caused by the climate change, we had quite a different approach, rather than considering the disaster mitigation as we discussed in the paper. This is to consider the disaster prevention. For this purpose, we should estimate changes in the frequency of event occurrence. However, Hansen et al (2014) suggests that it is difficult to estimate the change of likelihood with a high confidence level when the sample numbers are not high enough, especially for extreme events. To increase the sample number efficiently in cases where we are not able to conduct sufficient ensemble experiments, we have another option of introducing a stochastic hurricane model of the type used in Lin et al (2012). This method imposes a stronger limitation of the assumed hurricane structure and intensity than a dynamically downscaling approach. Thus, it is indispensable to carry out both types of approach jointly, to clarify the disaster risk of extreme events as a whole.

There are uncertainties associated with estimations of the effects of anthropogenic signals on SST (Christidis and Stott 2014, Shiogama et al 2014) and attribution analysis results may also be sensitive to models used. We showed here that very high resolution regional climate models are necessary to simulate a category 5 typhoon. It should be kept in mind that more accurate representation of the typhoon with higher spatial resolution does not necessary ensure more reliable estimates of the anthropogenic contribution. Inter-comparison studies of very high resolution regional climate models should be carried out in which different estimates of anthropogenic signals are used.

This work was supported by the SOUSEI programmes of the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. The Guiuan radar data were supplied by PAGASA, the Philippines, and Dr Kubota of Japan Agency for Marine-Earth Science and Technology (JAMSTEC) helped us to read the original data. We also acknowledge Dr Fujita of the Numerical Prediction Division of the Japan Meteorological Agency (JMA) and Dr Yamaguchi of the Meteorological Research Institute for providing us with information about JMA's operational weekly ensemble prediction system.