Hazard: Tropical cyclone surge from linear wind-surge relationship and a bathtub model

The TCSurgeBathtub class models surges generated by tropical cyclones. Given an elevation data set and a TropCyclone instance, it computes the surges for each historical and/or synthetic event at every centroid. TCSurgeBathtub inherits from Hazard and has an associated hazard type TCSurgeBathtub.

Model description

As a first approximation, the tropical cyclone’s wind field in each grid cell is used as an input to a simplified version of the wind-surge relationship in Xu (2010), which is based on pre-run SLOSH outputs.

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

# conversion factors
mph2ms = 0.44704;
f2m = 0.3048;

# the points read from the SLOSH graph
v0 = 60*mph2ms;
v1 = 140*mph2ms;
s0 = 6*f2m;
s1 = 18*f2m;

# the parameters for the linear function: a*(v-v0)+s0
a = (s1-s0)/(v1-v0)

# graphical representation
v = np.arange(20, 100)
vmph = np.arange(60, 141, 20)

plt.plot(v, a*(v-v0)+s0, label='approximation')
plt.scatter(vmph*mph2ms, a*(vmph*mph2ms-v0)+s0, label='SLOSH points')
plt.xlabel('wind speed (m/s)')
plt.ylabel('surge height (m)')
plt.legend()
plt.show()

The elevation of the centroids is then subtracted from the surge using the user-specified elevation data set. The elevation data set has to be given as a path to a GeoTIFF grid data file that covers the region affected by the tropical cyclone. A global data set is freely available as SRTM15+V2.0 (a sample of which is used in the example below). The improved-quality CoastalDEM data set is available on request from Climate Central.

In a final step, a decay of the surge height depending on the distance from the coastline by 0.2 meters per kilometer is implemented following Pielke and Pielke (1997).

Optionally, a user-specified sea level rise offset is added to the result.

Example

We compute the surges of Sidr 2007 and Roanu 2016 over Bangladesh as follows:

%matplotlib inline
# 1: tracks retrieval
from climada.hazard import TCTracks

tr_usa = TCTracks.from_ibtracs_netcdf(provider='usa', storm_id=['2007314N10093', '2016138N10081']) # SIDR 2007 and ROANU 2016
tr_usa.equal_timestep(0.5)
ax = tr_usa.plot()
ax.get_legend()._loc = 2 # correct legend location
ax.set_title('SIDR and ROANU'); # set title

2021-07-09 15:24:31,990 - climada.hazard.tc_tracks - INFO - Progress: 50%
2021-07-09 15:24:32,004 - climada.hazard.tc_tracks - INFO - Progress: 100%
2021-07-09 15:24:32,015 - climada.hazard.tc_tracks - INFO - Interpolating 2 tracks to 0.5h time steps.

# 2: wind gusts computation
from climada.hazard import TropCyclone, Centroids

# define centroids raster
min_lat, max_lat, min_lon, max_lon = 20, 27, 88.5, 92.5
cent_bang = Centroids.from_pnt_bounds((min_lon, min_lat, max_lon, max_lat), res=0.015)

tc_bang = TropCyclone.from_tracks(tr_usa, centroids=cent_bang)

2021-07-09 15:24:33,987 - climada.util.coordinates - INFO - Sampling from /home/tovogt/.climada/data/GMT_intermediate_coast_distance_01d.tif
2021-07-09 15:24:34,077 - climada.hazard.trop_cyclone - INFO - Mapping 2 tracks to 125424 coastal centroids.
2021-07-09 15:24:41,455 - climada.hazard.trop_cyclone - INFO - Progress: 100%

# 3: surge computation
from climada_petals.hazard import TCSurgeBathtub
from climada.util.constants import DEMO_DIR

# If you have the global SRTM15+V2.0 elevation data set, you can replace the following
# sample DEM data set by your SRTM15+V2.0.tiff:
topo_path = DEMO_DIR.joinpath('SRTM15+V2.0_sample.tiff')
ts_bang = TCSurgeBathtub.from_tc_winds(tc_bang, topo_path)

2021-07-09 15:24:41,466 - climada.util.coordinates - INFO - Sampling from /home/tovogt/.climada/demo/data/SRTM15+V2.0_sample.tiff

# plot elevation of the raster
elevation = ts_bang.centroids.get_elevation(topo_path)
ts_bang.centroids.plot(c=elevation, vmin=0, vmax=10)

2021-07-09 15:24:41,544 - climada.util.coordinates - INFO - Sampling from /home/tovogt/.climada/demo/data/SRTM15+V2.0_sample.tiff

<GeoAxesSubplot:>

# plot wind and surge SIDR
ax = tc_bang.plot_intensity(1)
ax.set_title('Wind SIDR 2007')
ax = ts_bang.plot_intensity(1)
ax.set_title('Surge SIDR 2007');

# plot wind and surge ROANU
ax = tc_bang.plot_intensity(2)
ax.set_title('Wind ROANU 2016')
ax = ts_bang.plot_intensity(2)
ax.set_title('Surge ROANU 2016');