Making hazard curves for PTHA¶
The 36 See Cascadia CoPes Hub Ground Motions and Tsunami Sources were not originally designed for Probabilistic Tsunami Hazard Assessment (PTHA), and may not cover the full range of CSZ earthquakes. Moreover, this is a relatively small set of sources.
However, the relative likelihood of these particular events have been assigned as weights in the logic tree that is shown and discussed in the webpage linked above. The notebook CompareGaugeMaxima shows how these weights might be used in computing e.g. the median event at a particular location (from this particular set).
The 36 weights sum to 1, and might be viewed as conditional probabilities (i.e. given that a major CSZ earthquake occurs, the weight is the relative probability of each). These weights can be combined with a presumed annual probability that such an event occurs in order to obtain an annual probability of each event. These can then be used (together with tsunami simulations of each event) to perform PTHA (with the caveats mentioned above).
One typical product of a PTHA analysis is a "hazard curve" that plots the annual probability that some quantity of interest exceeds a threshold as a function of that threshold. This notebook illustrates how to compute the hazard curve for the maximum water depth at one of the onshore gauges for the Lagoon Creek test problem.
We also show how to make plots that illustrate the how the inundation depth varies over the 36 events after sorting them, which helps to see how much variation there is between these events, and to pick out the events that give the greatest depths. We can also disaggregate, to clearly see which events lead to the the 2475-year hazard, for example.
%matplotlib inline
from pylab import *
from clawpack.pyclaw.gauges import GaugeSolution
Specify the set of events and the annual probability of each¶
We use the weights from the logic tree shown in Cascadia CoPes Hub Ground Motions and Tsunami Sources.
depths = ['D','M','S']
# buried_locking events:
all_events = [f'BL10{depth}' for depth in depths] \
+ [f'BL13{depth}' for depth in depths] \
+ [f'BL16{depth}' for depth in depths] \
# add random events:
all_events += [e.replace('L','R') for e in all_events]
# add Frontal Thrust events:
all_events += [e.replace('B','F') for e in all_events]
all_events.sort()
print(f'all_events contains {len(all_events)} events')
all_events contains 36 events
event_weights = {}
for event in all_events:
w = 1/3 * 0.5
if 'B' in event:
w *= 0.75
else:
w *= 0.25
if 'D' in event:
w *= 0.3
elif 'M' in event:
w *= 0.5
else:
w*=0.2
event_weights[event] = w
#print(f'{event}: {w:.5f}')
weights = array([event_weights[e] for e in all_events])
print(f'The weights sum to {weights.sum():.3f}')
The weights sum to 1.000
Probabilities for a Poisson process¶
We make the common assumption that earthquakes are a Poisson process with some rate $\lambda$, meaning that the probability of an earthquake in any small time interval of length $\Delta t$ is $1 - \exp(-\lambda\Delta t)$. If $t$ is measured in years then the probability of an earthquake in 1 year ($\Delta t =1$) is $$ 1 - e^{-\lambda} = \lambda - \frac 1 2 \lambda^2 + \cal{O}(\lambda^3) \approx \lambda $$ for small $\lambda$, so $\lambda$ is roughly the annual probability of an earthquake.
Moreover it can be shown that the expected time to the next earthquake is $1/\lambda$, so this quantity is the "return time" or "recurrence interval".
Note that the Poisson process is independent of the time since the last event, so the annual probability is the same in 2026 as it was in 1701, just after the last big earthquake.
We use CSZ_return_time = 526 years for illustration below, the same value as used in the National Seismic Maps, but one could argue that a smaller value should be used at this point in time, based on the time since the last big one.
The annual probability of an event with weight $w$ from the logic tree is then $1 - \exp(-w/526)$. For small $w$ this is approximately equal to $w/526$, but not exactly.
Table of event weights and probabilities¶
The table below shows the weight for each of the 36 events and the annual probability of that event, along with the return time for that event (which are very long since the probability of each single event is very small).
CSZ_return_time = 526 # years
# annual probability of each event occuring:
#event_probs = weights / CSZ_return_time # this is WRONG in quadratic terms
event_probs = 1 - exp(-weights/CSZ_return_time) # correct for Poisson process
print(f'\nAssuming CSZ_return_time = {CSZ_return_time} years:\n')
print(f'event weight annual probability return time (years)')
for k,event in enumerate(all_events):
print(f'{event}: {weights[k]:.5f} {event_probs[k]:.10f} {CSZ_return_time/weights[k]: 8.0f}')
Assuming CSZ_return_time = 526 years: event weight annual probability return time (years) BL10D: 0.03750 0.0000712902 14027 BL10M: 0.06250 0.0001188142 8416 BL10S: 0.02500 0.0000475274 21040 BL13D: 0.03750 0.0000712902 14027 BL13M: 0.06250 0.0001188142 8416 BL13S: 0.02500 0.0000475274 21040 BL16D: 0.03750 0.0000712902 14027 BL16M: 0.06250 0.0001188142 8416 BL16S: 0.02500 0.0000475274 21040 BR10D: 0.03750 0.0000712902 14027 BR10M: 0.06250 0.0001188142 8416 BR10S: 0.02500 0.0000475274 21040 BR13D: 0.03750 0.0000712902 14027 BR13M: 0.06250 0.0001188142 8416 BR13S: 0.02500 0.0000475274 21040 BR16D: 0.03750 0.0000712902 14027 BR16M: 0.06250 0.0001188142 8416 BR16S: 0.02500 0.0000475274 21040 FL10D: 0.01250 0.0000237640 42080 FL10M: 0.02083 0.0000396063 25248 FL10S: 0.00833 0.0000158427 63120 FL13D: 0.01250 0.0000237640 42080 FL13M: 0.02083 0.0000396063 25248 FL13S: 0.00833 0.0000158427 63120 FL16D: 0.01250 0.0000237640 42080 FL16M: 0.02083 0.0000396063 25248 FL16S: 0.00833 0.0000158427 63120 FR10D: 0.01250 0.0000237640 42080 FR10M: 0.02083 0.0000396063 25248 FR10S: 0.00833 0.0000158427 63120 FR13D: 0.01250 0.0000237640 42080 FR13M: 0.02083 0.0000396063 25248 FR13S: 0.00833 0.0000158427 63120 FR16D: 0.01250 0.0000237640 42080 FR16M: 0.02083 0.0000396063 25248 FR16S: 0.00833 0.0000158427 63120
Combining probabilities¶
Note that if one earthquake event has annual probability $p_1$ and a different event has probability $p_2$ then the probability that either one of those events (or both) happen in a year is 1 minus the probability that neither event happens, i.e. $$ p = 1 - (1-p_1)(1-p_2) = p_1 + p_2 - p_1p_2. $$ This is approximately equal to $p_1 + p_2$ if both probabilities are very small, but not at all the same for large probabilities.
Also note that if the events are Poisson with rates $\lambda_1$ and $\lambda_2$ then $p_j = 1 - \exp(-\lambda_j)$ and so the formula above gives $$ p = 1 - \exp(-\lambda_1)\exp(-\lambda_2) = 1 - \exp(-(\lambda_1+\lambda_2)), $$ so the rates simply add together, i.e. $p = 1 - \exp(-\lambda)$ with $\lambda = \lambda_1 + \lambda_2$.
These formulas and variants of it are used repeatedly below in computing probabilities for combinations of events, e.g. all events exceeding some threshold.
Load the gauge data for all events¶
This assumes the script make_gauge_upload.py has already been run to make .txt files for each event/gauge combination in the format used by the Code Verification Platform.
gauges_dir = 'geoclaw_gauges_to_upload'
gaugeno = 2
hmax = []
events = array(all_events)
for event in events:
fname_gauge = f'{gauges_dir}/{event}_gauge{gaugeno:05d}.txt'
gauge_data = loadtxt(fname_gauge, comments='#', skiprows=15)
h = gauge_data[:,1]
hmax_k = h.max()
#print(f'{event}: {hmax_k:.3f}')
hmax.append(hmax_k)
hmax = array(hmax)
figure(figsize=(8,10))
ievents = range(len(events))
plot(hmax, ievents)
yticks(ievents, events);
grid(True)
xlabel('meters')
title(f'Maximum depth h');
Sort from smallest to largest hmax¶
This gives a nicer plot and is needed for creating hazard curves.
ind = argsort(hmax)
hmax_sort = hmax[ind]
events_sort = events[ind]
probs_sort = event_probs[ind]
weights_sort = weights[ind]
Plot as a bar chart within indication of conditional probabilities¶
heights_sort = 14 * weights_sort # scale to get reasonable bar widths
fig,ax = subplots(figsize=(8,9))
# plot short bars to left of 0 to show width also for events with hmax=0:
left_offset = -0.5
ax.barh(events_sort, left_offset, left=left_offset,
height=heights_sort, color='b', alpha=0.2)
ax.set_xlim(2*left_offset, 1.05*hmax_sort.max())
ax.barh(events_sort, hmax_sort,
height=heights_sort, color='b')
ax.invert_yaxis()
ax.grid(axis='x')
ax.set_title(f'Lagoon Creek Gauge {gaugeno}\n'
'hmax for each event, bar width proportional to weight');
fname = f'LagoonCreek_Gauge{gaugeno:05d}_depth_bars.png'
savefig(fname)
print('Created ',fname)
Created LagoonCreek_Gauge00002_depth_bars.png
Make hazard curve¶
Compute the cumulative weights pcum[k] is the probability of at least one event happening that meets or exceeds the depth hmax_sorted[k].
pcum = zeros(len(hmax_sort)+1) # with one extra 0 to start sum at end, and for plotting
for k in range(len(events_sort)-1,-1,-1):
pcum[k] = pcum[k+1] + probs_sort[k] - pcum[k+1] * probs_sort[k]
# add one more data point with h > hmax.max() and probability 0, for tail of step function:
hmax_step = hstack((hmax_sort, hmax_sort.max()+1))
Functions mapping probability p to return time rt or vice versa:¶
The 1e-9 floor is to avoid divide by 0 when using to plot a second axis.
#p2rt = lambda p: -1/log(1-p)
#rt2p = lambda rt: 1 - exp(-1/rt)
p2rt = lambda p: -1/log(1-maximum(p,1e-9))
rt2p = lambda rt: 1 - exp(-1/maximum(rt,1e-9))
Figure out events contributing to 975- or 2475-year hazard:¶
Events with return times of 975 or 2475 are sometimes called 1000-year or 2500-year events. The modified values come from the fact that for a Poisson process, a return time of $1/\lambda = 975$ results in a 5% chance that it will happen in the next 50 years, whereas $1/\lambda = 2475$ gives a 2% chance in the next 50 years, and this is how the hazard is more properly defined:
1 - exp(-50/975), 1 - exp(-50/2475)
(0.04998931898973191, 0.01999932662657944)
rt = 2475
p = rt2p(rt)
k = where(pcum >= p)[0].max()
print(f'{len(pcum) - k} events contribute to {rt}-year hazard:')
events_rt2475 = events_sort[k-1:]
print(events_rt2475)
h2475 = hmax_sort[k]
print(f'2475-year inundation depth = {h2475:.2f} meters\n')
rt = 975
p = rt2p(rt)
k = where(pcum >= p)[0].max()
print(f'{len(pcum) - k} events contribute to {rt}-year hazard:')
events_rt975 = events_sort[k-1:]
print(events_rt975)
h975 = hmax_sort[k]
print(f'975-year inundation depth = {h975:.2f} meters')
10 events contribute to 2475-year hazard: ['FL16D' 'BL10D' 'FR10S' 'FR10D' 'BR10S' 'BR10D' 'FL13S' 'BL10M' 'FL10M' 'FL10D'] 2475-year inundation depth = 10.16 meters 19 events contribute to 975-year hazard: ['BL16M' 'BR13M' 'BR13S' 'BR10M' 'FR10M' 'BL16D' 'BL13M' 'BL13S' 'FR13D' 'FL16D' 'BL10D' 'FR10S' 'FR10D' 'BR10S' 'BR10D' 'FL13S' 'BL10M' 'FL10M' 'FL10D'] 975-year inundation depth = 5.41 meters
Plot hazard curve as step function¶
fig,ax = subplots(figsize=(8,5))
step(hmax_step, pcum, 'b', where='pre')
xlim(-0.1, hmax_step.max())
ylim(2e-6, .002)
grid(True)
ax.set_xlabel('Exceedence values for max water depth')
ax.set_ylabel('Annual probability of at least one event exceeding')
# Add second axis label for return time:
secax = gca().secondary_yaxis('right', functions=(p2rt,rt2p))
secax.set_yticks([10000, 5000, 3000, 2000, 1000, 600])
secax.set_ylabel('return time (years)');
# plot dashed lines at specific return times:
label = f'975-year depth = {h975: .2f}m'
ax.plot(hmax_step, rt2p(975)*ones(hmax_step.shape), 'g--', label=label)
label = f'2475-year depth = {h2475: .2f}m'
ax.plot(hmax_step, rt2p(2475)*ones(hmax_step.shape), 'r--', label=label)
# plot vertical dashed lines showing depths:
ax.plot([h2475,h2475], [0,rt2p(2475)], 'r--')
ax.plot([h975,h975], [0,rt2p(975)], 'g--')
ax.legend(loc='upper right', framealpha=1, fontsize=10)
title(f'Hazard curve for {len(events)} events');
Check that hazard curve approaches expected probability at y-axis:¶
The way we defined the cumulative sum pcum of events above, the value pcum[0] should always be the annual probability of any of the 36 events happening, which is the value
we took for CSZ_return_time.
print(f'total probability = pcum[0] = {pcum[0]:.9f}')
print(f'should agree with 1 - exp(-1/526) = {1 - exp(-1/526):.9f}')
total probability = pcum[0] = 0.001899335 should agree with 1 - exp(-1/526) = 0.001899335
Plot color-coded events for these return periods:¶
heights_sort = 7500 * probs_sort # scale to get reasonable height plot
fig,ax = subplots(figsize=(8,9))
colors = len(events)*['b']
for k,event in enumerate(events_sort):
if event in events_rt2475:
colors[k] = 'r'
elif event in events_rt975:
colors[k] = 'g'
ax.barh(events_sort,hmax_sort,height=heights_sort,color=colors)
left_offset = -0.5
ax.barh(events_sort, left_offset, left=left_offset,
height=heights_sort, color=colors, alpha=0.2)
ax.set_xlim(2*left_offset,30);
ax.invert_yaxis()
ax.set_title('hmax for each event, bar width proportional to probability\n'
'red events contribute to 2475-year hazard\n'
'red + green events to 975-year hazard');
ax.grid(axis='x')
Make side by side plots¶
fig,axs = subplots(1,2,figsize=(12,7))
ax = axs[0]
ax.barh(events_sort,hmax_sort,height=heights_sort,color=colors)
left_offset = -0.5
ax.barh(events_sort, left_offset, left=left_offset,
height=heights_sort, color=colors, alpha=0.2)
ax.set_xlim(2*left_offset, 1.05*hmax_sort.max());
ax.invert_yaxis()
ax.grid(axis='x')
ax.set_xlabel('Max water depth of each event at specified location')
ax.set_title('hmax for each event, bar width proportional to probability\n'
'red events contribute to 2475-year hazard\n'
'red + green events to 975-year hazard');
ax = axs[1]
ax.step(hmax_step, pcum, 'b', where='pre')
ax.set_ylim(2e-6, 0.002)
ax.set_xlabel('Exceedence values for max water depth')
ax.set_ylabel('Annual probability of at least one event exceeding')
ax.grid(True)
ax.set_title(f'Lagoon Creek Gauge {gaugeno}\nHazard curve for {len(events)} events');
# Add second axis label for return time:
secax = gca().secondary_yaxis('right', functions=(p2rt,rt2p))
secax.set_yticks([10000, 5000, 3000, 2000, 1000, 600])
secax.set_ylabel('return time (years)');
# plot dashed lines at specific return times:
label = f'975-year depth = {h975: .2f}m'
ax.plot(hmax_step, rt2p(975)*ones(hmax_step.shape), 'g--', label=label)
label = f'2475-year depth = {h2475: .2f}m'
ax.plot(hmax_step, rt2p(2475)*ones(hmax_step.shape), 'r--', label=label)
ax.legend(loc='upper right', framealpha=1, fontsize=10);
# plot vertical dashed lines showing depths:
ax.plot([h2475,h2475], [0,rt2p(2475)], 'r--')
ax.plot([h975,h975], [0,rt2p(975)], 'g--')
ax.legend(loc='upper right', framealpha=1, fontsize=10)
tight_layout()
fname = f'LagoonCreek_Gauge{gaugeno:05d}_depth_HazardCurve.png'
savefig(fname)
print('Created ',fname)
Created LagoonCreek_Gauge00002_depth_HazardCurve.png
