Compare Gauge Maxima over all Sources¶
This notebook illustrates one approach to comparing the 36 Cascadia CoPes Hub Ground Motions and Tsunami Sources at one of the Lagoon Creek gauge locations.
The maximum water depth is plotted for each of the 36 events, both in the usual alphabetical order of event names and also after ordering from smallest to largest values.
The sources have been assigned weights in the logic tree that is shown and discussed in the webpage linked above. These 36 weights sum to 1, and might be viewed as conditional probabilities, given that a single major CSZ earthquake occurs. These weights can be used to determine, for example, the median inundation depth over all the events, or the depth corresponding to any quantile.
Note that this is not a full Probabilistic Tsunami Hazard Assessment (PTHA), since it is based on the assumption that exactly one CSZ earthquake has occurred, and does not take into account the annual probability that this happens, or the possibility that more than one earthquake could occur in a given time period. The notebook HazardCurveDemo discusses how PTHA might be performed, although we note that these particular sources were not originally designed for PTHA purposes.
However, the analysis in this notebook might be useful to those studying the paleoseismic evidence, for example, where the question is more focused on the chance that a past CSZ tsunami could have reached a particular point, rather than hazard mitigation.
%matplotlib inline
from pylab import *
from clawpack.pyclaw.gauges import GaugeSolution
Specify the set of events and the logic tree weight 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
print('Event weights from logic tree:')
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}')
Event weights from logic tree:
BL10D: 0.03750
BL10M: 0.06250
BL10S: 0.02500
BL13D: 0.03750
BL13M: 0.06250
BL13S: 0.02500
BL16D: 0.03750
BL16M: 0.06250
BL16S: 0.02500
BR10D: 0.03750
BR10M: 0.06250
BR10S: 0.02500
BR13D: 0.03750
BR13M: 0.06250
BR13S: 0.02500
BR16D: 0.03750
BR16M: 0.06250
BR16S: 0.02500
FL10D: 0.01250
FL10M: 0.02083
FL10S: 0.00833
FL13D: 0.01250
FL13M: 0.02083
FL13S: 0.00833
FL16D: 0.01250
FL16M: 0.02083
FL16S: 0.00833
FR10D: 0.01250
FR10M: 0.02083
FR10S: 0.00833
FR13D: 0.01250
FR13M: 0.02083
FR13S: 0.00833
FR16D: 0.01250
FR16M: 0.02083
FR16S: 0.00833
The weights sum to 1.000
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.
Set gaugeno to 1, 2, or 3 for one of the onshore gauges. Note that at offshore gauges a clearer comparison might be to plot the maximum sea surface elevation rather than water depth.
#gauges_dir = 'geoclaw_gauges_to_upload'
gauges_dir = 'geoclaw_gauges_all_events'
gaugeno = 2
hmax = []
events = array(all_events)
print(f'Maximum water depth at Gauge {gaugeno} for each event:')
for event in events:
fname_gauge = f'{gauges_dir}/{event}_gauge{gaugeno:05d}.txt'
#print(fname_gauge)
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)
Maximum water depth at Gauge 2 for each event: BL10D: 10.162 BL10M: 13.501 BL10S: 3.581 BL13D: 3.854 BL13M: 8.365 BL13S: 8.745 BL16D: 6.998 BL16M: 5.155 BL16S: 3.731 BR10D: 12.482 BR10M: 6.072 BR10S: 12.225 BR13D: 0.289 BR13M: 5.406 BR13S: 5.431 BR16D: 0.000 BR16M: 0.000 BR16S: 0.657 FL10D: 28.884 FL10M: 14.991 FL10S: 3.823 FL13D: 0.553 FL13M: 3.296 FL13S: 13.467 FL16D: 9.980 FL16M: 2.462 FL16S: 4.882 FR10D: 10.877 FR10M: 6.129 FR10S: 10.613 FR13D: 9.018 FR13M: 1.584 FR13S: 3.311 FR16D: 4.699 FR16M: 2.119 FR16S: 3.119
Plot these values for each event¶
figure(figsize=(8,10))
ievents = range(len(events))
plot(hmax, ievents)
yticks(ievents, events)
gca().invert_yaxis() # so first event is on top
grid(True)
xlabel('meters')
title(f'Maximum depth h');
Sort from smallest to largest hmax¶
This gives a nicer plot and is needed for computing quantiles.
ind = argsort(hmax)
hmax_sort = hmax[ind]
events_sort = events[ind]
weights_sort = weights[ind]
figure(figsize=(8,9))
ievents = range(len(events)-1,-1,-1)
plot(hmax_sort, ievents)
yticks(ievents, events_sort);
grid(True)
xlabel('meters')
title(f'Maximum depth h (after sorting events by h)');
Plot as a bar chart with an indication of weights¶
The width of each bar in the bar chart will be proportional to the weight of that event from the logic tree.
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
Compute cumulative weights and quantiles:¶
wcum = 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):
wcum[k] = wcum[k+1] + weights_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))
quantiles = [50,80]
wq = {}
hq = {}
events_wq = {}
for q in quantiles:
wq[q] = 1 - 0.01*q
k = where(wcum >= wq[q])[0].max()
print(f'{len(wcum) - k} events contribute to the {q}% quantile depth:')
events_wq[q] = events_sort[k-1:]
print(events_wq[q])
hq[q] = hmax_sort[k]
print(f'with {q}% quantile inundation depth = {hq[q]:.2f} meters\n')
19 events contribute to the 50% quantile depth: ['BL16M' 'BR13M' 'BR13S' 'BR10M' 'FR10M' 'BL16D' 'BL13M' 'BL13S' 'FR13D' 'FL16D' 'BL10D' 'FR10S' 'FR10D' 'BR10S' 'BR10D' 'FL13S' 'BL10M' 'FL10M' 'FL10D'] with 50% quantile inundation depth = 5.41 meters 10 events contribute to the 80% quantile depth: ['FL16D' 'BL10D' 'FR10S' 'FR10D' 'BR10S' 'BR10D' 'FL13S' 'BL10M' 'FL10M' 'FL10D'] with 80% quantile inundation depth = 10.16 meters
fig,ax = subplots(figsize=(8,5))
step(hmax_step, wcum, 'b', where='pre')
xlim(-0.1, hmax_step.max())
ylim(0, 1)
grid(True)
ax.set_xlabel('Exceedence values for max water depth')
ax.set_ylabel('Probability of single event exceeding')
# plot dashed lines at specific quantiles:
qcolors = ['g','r','m']
for k,q in enumerate(quantiles):
label = f'{q}% quantile depth = {hq[q]: .2f}m'
ax.plot(hmax_step, wq[q]*ones(hmax_step.shape), '--', color=qcolors[k], label=label)
# plot vertical dashed lines showing depths:
ax.plot([hq[q],hq[q]], [0,wq[q]], '--', color=qcolors[k])
ax.legend(loc='upper right', framealpha=1, fontsize=10)
title(f'Cumulative distribution of depths for {len(events)} events');
Summary¶
The plot above shows the distribution of maximum inundation depth at one particular gauge location. The probability on the vertical axis goes from 0 to 1 since it is based on the assumption that a single CSZ event happens, chosen from the 36 CoPes Hub sources with probability corresponding to its weight in the logic tree, and these weights sum to 1. The plot is a step function with steps corresponding to each event, and the vertical jump at each step is the weight of that event. It is in fact simply the upper envelope of the bar chart shown earlier.
This curve is similar to a hazard curve that might come out of a PTHA analysis, but that curve would show the annual proability of exceeding each water depth. It is more complicated to compute, since one must also allow for the possibility that no event occurs or more than one event occurs in a given year. One must also make assumptions on the annual probability of any major CSZ event occuring and about the probability distribution of a major event.
For more discussion and an illustration of how a hazard curve could be computed, see
HazardCurveDemo, rendered from the notebook geoclaw_multirun/HazardCurveDemo.ipynb in the repository.