Cyclone generation Algorithm including a THERmodynamic module for Integrated National damage Assessment (CATHERINA 1.0) compatible with Coupled Model Intercomparison Project (CMIP) climate data

Tropical cyclones are responsible for a large share of global damage resulting from natural disasters, and estimating cyclone-related damage at a national level is a challenge attracting growing interest in the context of climate change. The global climate models, whose outputs are available from the Coupled Model Intercomparison Project (CMIP), do not resolve tropical cyclones. The Cyclone generation Algorithm including a THERmodynamic module for Integrated National damage Assessment (CATHERINA), presented in this paper, couples statistical and thermodynamic relationships to generate synthetic tracks sensitive to local climate conditions and estimates the damage induced by tropical cyclones at a national level. The framework is designed to be compatible with the data from CMIP models offering a reliable solution to resolve tropical cyclones in climate projections. We illustrate this by producing damage projections in representative concentration pathways (RCPs) at the global level and for individual countries. The algorithm contains a module to correct biases in climate models based on the distributions of the climate variables in the reanalyses. This model was primary developed to provide the economic and financial community with reliable signals allowing for a better quantification of physical risks in the long term, to estimate, for example, the impact on sovereign debt.

1 Introduction

Climate-related physical risks pose a growing threat to humanity, and the design and implementation of adequate adaptation and mitigation measures require assessment of future physical risks at a national and global scale. The projections of global climate models are an important source of information about the future climate. However, the spatial resolution of these models is unfortunately still not sufficient to fully resolve extreme events, particularly tropical cyclones. On the other side of the spectrum, integrated assessment models (IAMs) directly assess the impact of the future climate on economic activity. Although these models are used to calibrate optimal mitigation or adaptation pathways, most of them only embed a very simple physical damage module usually limited to a generic temperature damage function. The objective of this paper is to fill the gap between climate models and integrated assessments: we build synthetic cyclones based on the climate data produced by the global climate models and evaluate the economic damage of these synthetic cyclones under various assumptions regarding the socioeconomic scenarios. Tools to assess the impact of future cyclones in shared socioeconomic pathways are starting to appear in the literature. For example, Geiger et al. (2021) evaluate the population exposure. Our study instead focuses on the damage costs of tropical cyclones with the aim of including these advanced signals in integrated economic modeling.

The physics of tropical cyclones has been thoroughly studied in the literature. The thermodynamic cyclone theory builds upon the seminal contributions by Emanuel (1988) and followed by Holland (1997) and Emanuel (1999). Concerning the impact of climate change, it is well known that the presence of greenhouse gases in the atmosphere increases the radiative forcing, which leads to a progressive warming of the atmosphere (Butchart et al., 2000) and the rise in sea surface temperatures (Solomon et al., 2007; Pachauri et al., 2014). This phenomenon increases the amount of energy available for cyclones to grow in intensity (Emanuel, 1991), and this growth is already starting to be measurable (Emanuel, 2005). A 2010 review concluded that it is unclear whether past changes in tropical cyclone (TC) activity have exceeded the variability expected from natural causes (Knutson et al., 2010), while a more recent review (Walsh et al., 2016) confirms that studies are beginning to suggest climate trends of various kinds in TC data over the past few decades. This review states that no trend in average cyclone frequency or intensity has been detected globally, but a substantial increase in the proportion of intense cyclones has been demonstrated, both globally and in individual basins, with the exception of the North Pacific.

Models relying on statistical relationships (James and Mason, 2005; DeMaria and Kaplan, 1994; Kaplan and DeMaria, 1995) available in the literature produce synthetic cyclones with properties closely resembling those of observed cyclones. Recently, Bloemendaal et al. (2020) developed a modeling framework to simulate realistic synthetic tropical cyclone tracks: the Synthetic Tropical cyclOnes geneRation Model (STORM). This model computes the maximum pressure intensity (MPI) associated with the sea surface temperature (SST) and uses this potential as a predictor in the central pressure dynamics (James and Mason, 2005). In line with Merrill (1987), we find that, although the SST plays a major role, this variable alone is not a reliable predictor of whether a given storm will intensify. Thus, we prefer to rely on the formulation of Holland (1997) and to model the effect of climate change on the maximum potentials in the different scenarios through a better description of the underlying thermodynamic phenomenon, well described by Emanuel (1988), Holland (1997) or Emanuel (1999). We therefore develop an alternative to STORM by adding a thermodynamic module from the perspective of producing cyclones in different climate scenarios. In particular, we retrieve two additional variables (relative humidity and tropopause temperature) from climate models to bridge the gap between data from global circulation models (GCMs) and the theory of the intensification of tropical depressions.

Risk assessments have been developed based on hurricane potential intensity maps to assess the damage in the US and around the world (Emanuel et al., 2008; Emanuel, 2011; Mendelsohn et al., 2012). On the damage modeling side, Bresch (2017), Lüthi (2019), and Aznar Siguan and Bresch (2019) set up a platform for physical risk estimation (CLIMADA), coupled with a database of estimated values of local assets (Eberenz et al., 2020a, b). The asset resolution (30 arcsec) and geospatial description of extreme events are particularly advanced. We propose a simplification of this work that applies in the context of national level assessment. The CLIMADA framework focuses on damage modeling based on global aggregated temperature projections and does not make use of the climate data produced by atmosphere–ocean general circulation models (AOGCMs). Coupling CLIMADA and STORM methodology with an extended thermodynamic module fitted on four climate variables, our approach provides a novel long-term tail risk assessment at a national level, offering an adaptive framework to estimate investments required to mitigate this risk.

This paper makes three contributions. First, we provide an algorithm to generate synthetic cyclones from climate data inspired by Bloemendaal et al. (2020) fitted on four physical variables extracted from ERA5 reanalysis and including a thermodynamic module to better describe cyclone physics. Second, we build an algorithm generating synthetic tracks directly from CMIP models, expose the biases in CMIP5 datasets and propose a correction module based on Vrac et al. (2012). Third, we bridge the gap between climate data and damage modeling by using the physical asset values from Eberenz et al. (2020b) and computing the damage along cyclone tracks using the region-specific damage functions designed in the CLIMADA project. Combining open data sources and methodologies allows us to propose a complete bottom-up integrated physical risk assessment model for tropical cyclones presented in Fig. 1: the Cyclone generation Algorithm including a THERmodynamic module for Integrated National damage Assessment (CATHERINA).

Figure 1 Cyclone generation Algorithm including a THERmodynamic module for Integrated National damage Assessment (CATHERINA) framework.

The process is the following. The cyclones are initiated with spatial and seasonal distribution estimated on the International Best Track Archive for Climate Stewardship (IBTrACS) database, similarly to Bloemendaal et al. (2020). Their movement is described with a simple autoregressive stochastic process (James and Mason, 2005; Bloemendaal et al., 2020). Along the simulated cyclone tracks, we retrieve climate variables from climate models and compute locally the maximum potential intensity based on the simplified expression in Holland (1997). Some controls such as the maximum pressure drop observed for the corresponding temperature and the decay relationship for cyclones evolving over land are also fitted for each basin and applied in the synthetic track-generation algorithm. Extracting the climate variables from different models allows us to correct the biases and evaluate model uncertainty. Next, we use the physical asset values (Eberenz et al., 2020b) and regional damage functions from CLIMADA to evaluate the cyclone-related damage at a national level. This step requires us to extract the local physical asset values and aggregate them on tiles of length defined proportionally to the average radius to the maximum wind (approximately 50 km) along the cyclone path. Summing the losses along tracks for each year and for each country allows us to establish a national assessment of the damage generated by tropical cyclones under various assumptions concerning future adaptation to cyclone risk.

CATHERINA aims to provide country level estimates for future damage from tropical cyclones, consistent with climate model projections, with mainly economic and financial applications in mind. Examples of applications include estimating the impact of cyclones on the creditworthiness of government debt, providing a physical risk module for integrated assessment models, creating physical risk scenarios for stress testing the resiliency of the financial system at a country level and at a global level, etc. Given this aggregate country level analysis objective, our model is certainly a simplification compared to state-of-the-art cyclone dynamics models, does not aim for precise prediction of individual cyclone tracks and does not integrate a bottom-up description of damage to individual assets. At the same time, our paper improves earlier studies of cyclone damage on the aggregate level. For example, compared to Mendelsohn et al. (2012), our model uses several climate scenarios with state-of-the-art bias correction as well as shared socioeconomic pathways (SSPs) to project population and regional domestic products and is based on precise physical asset value distribution.

The paper is organized as follows. Section 2 first describes the datasets used to fit the model on ERA5 and to generate synthetic cyclones based on both the ERA5 and CMIP5 model datasets (Sect. 2.1). Next, we present the statistical calibration process and the details of the thermodynamic instrumental variables (Sect. 2.2 to 2.4). Section 3 presents the dataset used to represent physical asset exposure and the SSP framework used to project their future value. Section 4 recalls the calibration methods implemented in the CLIMADA environment to fit the regional damage functions, defines the parameters of these functions using Eberenz et al. (2021) and applies them along the synthetic tracks to study the distribution of national annual damages. Section 5 explores the properties of the synthetic tracks produced with the ERA5 and seven CMIP models, assesses climate uncertainty and introduces the bias-correction module in the context of changing climate conditions. To close this section, we present the global and regional projections of cyclone damage between 2070 and 2100 obtained with CATHERINA and illustrate briefly how these damage estimates may be used to quantify the future impact of tropical cyclones on sovereign default risk.

2 Modeling tropical cyclones

Our model structure follows Bloemendaal et al. (2020) with three main modeling steps: genesis, displacement of the eye and calibration of the cyclone properties. The STORM model relies on statistical relationships (James and Mason, 2005; DeMaria and Kaplan, 1994; Kaplan and DeMaria, 1995). This simulation method differs from the purely thermodynamic approach developed by Kerry Emmanuel (Emanuel, 1999; Emanuel et al., 2008).

The major change in our specification compared to Bloemendaal et al. (2020) is that we use the local definition of the available thermodynamic potential based on climate data. In particular, we use relative humidity (RH) and tropopause temperature (at 50 hPa, T_tropo) for better theoretical representation of the physics underlying the intensification process. In this section, we present the process of generating synthetic tracks, characterized by the maximum wind (V_t) and central pressure ($P_t^{\mathrm{c}}$) at each time step, given the climate conditions extracted from climate models.

2.1 Input data

In the CATHERINA framework, we fit the properties of historical cyclones (IBTrACS database) on past climate reanalyses (ERA5) from the perspective of describing future cyclones based on global climate model outputs (CMIP) with a lower spatial and temporal resolution. This perspective constrains us to use monthly data.

2.1.1 Climate model data (CMIP)

CATHERINA aims at generating cyclone tracks with properties drawn from climate models to enable national damage assessments, bridging the gap between AOGCM outputs and damage assessments. To reduce the bias in the variables produced by climate models and to evaluate the performance of CATHERINA on past data by comparing the simulated cyclone damages to the realized ones, we use historical simulations (as opposed to future climate projections) from the Coupled Model Inter-comparison Project (Phase 5) models (Taylor et al., 2012). We use the historical climate simulations at the monthly frequency for the RH at 2 m, SST, mean sea level pressure (MSLP) and tropopause temperature (T_tropo) (at a pressure level of 50 hPa) from the Copernicus climate data store.¹ CMIP5 data are used in the fifth assessment report of the Intergovernmental Panel on Climate Change (IPCC). The latest synthesis report in 2022 (IPCC AR6) uses CMIP6 datasets, but in the present paper, we use CMIP5 data because of the broader availability of climate variables.

We use models from the following climate centers: NASA and the Goddard Institute for Space Studies (GISS-E2-H, USA), Institut Pierre Simon Laplace (IPSL-CM5A-NR, France), Bureau of Meteorology – Commonwealth Scientific and Industrial Research Organisation (ACCESS1-0, BoM-CSIRO, Australia), Beijing Climate Center (bcc-csm1-1-m, China), Institute of Numerical Mathematics (INMCM4, Russia), Norwegian Climate Centre (NorESM1-ME, Norway), and Canadian Centre for Climate Modelling and Analysis (CanESM2, Canada). The spatial resolution goes from 0.75 to 2.5^∘ depending on the model (see Table 1). Each climate model produces a potentially biased estimate of multiple climate variables at the spatial resolution given in Table 1 and on a monthly basis. The choice of the GCMs was driven by the availability of the variables of interest in the CDS in the representative concentration pathways used in the exercise (RCPs 2.6, 4.5 and 8.5 W m⁻²) in both single level and multiple pressure level monthly data in the same ensemble (r1i1p1). We also aimed at having multiple regions represented.

To reduce the influence of model bias, we use a large number of models and consider the distribution of results provided by all the models. Then we correct, variable by variable and for each basin and each model, the biases with respect to the reanalysis along the same tracks.

Table 1 Climate data resolution.

2.1.2 ERA5 reanalysis

Climate reanalyses describe the historical climate conditions, obtained by assimilating all available observations into the models. They provide numerical estimates of atmospheric parameters (e.g., air temperature, pressure and wind) at different altitudes/pressure levels and surface parameters (such as rainfall, soil moisture content, ocean-wave height and sea surface temperature) on a single level. We use reanalyses to calibrate the cyclone generation algorithm based on the most realistic available estimates of climate variables.

We use the ERA5 reanalysis from the European Centre for Medium-Range Weather Forecasts (ECMWF) to fit the CATHERINA model² (Hersbach et al., 2020). This dataset covers the Earth on a 30 km grid (∼0.25^∘) and resolves the atmosphere using 137 levels from the surface up to a height of 80 km. In this paper, to ensure compatibility with CMIP5 models, we extract MSLP, SST, sea level RH and tropopause temperature (T_tropo) at the monthly frequency. Because ERA5 better resolves past tropical cyclones than climate models, the historical mean sea level pressure values in ERA5 are influenced by their presence. Consequently, we retrieve the mean sea level pressure and environmental relative humidity 500 km ($\sim \mathrm{5}{}^{\circ }$ longitude) away from the storm center to extract a value for P_env, which is meant to represent the pressure – at a given latitude and season – in normal environmental conditions.³

2.1.3 Historical cyclone tracks (IBTrACS)

We use the IBTrACS database (Knapp et al., 2010).⁴ This database provides information on past cyclone tracks at 3 h frequency. We remove the events classified as disturbance or extratropical and do not consider the South Atlantic basin (see Fig. A1 for more information). Climate reanalysis availability requires us to focus on post-1980 cyclones, which reduces the database to 4574 cyclones. In the context of an integrated damage assessment, to focus on the events that have a potentially substantial impact on assets, we select only tropical cyclones with maximum wind speed exceeding 35 m s⁻¹, obtaining 2966 on the full database and 1451 focusing on tropical cyclones between 1980 and April 2020. In Fig. 2, we plot the central pressure along each cyclone life. This graph suggests that the cyclone phases are fully represented in the database, i.e., from the genesis to dissipation. The northern Indian basin has the lowest number of reported events with wind speed above 35 m s⁻¹ (50 compared to 291 for the eastern Pacific, 185 for the North Atlantic, 305 for southern India, 158 for the South Pacific and 515 for the western Pacific) with variable reporting quality, which explains the more erratic shapes of the central pressure. For example, the return to normal of some events does not seem to be completely reported for this basin, as indicated in Fig. 2. We extract the maximum wind speed, cyclone eye pressure and coordinate variations of the eye from this database.

Figure 2 Normalized evolution of the central pressure (hPa) during a cyclone's life. The central pressure values are retrieved from IBTrACS. For the western Pacific we reduced the sample to events between 1995 and 2020 for visualization purposes.

2.2 Cyclone genesis

The scientific consensus is that climate change will induce a reduction in tropical cyclone frequency: “Existing modeling studies also consistently project decreases in the globally averaged frequency of tropical cyclones, by 6 to 34 %. Balanced against this, higher resolution modeling studies typically project substantial increases in the frequency of the most intense cyclones” (Knutson et al., 2020, p. 1). Although thermodynamic descriptions of cyclone genesis exist in the literature (Gray, 1975; DeMaria et al., 2001), we choose to rely on a simple statistical model based on past frequencies for the genesis.

The number of synthetic cyclones each year is determined by the Poisson distribution in each basin, with the intensity parameter defined as the average number of cyclones per year in the historical data. We use the parameters given in Bloemendaal et al. (2020), i.e., ${\widehat{\mathit{\lambda }}}_{\mathrm{EP}}=\mathrm{14.5}$ for the eastern Pacific, ${\widehat{\mathit{\lambda }}}_{\mathrm{NA}}=\mathrm{10.8}$ for the North Atlantic, ${\widehat{\mathit{\lambda }}}_{\mathrm{NI}}=\mathrm{2.0}$ for northern India, ${\widehat{\mathit{\lambda }}}_{\mathrm{SI}}=\mathrm{12.3}$ for southern India, ${\widehat{\mathit{\lambda }}}_{\mathrm{SP}}=\mathrm{9.3}$ for the South Pacific and ${\widehat{\mathit{\lambda }}}_{\mathrm{WP}}=\mathrm{22.5}$ for the western Pacific. The parameters would have been smaller if estimated using our filtered database of tropical cyclones with wind speeds above 35 m s⁻¹: ${\widehat{\mathit{\lambda }}}_{\mathrm{EP}}^{\mathrm{35}}=\mathrm{7.31}$, ${\widehat{\mathit{\lambda }}}_{\mathrm{NA}}^{\mathrm{35}}=\mathrm{4.43}$, ${\widehat{\mathit{\lambda }}}_{\mathrm{NI}}^{\mathrm{35}}=\mathrm{1.6}$, ${\widehat{\mathit{\lambda }}}_{\mathrm{SI}}^{\mathrm{35}}=\mathrm{6.81}$, ${\widehat{\mathit{\lambda }}}_{\mathrm{SP}}^{\mathrm{35}}=\mathrm{4.00}$ and ${\widehat{\mathit{\lambda }}}_{\mathrm{WP}}^{\mathrm{35}}=\mathrm{11.86}$. However, we maintain the original parameters used in STORM to take into consideration the fact that some events will be generated in conditions not favorable for the development of cyclones and cleared out of the database. The number of cyclones making landfall being critical for damage estimation (Hall and Jewson, 2007; Lee et al., 2018; Arthur, 2021), we ensured that the simulated landfall rates are in line with the observations over the historical period: Fig. 3 shows that the parameters of Bloemendaal et al. (2020) lead to approximately the same number of intense cyclones making landfall per year in each basin in the historical data, in our simulations based on reanalysis data and in simulations based on GCMs. We note however that the framework can be further improved by choosing the intensity parameter to match exactly the average historical number of cyclones making landfall.

Figure 3 Number of tropical cyclones over 35 m s⁻¹ making landfall simulated using Bloemendaal et al. (2020) parameters over the historical period (1980–2010).

Similarly, the temporal and spatial initial positions of synthetic future cyclones (longitude x, latitude y and starting month m) are generated by independent sampling from the historical distribution of these variables. Figure 4 shows the geographical distribution of cyclones retrieved from IBTrACS (i.e., x_obs,y_obs), and the histograms in Fig. 5 show the monthly distribution (m_obs) of cyclones in each basin.

Figure 4 Spatial distribution of genesis points (tropical cyclones over 35 m s⁻¹ in IBTrACS). The color scale corresponds to the count of cyclones per 5×5^∘ box (truncated to 25 for scaling reasons). The genesis location corresponds to the first reported point of each track.

Figure 5 Relative monthly distribution (%) of cyclones in each basin defined by potentially damaging cyclones (over 35 m s⁻¹) in IBTrACS. Each bar gives the probability of each cyclone being allocated to a given month.

2.3 Cyclone trajectories

A rich literature focuses on cyclone tracking algorithms; see, e.g., Neu et al. (2013) for a comprehensive review. Although more advanced definitions have been proposed (Hall and Jewson, 2007; Fabregat et al., 2016), we choose, in line with Bloemendaal et al. (2020), to implement a simple autoregressive model for cyclone coordinates. Following James and Mason (2005), the time evolution of the latitude and longitude of the cyclone center is described by the following stochastic dynamics:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{x}_t=a_{\mathrm{0}}+a_{\mathrm{1}}\mathrm{\Delta }{x}_{t-\mathrm{1}}+{\mathit{\epsilon }}_t^x,& {\displaystyle }{\mathit{\epsilon }}_t^x\sim \mathcal{N}(\mathrm{0},{\mathit{\sigma }}_x),\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{y}_t=b_{\mathrm{0}}+b_{\mathrm{1}}\mathrm{\Delta }{y}_{t-\mathrm{1}}+{\displaystyle \frac{b_{\mathrm{2}}}{y_t}}+{\mathit{\epsilon }}_t^y,& {\displaystyle }{\mathit{\epsilon }}_t^y\sim \mathcal{N}(\mathrm{0},{\mathit{\sigma }}_y).\end{array}$$

Here x_t and y_t are the latitude and longitude of the cyclone center sampled with a 3 h time step. $\mathrm{\Delta }{x}_t=x_t-x_{t-\mathrm{1}}$, $\mathrm{\Delta }{y}_t=y_t-y_{t-\mathrm{1}}$, ${\mathit{\epsilon }}_t^x$ and ${\mathit{\epsilon }}_t^y$ are i.i.d. (independent and identically distributed) noises independent of one another, and the constants a₀, a₁, b₀, b₁, b₂, σ_x and σ_y are fitted on IBTrACS data independently for each basin by least squares regression. The nonlinear term in the incremental variation of the latitude is justified by the tendency for cyclones to move away from the Equator, especially at very low latitudes (James and Mason, 2005, p. 183).

To take into account the dependency of the cyclone displacement on the location of the eye, and following Bloemendaal et al. (2020), we fitted these relationships locally using an additional grouping by 5^∘ longitude and latitude sections and months. Figure A2 illustrates the parameters a₀ and a₁, the adjusted R² and the number of observations averaged over months used to fit Eq. (1). We note that the trajectories in the northern Indian basin are less well captured in some areas, particularly longitudinal movements near the coast from Yemen, Oman, and the United Arab Emirates (see the maps of adjusted R² in Fig. A2).

This statistical definition of cyclone trajectories does not consider changes in track behavior. For example, observed trends in tropical cyclone translation speed (Kossin, 2018) and poleward migration of maximum intensity (Kossin et al., 2014) could be considered to improve the projections of tropical activity. This indeed has potential implications for tropical-cyclone-related risk in some areas where vulnerabilities are high and the present-day frequency of tropical cyclones is low (Bruyère et al., 2019).

2.4 Thermodynamic description of cyclone intensity

The intensity of cyclones in our model is defined through the following five steps described in detail in subsequent paragraphs. The wind–pressure relationship (WPR) links the central pressure to the maximum 10 min sustained wind speed. The local thermodynamic MPI is determined from local meteorological variables. The maximum pressure drop (MPD) is determined from historically observed pressures. The depression dynamics (DD) along tracks is defined using an autoregressive stochastic equation. When the cyclone arrives on land, a statistical decay relationship (SDR) dictates the evolution.

2.4.1 WPRs

We describe the cyclone intensity through its central pressure $P_t^{\mathrm{c}}$, which is linked to the maximum wind through an empirical relationship. Let V_t be the maximum 10 min sustained wind speed (m s⁻¹)⁵ of the cyclone at time t. This maximum wind is observed around 50 km away from the storm center on average⁶ and reported in the IBTrACS dataset for historically observed cyclones.

The WPR is calibrated separately for each basin and takes the following form:

$$\begin{array}{}\text{(3)}& V_t=a{\left(P_{\mathrm{env}}(x_t,y_t,t)-P_t^{\mathrm{c}}\right)}^b,\end{array}$$

where $P_{\mathrm{env}}(x,y,t)$ is the MSLP extracted 500 km away from the eye location at this time in ERA5 and $P_t^{\mathrm{c}}$ is the central pressure extracted from IBTrACS. This relationship is illustrated in Fig. 6, and the parameters a and b are fitted on ERA5 and IBTrACS data using nonlinear least squares.

Figure 6 Maximum wind (m s⁻¹) and pressure drop (hPa) values from IBTrACS together with the fit of Eq. (3). The coefficients a and b fitted per basin are provided in Table A2.

We acknowledge that most cyclone track data use WPRs to determine P_c. The conversion back to wind speed from the reported P_c using a basin-specific WPR still introduces errors, as different WPRs are used to operationally estimate P_c within each basin (Harper, 2002; Courtney, 2009; Courtney and Burton, 2018; Courtney et al., 2021). However, given the similarity of the relationships, we find that basin level estimations are a sufficient proxy in the context of this illustration of the CATHERINA framework.

2.4.2 Local thermodynamic MPI

We compute the local thermodynamic MPI following the thermodynamic relationships defined in Holland (1997). This is particularly relevant in the context of a national damage assessment under a changing climate. Greenhouse gas emissions not only warm up the oceans, but also cool down the lower stratosphere (Butchart et al., 2000; Forster et al., 2007; Ramaswamy et al., 2006). Thus, the tropopause temperature – that is, the temperature corresponding to a pressure of 50 hPa or to an altitude of approximately 20 km available in a multi-level CMIP dataset (T_tropo) – must be included in the modeling of the intensification process. Indeed, the thermodynamic efficiency factor ℰ_t proportional to the difference between the tropopause and sea surface temperatures plays an essential role in the determination of the central pressure of tropical cyclones. The relative humidity (RH) (which changes with climate change; see Sherwood et al., 2010) is also an influential parameter allowing for a better description of thermodynamic potential enabling cyclone intensification. Adding these two climate variables enables the CATHERINA model to better take into account the additional energy potential due to the widening of temperature differences between the sea surface and upper troposphere and variation in moist entropy.

Following the seminal formulation in Emanuel (1988) and integrating additional simplifications proposed in the subsequent paper (Emanuel, 1991) leads to the following framework summarized in Holland (1997).

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }{\mathrm{MPI}}_t=\mathrm{MSLP}(x_t,y_t,t)\cdot {\exp }^{-X_t},\text{(5)}& {\displaystyle }& {\displaystyle }{X}_t={\displaystyle \frac{{\mathcal{E}}_t\cdot \mathrm{SST}(x_t,y_t,t)\cdot \mathrm{\Delta }{S}_t^{\mathrm{m}}-\frac{f(y_t{)}^{\mathrm{2}}{r}_{\mathrm{env}}^{\mathrm{2}}}{\mathrm{4}}}{R_{\mathrm{d}}\cdot \mathrm{SST}(x_t,y_t,t)}},\text{(6)}& {\displaystyle }& {\displaystyle }{\mathcal{E}}_t={\displaystyle \frac{\mathrm{SST}(x_t,y_t,t)-T_{\mathrm{tropo}}(x_t,y_t,t)}{\mathrm{SST}(x_t,y_t,t)}},\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{S}_t^{\mathrm{m}}=R_{\mathrm{d}}\ln \left({\displaystyle \frac{\mathrm{MSLP}(x_t,y_t,t)}{P_{t-\mathrm{1}}^{\mathrm{c}}}}\right)+{\displaystyle \frac{L_{\mathit{\upsilon }}({q_{\mathrm{c}}^{\star }}_t-q_t^{\mathrm{env}})}{\mathrm{SST}(x_t,y_t,t)}},\end{array}$$

with moist entropy potential defined along track using specific humidity in the eye vs. at environmental conditions:

$$\begin{array}{ll}\text{(8)}& {\displaystyle }& {\displaystyle }{q_{\mathrm{c}}^{\star }}_t={\displaystyle \frac{\mathrm{3.08}}{P_{t-\mathrm{1}}^{\mathrm{c}}}}\exp \left({\displaystyle \frac{\left(\mathrm{SST}(x_t,y_t,t)-\mathrm{273.15}\right)}{\mathrm{SST}(x_t,y_t,t)-\mathrm{29.65}}}\right),{\displaystyle }& {\displaystyle }{q}_t^{\mathrm{env}}={\displaystyle \frac{\mathrm{3.08}\cdot \mathrm{RH}(x_t,y_t,t)}{\mathrm{MSLP}(x_t,y_t,t)}}\\ \text{(9)}& {\displaystyle }& {\displaystyle }\times \exp \left({\displaystyle \frac{\mathrm{17.67}\left(\mathrm{SST}(x_t,y_t,t)-\mathrm{273.15}\right)}{\mathrm{SST}(x_t,y_t,t)-\mathrm{29.65}}}\right).\end{array}$$

where $(x_t,y_t,t)$ are the coordinates of the eye defined in Eqs. (1) and (2), $\mathrm{SST}(x_t,y_t,t)$ and $T_{\mathrm{tropo}}(x_t,y_t,t)$ are, respectively, the sea surface and tropopause temperatures, R_d=287.058 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ is the specific gas constant for dry air, $\mathrm{MSLP}(x_t,y_t,t)$ is the mean local sea level pressure, $\mathrm{RH}(x_t,y_t,t)$ is the near-surface relative humidity at 2 m extracted from the monthly dataset of ERA5 climate reanalysis or CMIP climate models. f(y_t)=2ωsin (y_t) is a Coriolis parameter depending on the latitude, r_env is the distance between the eye and the area under regular conditions (fixed at 500 km), q_env and $q_{\mathrm{c}}^{*}$, respectively, are the specific humidity at environmental conditions and at saturation, i.e., for RH=100 %, in the eye. ΔS^m is the difference of moist entropy between the environment and the storm center and L_υ is the latent heat of vaporization. The distributions of the climate variables and the instrumental indicators computed from ERA5 climate variables along IBTrACS involved in this step are shown in Fig. 7.

Figure 7 Empirical distribution of physical variables and CATHERINA instrumental variables computed from ERA5 along tracks in IBTrACS.

2.4.3 MPD

Several papers, including Bloemendaal et al. (2020), link the sea surface temperature directly to the pressure drop, or equivalently the wind speed, via a statistical relationship (DeMaria and Kaplan, 1994). Merrill (1987) suggests that this predictor alone does not provide a good indication of whether a given storm will intensify. However, in line with Emanuel (1988) and DeMaria and Kaplan (1994), the sea surface temperature can be used to fix a limit for the pressure drop. Thus, in CATHERINA, to prevent the pressure drop from diverging in the projection, we cap it by the maximum observed pressure drop for the corresponding sea surface temperature:

$$P_t^{\mathrm{c}}:=max(P_t^{\mathrm{c}},\mathrm{MSLP}(x_t,y_t,t)-\mathrm{MPD}(\mathrm{SST}(x_t,y_t,t)\left)\right),$$

where the maximum pressure drop function is given by the following equation:

$$\begin{array}{}\text{(10)}& \mathrm{MPD}\left(\mathrm{SST}\right)=A+B\cdot e^{C\left(\mathrm{SST}(x_t,y_t,t)-T_{\mathrm{0}}\right)},\end{array}$$

with T₀ = 30 ^∘C. To fit this functional relationship, we first retrieve, for each basin, and for each value of the sea surface temperature SST, rounded to 0.1, the maximum observed value of the pressure drop in the basin. These values are shown as dots in Fig. 8. The coefficients A, B and C from the relationship (10) are then fitted to these values by nonlinear least squares. The resulting MPD function is shown in Fig. 8 as solid line.

The definition of the MPD is identical in STORM and CATHERINA, but the role of this quantity differs in the two models. Figure 8b gives a misleading idea about the strength of correlation between the sea surface temperature and the central pressure: indeed, fitting the distribution on the full sample (instead of just the maximum pressure drop for each temperature value) shows a much weaker influence of sea surface temperature alone (see Fig. 8a), even on a weekly basis (see Jien et al., 2017). However, this instrumental variable is essential for preventing CATHERINA from producing unrealistically low central pressure. On the other hand, limiting the maximum pressure drop in simulations to the parametric function given by Eq. (10) fitted by nonlinear least squares to the observed MPD leads to an excessively strong limitation on the intensity of the simulated cyclones, since many points in Fig. 8 (lower graph) are above the red curve. Therefore, we shift this parametric function upward to the highest observed point to relax this limitation.

Figure 8 (a) Full distribution of observed pressure drop (hPa) values for a given SST (^∘C). (b) Maximum observed pressure drop values for a given SST. The red line shows the least squares fit of Eq. (10), which corresponds to parameter values A=30.6, B=86.3 and C=0.19. The blue curve is the capping function used to prevent unrealistic pressure drops, obtained by shifting the red curve upwards.

2.4.4 Depression dynamics

The evolution of the central pressure depending on the intensification factor Y_t (which is defined differently in our approach and in STORM, as we explain below) is described by the following autoregressive stochastic equation (James and Mason, 2005):

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{P}_t^{\mathrm{c}}=c_{\mathrm{0}}+c_{\mathrm{1}}\mathrm{\Delta }{P}_{t-\mathrm{1}}^{\mathrm{c}}+c_{\mathrm{2}}{e}^{-c_{\mathrm{3}}\left[P_t^{\mathrm{c}}-Y_t\right]}+{\mathit{\epsilon }}_t^P,\text{(12)}& {\displaystyle }& {\displaystyle }{\mathit{\epsilon }}_t^P\sim \mathcal{N}(\mathrm{0},{\mathit{\sigma }}_{P^{\mathrm{c}}}).\end{array}$$

This relationship channels the effect of global warming, affecting the thermodynamic potential constructed with climate projections, on the cyclone intensity. Thus, the incremental variation of the central depression of the cyclone is linked to the difference between the central pressure at time t and the potential available in the environment.

The intensification module of CATHERINA is inspired by STORM (Bloemendaal et al., 2020). The main differences are the definition of the thermodynamic MPI used in Eq. (11) and the role played by the MPD. In Bloemendaal et al. (2020), the MPI is defined by subtracting the MPD from the normal environmental pressure (MSLP), where the MPD is defined as a function of the SST and Bister and Emanuel (2002) values are used to bound their values. On the other hand, we define the thermodynamic intensification factor following Holland (1997) with variables extracted along the synthetic tracks and use the SST–MPD relationship as a capping function (see Fig. 8). Table A3 summarizes the main differences of the two approaches.

As illustrated in Fig. 9, both methods produce a similar dependence of the intensification function (11) on the distance to the potential and maximum pressure drop. When P_c approaches the local maximum potential intensity (or maximum pressure drop), ΔP^c is more likely to be positive and to decrease the storm intensity. In other words, we can distinguish two phases (see Fig. 9, left graph): the intensification phase (in blue), when the central pressure is above the local MPI threshold, and a decay phase (in red), when the central pressure is below this local MPI.

In contrast to the MPD, the MPI does not represent the maximum achievable pressure and can be exceeded when accounting for additional external factors not reflected in climate data. For example, Chen et al. (2021) suggest that rapid intensification also depends on dynamical factors (e.g., upper divergence and wind shear). While the James and Mason (2005) formulation implicitly assumes that these factors are accounted for in the residual term of Eq. (12), it does not consider that the distance to the maximum potential thus defined can take negative values. Indeed, this specification originally suggests using the maximum achievable central pressure, i.e., that $P_{\mathrm{c}}-\mathrm{MPI}>\mathrm{0}$. As a result, the two intensification factors are not defined in the same domain (see Fig. 7). Using the local thermodynamic MPI, negative values, corresponding to a situation where the central pressure is below the MPI (P_c<MPI), are associated with a positive response of the pressure dynamic module and a decrease in storm intensity. Using the MPD, the response is likely to become positive when the distance to the MPD is below 40 to 60 hPa, depending on whether the function has been applied to the local sea surface temperature or the maximum per grid box ($\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }\times$ month). Given the similar response provided, using the local MPI instead of MPD as the intensification factor offers a better theoretical representation of the conditions affecting cyclone intensification in the cyclone dynamic module. The central pressure dynamics used for the fitting and dynamics of the synthetic tracks produced in the North Atlantic basin are illustrated in Fig. 10. The intensification of synthetic cyclones is in line with historical observations.

Figure 9 Response of the depression dynamics model to distance to thermodynamic potentials over intensification factor domains (hPa).

Figure 10 Individual historical vs. synthetic depression dynamics in the North Atlantic basin, and the confidence interval.

2.4.5 SDR over land

We model the evolution of the cyclone after landfall using an exponential decay function considering that tropical cyclone intensity decreases as a function of the time and distance the tropical cyclone has covered whilst being over land (Kaplan and DeMaria, 1995). Similarly to Bloemendaal et al. (2020), after three steps on land we suppose that the wind at time t_L follows:

$$\begin{array}{ll}\text{(13)}& {\displaystyle }{V}_{t_{\mathrm{L}}}& {\displaystyle }=V_b+(R\cdot V_{\mathrm{0}}-V_b)e^{-\mathit{\alpha }{t}_{\mathrm{L}}}-f_{\mathrm{1}}\left(t_{\mathrm{L}}\right)\left(\ln {\displaystyle \frac{D^l}{D_{\mathrm{0}}}}\right)+f_{\mathrm{2}}\left(t_{\mathrm{L}}\right){\displaystyle }& {\displaystyle }=V(t_{\mathrm{L}},D^l,V_{\mathrm{0}}),\end{array}$$

where D^l is the distance to coast computed using natural Earth coastlines (available at https://www.naturalearthdata.com/downloads/10m-physical-vectors/, last access: 14 October 2022), V₀ is the wind at landfall and t_L the time spent on land by the eye. This function was fitted on IBTrACS using nonlinear least squares. In our procedure, we use the global parameters: R=0.79, V_b=15 m s⁻¹, α=0.044 h⁻¹, and $f_{\mathrm{1}}\left(t_{\mathrm{L}}\right)={\stackrel{\mathrm{̃}}{c}}_{\mathrm{1}}{t}_{\mathrm{L}}(t_{\mathrm{0},\mathrm{L}}-t_{\mathrm{L}})$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{1}}=\mathrm{3.35}\times {\mathrm{10}}^{-\mathrm{4}}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{2}}$ , $t_{\mathrm{0},\mathrm{L}}=\mathrm{172}$ h, $f_{\mathrm{2}}\left(t_{\mathrm{L}}\right)=d_{\mathrm{1}}{t}_{\mathrm{L}}(t_{\mathrm{0},\mathrm{L}}-t_{\mathrm{L}})$, $d_{\mathrm{1}}=-\mathrm{0.00186}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{2}}$ and D₀=1 km. Kaplan and DeMaria (1995) introduced this function to model the decay of tropical cyclones over land in a simple way and showed that it provides an acceptable approximation for t_L⩾12 h. As each time step is 3 h, we let the TC intensity be driven by Eq. (11) and the first three steps and apply the decay function after three steps, that is, for t_L⩾12 h. A more sophisticated description could integrate for instance, cyclone physics, kinetic energy, and non-meteorological parameters such as ground topology. The SDR puts a strong constraint on the cyclone evolution after three steps. However, in the context of national damage assessment, we reiterate that reported damage costs are the combination of a series of various impacts including storm surge and not only extreme wind and that the most exposed area is at landfall. We consider therefore that the hypothesis of a rapid decay is acceptable and in line with observations.

2.5 Cyclone generation algorithm

The full cyclone generation procedure is presented in Algorithm 1. The cyclone wind speed is initiated at 20 m s⁻¹, and the initial pressure is determined from the WPR (Eq. 3). While the cyclone is over the sea, the pressure evolution ΔP_c is entirely determined from Eq. (11) based on the local thermodynamic potential. To prevent the model from producing an unrealistically low central pressure, we cap the MPD using Eq. (10). With this truncation, the lower bound for the pressure is given by the observed low-pressure values in similar sea surface temperature conditions. While the cyclone is over the sea, the wind is defined with the WPR (Eq. 3). When the cyclone arrives on land, the MPI is computed from the last known climate variables for the first three steps and the pressure still follows the relationship (11). After three steps (9 h) on land, we start applying the decay relationship (Eq. 13) to define the wind. The variations of longitude and latitude are always defined using Eqs. (1) and (2). We force cyclones to remain in their genesis basins in this exercise. For example, running the algorithm on IPSL climate projections between 2075 and 2100 in RCP85 produces the output plotted in Fig. 11.

Figure 11 Example of 100 representative years of synthetic tracks produced with CATHERINA on IPSL-CM5A-MR raw climate data in RCP85 between 2075 and 2100 (i.e., four runs over the 25-year period).

The cyclone intensification process used is inspired by the STORM model from Bloemendaal et al. (2020), which includes a single climate variable, and extended following Holland (1997) and Emanuel (1988) to encompass two more variables. We found that this extension provides statistically significant instrumental variables in the description of tropical cyclone intensification, which is the aim of the algorithm. Consequently, even if some thermodynamic processes have been simplified in our approach, it is still a step forward with respect to the existing state of the art in the field of integrated assessment modeling for climate impact analysis. Indeed, our approach is easy to implement, more sophisticated in terms of processes included than most existing IAMs, has low bias due to our state-of-the-art bias-correction module described below, and can integrate any CMIP simulation with a limited set of available variables (only a few vertical levels, some only available at a monthly timescale, with some variables not always available).

Algorithm 1 Cyclone generation algorithm.

V←20 m s⁻¹

$P_{\mathrm{c}}(s=\mathrm{0})\leftarrow \mathrm{MSLP}-{\left({\displaystyle \frac{V}{a}}\right)}^{\mathrm{1}/b}$

if $\mathrm{MSLP}-P_{\mathrm{c}}>\mathrm{0}$ & V>v_m & s_l<4 then

while the pressure is below normal, wind is above threshold and we are not on land, do

$$\begin{array}{ll}{\displaystyle }x\left(s\right)& {\displaystyle }\leftarrow x(s-\mathrm{1})+\mathrm{\Delta }x\left(s\right)\\ {\displaystyle }\text{ where }& {\displaystyle }\mathrm{\Delta }x\left(s\right)\propto \text{Eq}.\left(\mathrm{1}\right)\\ {\displaystyle }y\left(s\right)& {\displaystyle }\leftarrow y(s-\mathrm{1})+\mathrm{\Delta }y\left(s\right)\\ {\displaystyle }\text{ where }& {\displaystyle }\mathrm{\Delta }{y}_t\propto \text{Eq}.\left(\mathrm{2}\right)\\ {\displaystyle }\mathrm{MPI}\left(s\right)& {\displaystyle }\leftarrow f_{\mathrm{MPI}}\left(y\right(s),P_{\mathrm{c}}(s-\mathrm{1}),\mathrm{SST}(s),\\ {\displaystyle }& {\displaystyle }{T}_{\mathrm{tropo}}\left(s\right),\mathrm{MSLP}\left(s\right),\mathrm{RH}\left(s\right))\\ {\displaystyle }& {\displaystyle }{\mathrm{f}}_{\mathrm{MPI}}\propto \text{Eq. (4)}\\ {\displaystyle }{P}_{\mathrm{c}}\left(s\right)& {\displaystyle }\leftarrow \max \left(P_{\mathrm{c}}\right(s)+\mathrm{\Delta }{P}_{\mathrm{c}}(s),\mathrm{MSLP}(s)-\mathrm{MPD}(s\left)\right)\\ {\displaystyle }& {\displaystyle }\mathrm{MPD}\propto \text{Eq}.\left(\mathrm{10}\right)\mathit{\&}\\ {\displaystyle }& {\displaystyle }\mathrm{\Delta }{P}_{\mathrm{c}}\left(s\right)\propto \text{Eq}.\left(\mathrm{11}\right)\\ {\displaystyle }V\left(s\right)& {\displaystyle }\leftarrow a{\left(\mathrm{MSLP}-P_{\mathrm{c}}\left(s\right)\right)}^b\end{array}$$

if on land = TRUE then

$s_{\mathrm{l}}\leftarrow s_{\mathrm{l}}+\mathrm{1}$

end if

else

if s_l>4 then

Same functional for x and y but, compute distance to land D(s) from natural earth coastlines and do

$$\begin{array}{ll}{\displaystyle }V\left(s\right)& {\displaystyle }\leftarrow V_b+(R\cdot V_{\mathrm{0}}-V_b)e^{-\mathit{\alpha }{s}_{\mathrm{l}}}\\ {\displaystyle }& {\displaystyle }-f_{\mathrm{1}}\left(t_{\mathrm{L}}\right)\left(\ln {\displaystyle \frac{D}{D_{\mathrm{0}}}}\right)+f_{\mathrm{2}}\left(t_{\mathrm{L}}\right)\\ {\displaystyle }& {\displaystyle }V\left(s\right)\propto \text{Eq}.\left(\mathrm{13}\right)\end{array}$$

$s_{\mathrm{l}}\leftarrow s_{\mathrm{l}}+\mathrm{1}$

end if

end if

Note: this algorithm assumes step-wise extraction of climate data in the Monte Carlo process. Another way, closer to the framework suggested in Bloemendaal et al. (2020), would be to (i) compute the tracks without properties, (ii) retrieve all climate variables and (iii) determine the properties using the extracted climate conditions in the last step.

3 Exposure in the shared socioeconomic pathways

3.1 Physical asset exposure

Eberenz et al. (2020a) present a methodology to downscale physical asset values on a high-resolution grid using a combination of nightlight intensity, population data, and global country indicators and make their dataset fully available (Eberenz et al., 2020b). These estimates of physical asset values are based on the light intensity L_i – from nighttime lights of the Black Marble 2016 annual composite of the VIIRS day–night band (DNB) (Román et al., 2018) in 2016 at the 15 arcsec resolution – and the population P_pix per pixel – from the Gridded Population of the World (GPW) database (Center for International Earth Science Information Network (CIESIN), 2017) in 2015 at the 30 arcsec resolution for 224 countries and various additional sources⁷ allowing us to define the total asset (A_tot) for each country. This value is distributed to each grid cell proportionally to the light intensity L_i times the local population P_pix:⁸

$$\begin{array}{}\text{(14)}& A_{\mathrm{pix}}=A_{\mathrm{tot}}{\displaystyle \frac{L_{\mathrm{i}}\cdot P_{\mathrm{pix}}}{{\sum }_{{\mathrm{pix}}_i}^N\left(L_{\mathrm{i}}\cdot P_{{\mathrm{pix}}_i}\right)}}.\end{array}$$

The physical asset value is expressed in USD as of 2014. Using this dataset in the future requires correction (either simulated or reported damages) for inflation using, for instance, the consumer price index (available at https://fred.stlouisfed.org, last access: 14 October 2022).

This method for allocation of national assets has limitations. For example, the distribution of assets near the coast, industrial production sites or agricultural facilities may not be well represented. However, with this approach, asset values are defined on a uniform grid across countries and can be projected by multiplication by appropriate dynamic factors.

3.2 The SSP framework

Future exposure is sensitive to the scenarios of population growth and economic development. To take this into account, we use the framework of the SSPs introduced in O’Neill et al. (2014). These narratives are used in the IPCC development scenarios and provide a reference framework for risk assessment. A growing segment of the literature is dedicated to measuring the feasibility, costs and implications of achieving these scenarios (Riahi et al., 2017, 2021), and multiple IAMs were launched on assumptions based on these narratives (Riahi et al., 2017; Rogelj et al., 2018; Gidden et al., 2019).⁹ Figure 12 displays the projections of the main features used by CATHERINA, the global domestic production (GDP) and population in the five SSPs at the global level by the IIASA GDP model.¹⁰ CATHERINA uses these two indicators at the regional level (32 regions are available) to compute future local exposure (see Sect. 3.3).

Figure 12 SSP GDP and population variation until 2100, at the global scale, by the IIASA GDP model.

Let us recall the main assumptions underlying these narratives. The “middle road” pathway (SSP2) is used as the reference in most scenario analyses. It is a plausible baseline in terms of economic and social resiliency, in which the urbanization level is relatively high and GDP and population are constantly increasing following the observed historical trend. On the other hand, the “rocky road” or “national-rivalry” pathway (SSP3) presents totally different properties: relative stagnation of GDP with a strong increase in the population (see Fig. 12). SSP4 corresponds to a scenario with the highest inequality, and SSP5 is the most likely to lead to the higher concentration pathways (RCP8.5), with extensive use of fossil fuel reserves but higher economic development and global market integration.

Based on these storylines, it is clear that the physical exposure to tropical cyclones will be driven by different factors in different SSP scenarios. Scenarios with steady growth of GDP are generally associated with an increase in urbanization but a decrease in the global population by 2050. In general, all narratives, except SSP3, present relatively similar population dynamics at the global level. On the other hand, the scenario with the lowest economic growth (SSP3) assumes a sustained increase in the global population (cf. Fig. 12), in particular in rural areas. We can therefore expect, in the former case, the physical asset value exposure to be driven by the increase in regional wealth and mainly by the growth of the exposed population in the latter.

In this paper, we specify economic and climate parameters independently, while the literature generally associates specific SSPs and RCPs, in particular in the CMIP6 exercise. Indeed, integrated assessment modeling demonstrates that specific temperature targets can only be reached under certain socioeconomic conditions. The socioeconomic and representative concentration pathways are therefore intrinsically linked at the global scale. For example, CMIP6 refers to the following scenarios: SSP1–1.9, SSP1–2.6, SSP2–4.5, SSP3–7.0 and SSP5–8.5. On the other hand, this pairing is not straightforward at the regional level. In particular, the SSP2 scenario is associated with regional heterogeneity in socioeconomic pathways. As a result, although we can expect a convergence in the long run, we considered it relevant to use economic development scenarios independently of RCP in this exercise. When looking at the aggregated level, however, we consider only scenarios that are feasible (Rogelj et al., 2018). For instance, the RCP85 is achieved only in the conditions of the SSP5.

Figure 13 SSP matrix O’Neill et al. (2017) (based on O’Neill et al., 2014, Fig. 1, inspired by Kriegler et al., 2012, Fig. 3).

Lastly, we acknowledge another limitation of our approach: the vulnerability parameter (represented by the damage function parameter v_h in our framework) does not depend on the SSP, while we could expect a reduction in the vulnerability parameter in the scenarios where the adaptation challenges are limited, such as SSP5 (cf. Fig. 13). This question is left for further research.

3.3 Dynamic projection of local exposure in SSPs

To estimate future exposures along the cyclone tracks in each scenario, we use the downscaled estimation for the exposed wealth (Eberenz et al., 2020b) and the coefficients representing the change between the current state and the future scenario in the framework of the shared socioeconomic pathways (O’Neill et al., 2014, 2017; Jones and O'Neill, 2020). The local physical exposure at the coordinates (x,y) at time t in a region j in scenario k is defined as follows:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& \mathrm{\Phi }(x,y,j,k,t)=\\ & \underset{\text{Global macro factor}}{\underbrace{\left(F_{\text{GDP}}^{\mathrm{cap}}\right(j,k,t){)}^{{\mathit{\alpha }}_{\mathrm{1}}}}}\cdot \underset{\text{Local factor}}{\underbrace{\left(F_{\mathrm{pop}}\right(x,y,k,t){)}^{{\mathit{\alpha }}_{\mathrm{2}}}\cdot {\mathcal{L}}_P(x,y)}},\end{array}\end{array}$$

where ℒ_P(x,y) is the local population density from Eberenz et al. (2020b), and the factor $F_{\text{GDP}}^{\mathrm{cap}}$ is the projected GDP per capita growth for each region:

$$\begin{array}{}\text{(16)}& F_{\text{GDP}}^{\mathrm{cap}}(j,k,t)={\displaystyle \frac{\mathrm{GDP}(j,k,t)/\mathrm{GDP}(j,t=\mathrm{2020})}{P(j,k,t)/P(j,t=\mathrm{2020})}},\end{array}$$

where P is the total population of the region retrieved from the SSP database (Riahi et al., 2017) (https://tntcat.iiasa.ac.at/SspDb/, last access: 14 October 2022) and F_pop is the population exposure growth factor:

$$\begin{array}{}\text{(17)}& {\displaystyle }{F}_{\mathrm{pop}}(x,y,k,t)={\displaystyle \frac{p(x,y,k,t)}{p(x,y,t=\mathrm{2020})}},\end{array}$$

where $p(x,y,k,t)$ represents the local projections of population (Jones and O'Neill, 2020) illustrated in Fig. 15. Figure 14 displays the scenario-based projections of GDP and population in the five SSPs, at the regional level by the IIASA model.

We introduce the exponents α₁ and α₂ to disentangle the effects of increased cyclone intensity, GDP growth and population growth on the future damages as well as to account in a simple manner for possible future adaptation to tropical cyclone risk. Indeed, taking ${\mathit{\alpha }}_{\mathrm{1}}={\mathit{\alpha }}_{\mathrm{2}}=\mathrm{1}$ amounts to assuming that damages from cyclones of similar intensities grow proportionally to the local GDP, whereas newly accumulated wealth can be more resilient than the existing one and protected by additional adaptation measures. Indeed, Bakkensen and Mendelsohn (2019) test for evidence of adaptation in past cyclone damage data and find that α₁ is statistically different from 1 in all countries except the USA. At the global level, they estimate α₁ to be equal to 0.364, with a standard error of 0.175, which is significant at the 95 % confidence level.

Figure 14 Regional $F_{\mathrm{GDP}}^{\mathrm{cap}}$ factor variation in the SSP IIASA database (R32). Source: https://tntcat.iiasa.ac.at/SspDb/ (last access: 14 October 2022). The country mapping is provided in Table A5.

Figure 15 Local F_pop factor variation in all SSPs in 2100 (source: Jones and O'Neill, 2020, NASA Socio-Economic Data Application Center – SEDAC). The scenario-based population grid generation methodology is detailed by Jones and O'Neill (2020) with a last version downscaled at 1 km following Gao (2020). This population grid is available every 10 years. CATHERINA uses the closest value in the definition of the exposure.

4 Damage assessment at the national level

4.1 Damage modeling

The percentage of asset values destroyed by a tropical cyclone depends on multiple parameters. For example, to assess the vulnerability of specific infrastructures to tropical cyclones, precise descriptions of building vulnerability are provided in the Federal Emergency Management Agency (FEMA) reports. Unanwa et al. (2000) propose a series of wind-damage bands depending on building types. They show that the sensitivity is higher for commercial and institutional buildings than for residential and mid-rise buildings and that generalized damages occur above 43–60 m s⁻¹, while a sustained wind regime above 73 m s⁻¹ could lead to the destruction of the entire superstructure of most buildings (except for mid- and low-rise ones). These bottom-up approaches allow us to set the limits of the potential damage functions, but their use requires a complete inventory of assets and up-to-date values of numerous parameters (age, height, materials, etc.). Therefore, CATHERINA relies on regional damage functions calibrated by Eberenz et al. (2021) on wind speed along the cyclone track (IBTrACS) and reported damages in Guha-Sapir et al. (2018).

Damages provoked by tropical cyclones can be related to several sub-perils. While 40 % of cyclone damages are directly wind-related, another 40 % can generally be attributed to storm surge, and the rest of the damage is generated by heavy precipitation. However, CATHERINA does not distinguish these sub-perils but uses a statistical relationship to estimate the regional damage induced by a cyclone from a proxy variable given by the maximum wind speed. Indeed, the wind speed is the proxy used in the Saffir–Simpson hurricane wind scale to define the intensity of a cyclone. The damage function is fitted on multiple events from the total damage reported in the global disaster database EM-DAT (Guha-Sapir et al., 2018) (available at http://www.emdat.be (last access: 14 October 2022) This database, used in most studies on the topic, accounts for the total reported damage (sum of all sub-perils) and does not distinguish damages from sub-perils. Filtering the database by subtype “tropical cyclone” allows us to extract 1855 tropical cyclones in the period between 1980 and 2021, among which 1101 events have a reported total damage cost in USD (see Fig. 17). In terms of damage, tropical cyclones are, using the full set of observations from 1980 to 2021, the most damaging events reported (see Fig. 16). The database includes a start date field (day, month and year) allowing us to map 455 events, with the events reported in IBTrACS using the start year and month and country.¹¹ We use this database to validate our simplified estimation process.

Figure 16 Proportion of damage cost (total damage in USD) by disaster subtypes reported in EM-DAT. Using the number of people affected places tropical cyclones after riverine floods and droughts, and the number of deaths places ground movements in the first position (see Table A1 for details).

Figure 17 Number of tropical cyclones with reported damage cost in the EM-DAT database. The color scale indicates the average reported damage cost in USD in each country for each tropical cyclone.

4.2 An explicit damage function

Damage functions can take several different forms (Prahl et al., 2015), but the most common choice is a cubic functional of the wind speed. To estimate the fraction of loss from a storm with sustained wind speed V, Emanuel (2011) introduced the following formula:

$$\begin{array}{}\text{(18)}& f(V,v_h^j)={\displaystyle \frac{{\left(\max (V-v_{\mathrm{0}},\mathrm{0})\right)}^{\mathrm{3}}}{(v_h^j-v_{\mathrm{0}}{)}^{\mathrm{3}}+{\left(\max (V-v_{\mathrm{0}},\mathrm{0})\right)}^{\mathrm{3}}}},\end{array}$$

where f is the fraction of the property value lost, v₀=25.7 m s⁻¹, and $v_h^j$ is a parameter that needs to be calibrated for each region j. Figure 18 illustrates the shape of the damage function for different values of this parameter.

To account for adaptation, we could modify this function and introduce an additional threshold value. For instance, to account for local adaptation to the wind climate, Leckebusch et al. (2007) suggest scaling the wind value by the 98th percentile of the local wind speed distribution. However, the assumption that adaptation will always keep the damages from the 98th percentile wind at the same level is probably optimistic and would prevent us from using the model to estimate the required investment to balance future damages and economic shocks.

Figure 18 Fraction of property value lost as a function of wind speed (m s⁻¹), obtained using Eq. (18) with different values of v_h. Source: Emanuel (2011).

4.3 To a regional damage calibration

Using the reported damage estimates from EM-DAT crossed with cyclone tracks (IBTrACS) and geographical and socioeconomic information along these tracks, Lüthi (2019) refined the damage function approach using machine learning techniques introducing region-specific damage functions.

We recall the main steps of the methodology presented in Eberenz et al. (2021) to define the regional damage functions. The authors first defined the event damage ratio (EDR) as a ratio of simulated damage (SED) to normalized reported damage (NRD) for each cyclone:

$$\begin{array}{}\text{(19)}& \mathrm{EDR}(i,j)={\displaystyle \frac{\mathrm{SED}(i,v_h(j\left)\right)}{\mathrm{NRD}\left(i\right)}}.\end{array}$$

The total damage ratio (TDR) is then defined in each region by summing over events:

$$\begin{array}{}\text{(20)}& \mathrm{TDR}\left(j\right)={\displaystyle \frac{{\sum }_i\mathrm{SED}(i,v_h(j\left)\right)}{{\sum }_i\mathrm{NRD}\left(i\right)}}.\end{array}$$

For each event, there is a value for v_h allowing us to optimally calibrate the explicit damage function given in Eq. (18). The relatively wide distribution of v_h for the same country shows that there is a large uncertainty in the relationship between the wind speed and the corresponding fraction of losses. Figure A4 shows the uncertainty in regional damage functions depending on the optimization technique used, and Fig. A5 allows us to appreciate, for countries where more than five cyclones were reported, the spread of plausible damage functions.

Eberenz et al. (2021) propose two alternative optimization methodologies to find the value of $v_h^{\star }$ maximizing the prediction quality of the regional damages: root mean square fraction (RMSF), minimizing the spread of the EDRs,

$$\begin{array}{}\text{(21)}& {v_h}_{\mathrm{RMSF}}^{\star }\left(j\right)={\mathrm{argmin}}_j\exp \left(\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum {\left(\ln \left(\mathrm{EDR}\right(\mathrm{i}\left)\right)\right)}^{\mathrm{2}}}\right),\end{array}$$

and TDR, finding the value of $v_h^{\star }$, such that the ratio of total simulated damage – obtained summing over event damages – and total reported damage tends to 1.

$$\begin{array}{}\text{(22)}& {v_h}_{\mathrm{TDR}}^{\star }\left(j\right)={\mathrm{argmin}}_j\left|\mathrm{TDR}\left(j\right)-\mathrm{1}\right|\end{array}$$

The values of $v_h^{\star }$ obtained by Eberenz et al. (2021) with the two methods are given in Table 2. For most regions, the optimized curves are similar for the two optimization techniques, but the results diverge for the Philippines (WP2) and to a lesser extent for Mainland China (WP3) events. The case of the Philippines, for example, discussed in Eberenz et al. (2021), is explained by the large number of parameters involved in the damage estimation and emphasizes two main limitations of the model. First, this framework lacks an explicit representation of sub-perils that disrupt and damage several sectors and services. Second, differences in exposure and vulnerability between urban and rural areas exposed to tropical cyclones are likely to contribute to the large spread in EDR.

Table 2 Values of v_h (m s⁻¹) obtained using the TDR and RMSF methods for each region from the CLIMADA environment.

Coefficient from version 1.5 of the CLIMADA environment. Figure A4 also illustrates the shapes of the functions for the different optimization problems (RMSF vs. TDR) and version (1.0 vs. 1.5).

4.4 Simplified damage estimation along tracks

In the context of our national level assessment, we propose a simplified damage module. The simulated damage for a given cyclone – in both IBTrACS and our synthetic tracks – is computed using the following procedure for each individual cyclone. First, a uniform grid of physical asset values with steps given by the average cyclone radius is defined on the map of the affected area. The cyclone track is linearly interpolated, and the tiles affected by the cyclone (containing a part of the interpolated path) are identified (see Fig. 19).

Figure 19 Illustration of the high-resolution asset data (a) original physical asset value resolution. (b) Aggregation of the asset value over 0.25^∘ boxes to evaluate asset exposure along the interpolated cyclone track (here corresponding to the 2013 Hayan cyclone). (c) Illustration of damage calculation: damage is aggregated over the white boxes which correspond to cyclone locations over land.

Second, for each tile identified in the previous step, we retrieve the maximum wind speed V and compute the proportion of wealth lost $f(V,v_h^j)$ using the relation (18) with the total damage ratio parameter given in Eberenz et al. (2021). Then, we compute the total simulated damage by aggregating the physical asset exposure multiplied by the proportion of wealth lost on each tile over all tiles affected by the cyclone.

As a result of this procedure, we obtain the total simulated damage SED_i(j,t) caused by the ith cyclone in region j, simulated with climate variables for year t. Finally, the cyclone damage cost in region j and year t is simulated as follows:

$$\begin{array}{}\text{(23)}& \mathcal{D}(j,t)=\sum _i{\mathrm{SED}}_i(j,t),\end{array}$$

where the sum is taken over all cyclones occurring in a given year. This procedure can then be repeated many times to obtain the distribution of annual cyclone damages and to compute other statistics such as the mean and quantiles of this distribution.

The damage functions used in the second step are retrieved directly from the CLIMADA environment. These functions were fitted with the same physical asset value. However, in our case, we project these values onto a coarser grid (first step) in such a way that the extraction is simplified for a large number of synthetic tracks in the context of the present global exercise. To ensure that the estimated damages produced with this simplification are consistent with the historical records, we computed simulated damages over the historical tracks and compared the results to EM-DAT. We aggregate asset values on a 0.25 × 0.25^∘ grid. The spread between SED and NRD distribution remains important. To further reduce the errors, we thus divide the simulated damage by the average SED-to-NRD ratios in each region, which are computed using the RMSF damage function on IBTrACS and total reported damage from EM-DAT:

$$\begin{array}{}\text{(24)}& {\widehat{r}}_j={\displaystyle \frac{{\sum }_i\mathrm{SED}(i,v_h^{\mathrm{RMSF}}(j\left)\right)}{{\sum }_i\mathrm{NRD}}}.\end{array}$$

Figure 20 presents the estimated vs. actual damages computed using the RMSF damage function and the distribution of the re-scaled estimation by country using the intersection of IBTrACS with lands (762 events), crossed with EMDAT (606). Each dot represents the total damage over a year in a country (211 observations).

Figure 20 Estimated (RMSF) vs. reported damage (USD). The filled dots represent individual year–country pairs with reported damage in EM-DAT. This distribution is obtained using a 0.25 × 0.25 resolution projection. Crosses correspond to country averages (over all years).

5 Application to RCPs

5.1 Climate simulation bias correction for climate change application

The variables from climate model projections used by CATHERINA are subject to multiple biases. To reduce uncertainty caused by these biases, we use the cumulative distribution function transform (CDF-t) method developed in Vrac et al. (2012) and Michelangeli et al. (2009) to correct the distribution of each variable in each basin. Our bias-correction approach is the standard in the climate community (Navarro-Racines et al., 2020) (see http://ccafs-climate.org/bias_correction/, last access: 14 October 2022).

Consider a generic climate variable (denoted by χ) at a fixed location, which is available both from ERA5 reanalysis and from a given CMIP5 model. We are interested in two time periods: the historical period (covered by both the climate model and the reanalysis) and a future time period (covered only by the climate model). Let $F_{\mathrm{ERA}\mathrm{5}}^{\mathrm{h}}$ and $F_{\mathrm{CMIP}}^{\mathrm{h}}$ be the distribution functions of χ under reanalysis and under the climate model for the historical period, and $F_{\mathrm{CMIP}}^{\mathrm{f}}$ is the distribution function of χ under the climate model for the future period. The distribution function under the climate model is subject to much stronger biases than that under the reanalysis. The CDF-t method constructs the distribution function for χ with a reduced bias for the future time period, denoted by ${\widehat{F}}_{\mathrm{CMIP}}^{\mathrm{f}}$ and given by

$$\begin{array}{}\text{(25)}& {\widehat{F}}_{\mathrm{CMIP}}^{\mathrm{f}}(\cdot )=F_{\mathrm{ERA}\mathrm{5}}^{\mathrm{h}}\left(F_{\mathrm{CMIP}}^{\mathrm{h},-\mathrm{1}}\right(F_{\mathrm{CMIP}}^{\mathrm{f}}(\cdot )\left)\right),\end{array}$$

where $F_{\mathrm{CMIP}}^{\mathrm{h},-\mathrm{1}}$ is the inverse function of $F_{\mathrm{CMIP}}^{\mathrm{h}}$. For a given value ${\mathit{\chi }}_{\mathrm{CMIP}}^{\mathrm{f}}$ of the variable χ obtained for the future period from the climate model, the corresponding bias-corrected value ${\widehat{\mathit{\chi }}}_{\mathrm{CMIP}}^{\mathrm{f}}$ may then be computed via

$$\begin{array}{ll}{\displaystyle }{\widehat{\mathit{\chi }}}_{\mathrm{CMIP}}^{\mathrm{f}}& {\displaystyle }={\widehat{F}}_{\mathrm{CMIP}}^{\mathrm{f},-\mathrm{1}}\left(F_{\mathrm{CMIP}}^{\mathrm{f}}\right)\\ {\displaystyle }& {\displaystyle }=F_{\mathrm{CMIP}}^{\mathrm{f},-\mathrm{1}}\left(F_{\mathrm{CMIP}}^{\mathrm{h}}\right(F_{\mathrm{ERA}\mathrm{5}}^{\mathrm{h},-\mathrm{1}}\left(F_{\mathrm{CMIP}}^{\mathrm{f}}\right({\mathit{\chi }}_{\mathrm{CMIP}}^{\mathrm{f}}\left)\right)\left)\right).\end{array}$$

When the future period and the historical period coincide, the method reduces to the standard quantile transform:

$$\begin{array}{}\text{(26)}& {\displaystyle }{\widehat{\mathit{\chi }}}_{\mathrm{CMIP}}^{\mathrm{h}}=F_{\mathrm{ERA}\mathrm{5}}^{\mathrm{h},-\mathrm{1}}\left(F_{\mathrm{CMIP}}^{\mathrm{h}}\right({\mathit{\chi }}_{\mathrm{CMIP}}^{\mathrm{h}}\left)\right).\end{array}$$

First, we use the method on the historical period to compare the description of the thermodynamic potential and wind speed with and without correction, so Eq. (26) may be used directly. To extract the CDFs of the variables of interest, we generate synthetic track candidates from 1980 (beginning of ERA5) to 2010. We launch the simulation 10 times over these 30 years to obtain 300 representative years. By definition, for the genesis of the cyclones, the time of year and location are in line with historical cyclone data. However, in this module, the synthetic tracks are generated without climate constraints; i.e., cyclones are allowed to drift relatively far away from their genesis location (in the limits of their initial basin) and therefore can cover conditions which do not lead to the formation of tropical cyclones. At this stage, these tracks are not to be considered “tropical cyclone tracks” but “candidate” tracks. In the following stage, actual cyclone tracks will be generated from candidate tracks by filtering the ones where meteorological conditions for cyclone formation are satisfied. For each point in space and time along these synthetic tracks, we extract the values of the four climate variables from the reanalysis (ERA5) and from the historical simulations of the seven climate models. Then, by comparing the CDF of the climate variables estimated by the models with the reference CDF computed from the reanalysis, we determine the transformation allowing the values estimated by the models to better match those from ERA5.

The sea level pressure distributions are stable over basins and models. The tropospheric temperature and near-surface relative humidity distributions depend largely on the basin and display evidence of non-negligible model uncertainty. The sea surface temperature estimates along the same synthetic tracks in the historical period display much larger uncertainty. The northern Indian basin presents the widest uncertainty for all climate variables, which adds further uncertainty concerning the impact of climate change on tropical cyclones in this area.¹²

Individual variables entering the MPI computation are correlated as shown in Table 3. For example, sea surface temperature and tropopause temperature exhibit a negative correlation of 83 %. Therefore, applying bias correction to individual variables may lead to unrealistic combinations when evaluating the thermodynamic potentials. For example, extremely low tropopause temperatures associated with very high SSTs lead to unrealistically high lapse rates and therefore unrealistically large potentials. To overcome this issue we perform the bias correction on the MSLP, SST, thermodynamic efficiency factor ℰ, and relative humidity – which should not be strongly correlated to other variables according to the level of correlation present in the reanalysis – see Table 3. Figure 21 shows that this correction leads to similar distribution of the thermodynamic potentials in the models and the reanalysis.

Figure 21 Bias-correction module Q–Q plots (GCMs vs. ERA5) of maximum potential intensity (hPa). The Q–Q plots compare two distributions by plotting quantiles of one against the other.

Table 3 Correlation levels (global) of modeled variables affecting the MPI with their reference values in ERA5.

We apply the CDF-t correction technique along our historical synthetic tracks and compute the maximum potential intensity following Sect. 2.4.2. The pressure follows the dynamic process introduced in Sect. 2.4.4, and the corresponding wind is derived from the WPR (see Sect. 2.4.1). We define the model error as the relative error $\left(\frac{{\mathit{\chi }}_{\mathrm{CMIP}}-{\mathit{\chi }}_{\mathrm{ERA}\mathrm{5}}}{{\mathit{\chi }}_{\mathrm{ERA}\mathrm{5}}}\right)$ between the value produced by the model and the one produced by the ERA5 reanalysis. Figure 22 displays the average relative errors and shows that a 2 % relative error in the description of the maximum potential intensity can lead to more than 20 % error in the description of the implied wind compared to the result obtained with ERA5. This figure illustrates the efficiency of bias correction in the historical period. Indeed, on average, the CDF-t correction technique clearly reduces the error between the MPI estimated with climate reanalysis and the one computed from modeled climate data as well as (more importantly) the error in the description of the maximum wind speed. However, running CATHERINA on the different GCMs still produces a wide range of results in climate projections (cf. Fig. 23). Although the intensity of the storms increases with each model, underlying climate modeling uncertainty still strongly impacts the synthetic tracks produced.

Figure 22 Relative error at the level of maximum potential intensity (MPI) computed with ERA5 and climate data produced by the seven climate models for the historical period. The transparent bars represent original errors, and the color parts represent the residual relative error after CDF-t correction averaged over 30 years of cyclones.

Figure 23 Average and standard deviations of the maximum wind speed of tropical cyclones generated with the CMIP5 models after correction. Table A4 provides the average value of maximum wind speed (m s⁻¹).

5.2 Results in CMIP5 projections

The international climate modeling community introduced SSPs to translate varying narratives on the development of society in the long term. These projections impact the local physical asset value dynamics (Jones and O'Neill, 2020; Chen et al., 2020) and global macroeconomic variables (O’Neill et al., 2014, 2017). Under the assumptions of constant cyclone genesis frequency and a constant impact ratio (i.e., the damage functions remain the same), CATHERINA allows us to derive damage projections in varying climate and socioeconomic scenarios. Using bias-corrected climate variable projections from the seven climate models over the period 2070–2100, we provide an example of the application of the CATHERINA framework.¹³ As expected (cf. Mendelsohn et al., 2012, Fig. 3, for example), socioeconomic change leads to wider differences than climate change.

Figure 24 represents the expected value of future damage, and Table 4 provides the expected values (together with standard errors) along with the 50th, 66th and 95th percentiles of the global annual damage distribution in the simulations under various assumptions of future exposure growth (values of α₁ and α₂). In the discussion below, CATHERINA simulations of 300 representative years of synthetic tracks using ERA5 data over the 1980–2010 period are used as a baseline.

The choice ${\mathit{\alpha }}_{\mathrm{1}}={\mathit{\alpha }}_{\mathrm{2}}=\mathrm{1}$ assumes linear increase in exposure with respect to both GDP per capita and local exposed population. Using ${\mathit{\alpha }}_{\mathrm{1}}=\mathrm{1}/\mathrm{3}$, which is close to the value estimated in Bakkensen and Mendelsohn (2019), allows us to account for possible future adaptation. The configurations ${\mathit{\alpha }}_{\mathrm{1}}=\mathrm{1},{\mathit{\alpha }}_{\mathrm{2}}=\mathrm{0}$ and ${\mathit{\alpha }}_{\mathrm{1}}=\mathrm{0},{\mathit{\alpha }}_{\mathrm{2}}=\mathrm{1}$ allow us to decompose the risk contribution between GDP and exposed population. Finally, the choice ${\mathit{\alpha }}_{\mathrm{1}}=\mathrm{0},{\mathit{\alpha }}_{\mathrm{2}}=\mathrm{0}$ uses the current exposure value with no socioeconomic growth factor.

Assuming no future adaptation (${\mathit{\alpha }}_{\mathrm{1}}={\mathit{\alpha }}_{\mathrm{2}}=\mathrm{1}$), over the period 2070–2100, the RCP2.6 scenario, which is in line with the Paris Agreement and keeps global warming below 2 ^∘C by 2100, involves a growth of expected global annual financial losses from tropical cyclones by a factor of 4.2 on average. Ignoring socioeconomic and population growth factors (α₁=0 and α₂=0), our model suggests that the expected financial loss would grow by a factor of 1.6 due to increasing cyclone intensity. Taking into account adaptation, i.e., limiting the growth of damage with respect to GDP per capita (${\mathit{\alpha }}_{\mathrm{1}}=\mathrm{1}/\mathrm{3}$), the expected damage would grow by a factor of 2.6. In the case of SSP2–RCP4.5 (between 1.7 and 3.2 ^∘C warming by 2100) and SSP5–RCP8.5 (between 3.2 and 5.4 ^∘C warming by 2100), the average expected damage will be multiplied by 5.4 and 14.2, respectively, without adaptation. In RCP8.5 the expected damage will still grow by a factor of 2.8, ignoring the change in GDP per capita and population (α₁=0 and α₂=0). Interestingly, in SSP3, accounting for the population factor (that is, moving from α₂=0 to α₂=1) induces a decrease in global expected damage. This is due to a significant decrease in population in the USA, Australia, South Korea or Japan in the regions subject to tropical cyclones in this scenario. As these countries represent a large share of tropical-cyclone-related damages, the global expected damage is reduced.

Our expected damage estimates are subject to three types of uncertainty: internal climate variability, climate model uncertainty and socioeconomic uncertainty related to future exposure growth, adaptation measures and concentration pathways. The first two types of uncertainty are quantified by the standard errors in Table 4, while the last one may be evaluated by performing simulations under different assumptions about adaptation, SSP narratives and representative concentration pathways as illustrated in Fig. 24.

Figure 24 Expected value of global annual damage (in USD billion) in different SSP–RCP and exposure projection hypothesis configurations. The vertical dotted line corresponds to historical simulations with ERA5 data.

Table 4 Simulation statistics (in USD billion).

Note: ^* the couple SSP5–RCP26 exists in the integrated assessment modeling literature (Rogelj et al., 2018) but is not displayed because the SSP5 is more likely tied to high concentration scenarios.

Figure 25 Annualized regional expected damage (in USD billion) in SSP and RCP between 2070 and 2100 based on synthetic tracks produced with seven climate models (with bias correction) over 300 representative years launched independently. The scale differs between SSPs.

Figure 26 Annualized global damage (USD billion) in SSPs and RCPs between 2070 and 2100 based on synthetic tracks produced with seven climate models (with bias correction) over 250 representative years launched independently, assuming no adaptation (${\mathit{\alpha }}_{\mathrm{1}}={\mathit{\alpha }}_{\mathrm{2}}=\mathrm{1}$). Top range: 95 %–98 % (extremely unlikely losses), followed by 95 %–66 %, 50 %–66 % and 0 %–50 %. Synthetic historical damages are computed with 300 representative years of synthetic tracks generated with ERA5 between 1980 and 2010. The physical asset value is corrected by inflation between the cyclone year and 2014 (GDP reference year of Litpop). The global SSP–RCP configurations are presented as feasible in Rogelj et al. (2018).

Figure 25 displays the average annual damage per country, in the different shared socioeconomic pathways, assuming no future adaptation. We can see that the distributions across countries are slightly different from one SSP to another. Indeed, the distributions in SSP2 and SSP5 are similar, with higher expected damage in SSP5 because of the growth hypothesis this scenario relies on. However, SSP3 (rocky road) and SSP4 (inequality) are distributed differently. The scenario emphasizing inequalities – and its interpretation by scientists in terms of (i) socioeconomic developments (Riahi et al., 2017) and (ii) population distribution (Jones and O'Neill, 2016) – increases damage concentration in the United States. On the other hand, the rocky-road scenario, linked to a larger and more rural population, lower GDP and national rivalry, sees the damage more equally distributed among other nations.

Looking at the expected damage values is instructive, but because the aim of the model was also to stress test the resiliency of the financial and economic systems, it is important to examine the risk of extreme events, corresponding to higher quantiles of the loss distribution. Figure 26 and Table 4 show that the 95 % quantile of the loss distribution, corresponding to losses which arise, on average, once in 20 years, may be as high as 4–5 times the expected loss. This observation is in line with Coronese et al. (2019), who show that the impact of climate change is particularly striking for extreme events (see, for example, Fig. 2a in this reference).

5.3 Impact on sovereign bond spreads in RCPs

In this section, we briefly discuss the impact of cyclones on the sovereign bonds of the exposed economies. Particular attention will be given to emerging economies because of their higher vulnerability. We use an econometric model to relate the credit spread of sovereign bonds to the scenario-based distributions of damage developed in previous sections.

5.3.1 Sovereign credit spread model

The sovereign credit spread is a financial measure of the creditworthiness of a country, which is defined as the difference of yield between the country's sovereign bond and a risk-free asset, i.e., an AAA-rated bond. In practice, we use the spread with respect to the US 10-year treasury bond. Following Hilscher and Nosbusch (2010), we calibrate a cross-sectional econometric model for the credit spread based on annual end-of-year data:

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\mathrm{CS}}_t& =\mathit{\alpha }+{\mathit{\beta }}_{\mathrm{1}}\mathrm{\Delta }{\mathcal{C}}_t+{\mathit{\beta }}_{\mathrm{2}}{D}_t+{\mathit{\beta }}_{\mathrm{3}}{\mathrm{VIX}}_t+{\mathit{\beta }}_{\mathrm{4}}{r}_t^{\mathrm{10}\mathrm{Y}}\\ & +{\mathit{\beta }}_{\mathrm{5}}{\mathrm{TED}}_t+{\mathit{\beta }}_{\mathrm{6}}{\displaystyle \frac{L_t}{{\mathrm{GDP}}_t}}+{\mathit{\beta }}_{\mathrm{7}}{\displaystyle \frac{{\mathrm{reserves}}_t}{{\mathrm{GDP}}_t}}\\ & +{\mathit{\beta }}_{\mathrm{8}}{\mathcal{R}}_t+{\mathit{\beta }}_{\mathrm{9}}{\mathcal{K}}_t+{\mathit{\epsilon }}_t,\end{array}\end{array}$$

where CS is the end-of-year credit spread from the JP Morgan EMBI position report in BarraOne, 𝒞 is the commodity price index from Reuters Refinitiv, D is the average duration of the bonds from the JP Morgan EMBI position report in BarraOne, VIX is the CBOE volatility index from Reuters Refinitiv, r^10Y is the 10-year U.S. Treasury bond rate from Reuters Refinitiv, TED is the difference between the 3-month treasury bill and the 3-month LIBOR from Reuters Refinitiv, L refers to the total external debt stocks in USD from the World Bank, GDP refers to the end-of-year GDP from the World Bank, reserves are the total reserves including gold in USD, ℛ is the credit rating dummy variable, and 𝒦 is the country dummy variable. The model covers 74 countries between 2010 and 2020.

We chose to use all the bonds constituting the index, i.e., including multiple bonds per country per date. Table A7 shows the results of estimation of several submodels of our model. The model (1) assesses the effect of a variation of the debt $/$ GDP ratio only (this is the most important variable that we will use to quantify cyclone damages), the model (2) introduces one rating dummy (it equals zero for all bonds with ratings below B− and one otherwise), the model (3) introduces additional macro variables and the model (4) adds country effect dummies and bond duration.

In line with Edwards (1986), we find that the debt $/$ GDP ratio is significant in all the submodels. More importantly, we note that β₆ is relatively stable over the modeling frameworks, implying that the sensitivity of the spread to a sudden increase in the debt ratio does not strongly depend on external or idiosyncratic parameters. Including country effects in (4), we obtain a model reaching 71.5 % of the adjusted R².

5.3.2 Cyclone impact on sovereign spreads under representative concentration pathways

We now use the econometric model developed in the previous paragraph to assess the impact of tropical cyclones on emerging country bond spreads under relative concentration pathways. To this end, we make the following simplifying assumptions: (i) the bond spread model parameters remain stable over time, (ii) the cyclone damages are financed directly by the government by issuing new debt, and other variables of the model are not affected by cyclones, and (iii) only direct impact of cyclones is taken into account, and not the total economic costs. For each country j, we assess the annual bond spread variation due to cyclone damage in scenario k for the year t using the following formula:

$$\begin{array}{}\text{(28)}& {\mathrm{\Delta }}_{k,t}\mathrm{CS}(j,k,t)={\mathit{\beta }}_{\mathrm{6}}\times {\displaystyle \frac{\mathcal{D}(j,k,t)}{F_{\mathrm{GDP}}(j,k,t)\mathrm{GDP}(j,\mathrm{2020})}},\end{array}$$

where we recall that 𝒟 stands for annualized cyclone damage, and F_GDP is the GDP growth factor for the specified country/scenario. For example, an event simulated in 2077 with GISS-E2-H (NASA, USA) in RCP85 generated USD 158 billion damage in Australia. Normalized by the future projected domestic product in the scenario (e.g., reference value USD 1.3 trillion with growth coefficient 3.1, leading to 4.4 trillion), this damage represents about 7 % of the country's GDP in the year concerned. This would lead to a shock in the Australian sovereign bond yield of 27 basis points (one basis points equals 0.01 %), implying a considerable increase in financing cost.

Figure 27 Damage cost channeled to excess spread (cost of debt) of different countries, measured in basis points (one basis point =0.01 %).

We compare the spread variation defined from damages of the RCP26 baseline, RCP45 and RCP85 to obtain an annualized financial valuation of the cyclone-related physical climate risk. The impact of average annual cyclone damage on the spread of sovereign bonds of larger emerging countries, channeled by the impact on the debt-to-GDP ratio, is small. This is consistent with the historical data, as countries are not equally affected each year by tropical cyclones, and when they are, they develop a resiliency that is reflected in their damage functions. Therefore, we focus on the tail risk, that is, on the extreme event quantiles. Figure 27 represents the annual excess spread for a sample of countries, when the annual cyclone damages are given by the 90th, 95th and 98th percentiles of the distribution as well as by the maximum cyclone-related shock observed in the simulations. Recall that we have drawn 300 representative years from seven models with independent launches to increase the number of cyclone trajectories. On average, larger countries' excess spreads remain limited, but the maximum observed damages induce a significant excess cost of debt.

6 Conclusions

This paper proposes a structural framework to generate synthetic storms based on large-scale climate data and to produce simulations of annual cyclone damages per country. We show that, when used with reanalysis data and CMIP5 models over the historical period, our method produces tracks consistent with historical observations.

The synthetic tracks generated with our model have several applications. The first one is in natural disaster risk management to calibrate adaptation measures. For this purpose, the track-generation algorithm may be enhanced, for instance, by including dependency in the latitudinal and longitudinal incremental displacement, coupling with a meteorological forecasting model, or including ground topography to model the cyclone displacement over land. Another major field of application is to climate financial risk management, where this scenario-based event database can be used to evaluate physical risks and compute portfolio exposures. This would require us to better define asset level vulnerabilities.

The dataset used for the examples of this paper, based on low-resolution data and including a limited number of simulations, may not be accurate enough to properly calibrate adaptation measures. However, we believe that the framework presented here may be used to project a dense set of trajectories, compute expected damage and damage percentiles over the next decades and measure the investment required for adaptation and mitigation measures in the next 50 years. This work also reflects a practical exercise not carried out until now, of cross-referencing the latest datasets developed, putting into perspective both the socioeconomic and climatic development hypotheses and carrying out a bottom-up, rather than top-down, damage calculation.

Appendix A

Figure A1 Information in IBTrACS. The cyclones are split by nature (DS: disturbance; TS: tropical; TE: extratropical; ST: subtropical; NR: not reported; MX: mixture or contradicting nature reports from different agencies). We remove extratropical cyclones (75) and disturbance (205). The cyclones are reported in the following basins: eastern Pacific (EP); North America (NA); northern India (NI); South Atlantic (SA); southern India (SI); South Pacific (SP); western Pacific (WP). One can also explore sub-basin decomposition (Eberenz et al., 2021) (MM: missing – no sub-basin for this basin (no sub-basins provided for WP, SI); CS: Caribbean Sea; GM: Gulf of Mexico; CP: central Pacific; BB: Bay of Bengal; AS: Arabian Sea; WA: western Australia; EA: eastern Australia), but we chose to use basins. The different agencies worldwide report central pressure and maximum wind speed but sometimes use sometimes different standards. In particular, the reporting can vary in terms of sustained wind speed. According to the dataset documentation, the North Atlantic – US Miami (NOAA NHC) bureau (hurdat/atcf) gives the 1 min winds speed, while Tokyo, i.e., RSMC Tokyo (JMA), provides directly the 10 min sustained wind speed. (Similarly, newdelhi corresponding to RSMC New Delhi (IMD) gives the 3 min wind speed, reunion RSMC La Reunion (MFLR), the Australian TCWCs (TCWC Perth, Darwin, Brisbane) (BOM), and the RSMC Nadi (FMS); TCWC Wellington (NZMS) provides the 10 min sustained wind speed and (CMA) 2 min sustained wind.) The lack of reporting standards between agencies is a source of uncertainty in the input data.

Figure A2 Map of longitude variation mean coefficients fitted on a 5×5 grid grouped per month.

Figure A3 Map of latitude variation mean coefficients fitted on a 5×5 grid grouped per month.

Figure A4 Regional optimal damage functions of the CLIMADA package (Eberenz et al., 2021).

Figure A5 Country calibration of damage functions (when the countries have been hit more than five times).

Figure A6 Evidence of model uncertainty. A comparison of the properties obtained by generating 30 years of tropical cyclones with different models from the CMIP5 raw climate data vs. ERA5. The average number of cyclones making landfall is computed by averaging the number of events with a number of steps on land positive (s_l>0). Maximum pressure drop and wind are, respectively, computed by extracting the maximum value for the corresponding variables of the cyclone tracks. The light blue bars represent the mean of the variable of interest among the IBTrACS-filtered dataset (with the confidence interval drawn from the standard deviation of the distribution). The yellow bars represent the same variable extracted from the synthetic data generated by Algorithm 1 using ERA5 data. In terms of average values, the models produce consistent tracks in every basin. Then, we compare the output of this algorithm with different climate data produced on the historical period by the climate models taking part in the fifth phase of the CMIP. A general observation is the poorer performance of the model in the northern Indian basin. This could be due to the smaller number of intense cyclones remaining in the sample after the filtration by intensity (1.6 storms per year with winds exceeding 35 m s⁻¹).

Figure A7 Cumulative distribution functions per basin for the variables of interest along synthetic tracks produced with ERA5 and extracted (at the same location) from climate data produced by the seven climate models. NA: North America, EP: eastern Pacific, NI: northern India, SI: southern India, SP: South Pacific, WP: western Pacific. The bias-correction module is indeed fitted on a larger range of climate conditions. By definition, for the genesis of the cyclones, the time of year and the location are in line with historical cyclone data. However, in the bias-correction module, the synthetic tracks are generated without climate constraints; i.e., cyclones are allowed to drift relatively far away from their genesis location (in the limits of their initial basin) and therefore can cover conditions which do not lead to the formation of tropical cyclones. At this stage, these tracks are not to be considered “TC tracks” but “candidate” tracks. In the following stage, TC tracks will be generated from candidate tracks by filtering those ones where meteorological conditions for cyclone formation are satisfied.

Table A1 EM-DAT reporting proportions^a. Tropical cyclone-related statistics were emphasized.

Notes: ^a proportion excludes damages whose subtypes are not reported. ^b Reconstruction costs are not well reported over disaster subtypes. ^c Brush, bush, pasture. Among disaster subtypes, tropical cyclones present a higher quality in reporting and represent a large share of total damages.

Table A2 Wind pressure relationship coefficients (i.e., V=a (ΔP)^b, a in ${\mathrm{hPa}}^{-\mathrm{1}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$).

Holland (1997)Bister and Emanuel (2002)

Table A3 Comparison STORM–CATHERINA cyclone intensification module.

Table A4 Mean of the maximum wind speed of synthetic cyclones in representative representation pathways (m s⁻¹).

Table A5 Mapping countries to IIASA regions.

Table A6 Examples of countries' damage statistics in USD billion (α₁=1, α₂=1).

Table A7 Simple credit spread model.

Note: ${}^{*}p<\mathrm{0}$.1; ${}^{**}p<\mathrm{0}$.05; ${}^{***}p<$0.01. Dropping Argentina in the constant in (4).