Generating dtopo files for CSZ fakequakes¶
This notebooks demonstrates how the GeoClaw dtopotools module can be used to generate the dtopo file for a kinematic rupture specified on a set of triangular subfaults (available starting in Clawpack Version 5.5.0).
This uses one of the 1300 "fakequake" realizations from the paper
- Kinematic rupture scenarios and synthetic displacement data: An example application to the Cascadia subduction zone by Diego Melgar, R. J. LeVeque, Douglas S. Dreger, Richard M. Allen, J. Geophys. Res. -- Solid Earth 121 (2016), p. 6658. doi:10.1002/2016JB013314.
This requires cascadia30.mshout containing the geometry of the triangulated fault surface, from
https://github.com/dmelgarm/MudPy/blob/master/examples/fakequakes/3D/cascadia30.mshout
It also requires a rupture scenario in the form of a .rupt file from the collection of fakequakes archived at https://zenodo.org/record/59943#.WgHuahNSxE4.
This sample uses one rupture scenario extracted from data/cascadia.001297.
Version¶
Animation revised 2020-07-26 to run with v5.7.0
%matplotlib inline
from clawpack.geoclaw import dtopotools
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from copy import copy
from clawpack.visclaw import animation_tools
import os
from IPython.display import HTML
Set up CSZ geometry¶
fault_geometry_file = './cascadia30.mshout'
print('Reading fault geometry from %s' % fault_geometry_file)
print('\nHeader:\n')
print(open(fault_geometry_file).readline())
# read in .mshout (CSZ geoemetry)
cascadia = np.loadtxt(fault_geometry_file,skiprows=1)
cascadia[:,[3,6,9,12]] = 1e3*abs(cascadia[:,[3,6,9,12]])
print('Loaded geometry for %i triangular subfaults' % cascadia.shape[0])
For example, the first triangular fault in the given geometry of CSZ has the nodes
print(cascadia[0,4:7])
print(cascadia[0,7:10])
print(cascadia[0,10:13])
Plot the subfaults:¶
Set up a fault model with these subfaults, without yet specifying a particular earthquake scenario.
fault0 = dtopotools.Fault()
fault0.subfaults = []
nsubfaults = cascadia.shape[0]
for j in range(nsubfaults):
subfault0 = dtopotools.SubFault()
node1 = cascadia[j,4:7].tolist()
node2 = cascadia[j,7:10].tolist()
node3 = cascadia[j,10:13].tolist()
node_list = [node1,node2,node3]
subfault0.set_corners(node_list,projection_zone='10T')
fault0.subfaults.append(subfault0)
Now we can plot the triangular subplots:
fig = plt.figure(figsize=(15,10))
#ax = fig.add_subplot(121, projection='3d')
ax = fig.add_axes([.05,.05,.65,.9], projection='3d')
for s in fault0.subfaults:
c = s.corners
c.append(c[0])
c = np.array(c)
ax.plot(c[:,0],c[:,1],-c[:,2]/1000.,color='b')
ax.view_init(10,60)
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_zlabel('Depth (km)')
ax.set_title('Triangular subfaults')
#ax = fig.add_subplot(122)
ax = fig.add_axes([.75,.05,.2,.9])
for s in fault0.subfaults:
c = s.corners
c.append(c[0])
c = np.array(c)
ax.plot(c[:,0],c[:,1], 'b')
ax.set_aspect(1./np.cos(45*np.pi/180.))
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_title('Plan view')
Rupture scenario¶
We now read in rupture scenario, using data from [https://zenodo.org/record/59943#.WgHuahNSxE4].
rupt_fname = '_cascadia.001297.rupt'
print("Reading earthquake data from %s" % rupt_fname)
rupture_parameters = np.loadtxt(rupt_fname,skiprows=1)
This data is used to set the slip and rake on each of the subfaults loaded above. Since this is a dynamic rupture, we also set the rupture_time and rise_time of each subfault.
fault0 = dtopotools.Fault()
fault0.subfaults = []
fault0.rupture_type = 'kinematic'
rake = 90. # assume same rake for all subfaults
J = int(np.floor(cascadia.shape[0]))
for j in range(J):
subfault0 = dtopotools.SubFault()
node1 = cascadia[j,4:7].tolist()
node2 = cascadia[j,7:10].tolist()
node3 = cascadia[j,10:13].tolist()
node_list = [node1,node2,node3]
ss_slip = rupture_parameters[j,8]
ds_slip = rupture_parameters[j,9]
rake = np.rad2deg(np.arctan2(ds_slip, ss_slip))
subfault0.set_corners(node_list,projection_zone='10T')
subfault0.rupture_time = rupture_parameters[j,12]
subfault0.rise_time = rupture_parameters[j,7]
subfault0.rake = rake
slip = np.sqrt(ds_slip ** 2 + ss_slip ** 2)
subfault0.slip = slip
fault0.subfaults.append(subfault0)
Compute seafloor deformations with GeoClaw¶
We now run the create_dtopography routine to generate dynamic seafloor deformations at a given set of times. This applies the Okada model to each of the subfaults and evaluates the surface displacement on the grid given by x,y, at each time. These are summed up over all subfaults to compute the total deformation.
x,y = fault0.create_dtopo_xy(dx = 4/60.)
print('Will create dtopo on arrays of shape %i by %i' % (len(x),len(y)))
tfinal = max([subfault1.rupture_time + subfault1.rise_time for subfault1 in fault0.subfaults])
times0 = np.linspace(0.,tfinal,100)
dtopo0 = fault0.create_dtopography(x,y,times=times0,verbose=True);
x.shape,y.shape, dtopo0.dZ.shape
fig,(ax0,ax1,ax2, ax3) = plt.subplots(ncols=4,nrows=1,figsize=(16,6))
fault0.plot_subfaults(axes=ax0,slip_color=True,plot_box=False);
ax0.set_title('Slip on Fault');
X = dtopo0.X; Y = dtopo0.Y; dZ_at_t = dtopo0.dZ_at_t
dz_max = dtopo0.dZ.max()
t0 = 0.25*tfinal # time to plot deformation
dtopotools.plot_dZ_colors(X,Y,dZ_at_t(t0),axes=ax1,
cmax_dZ = dz_max, add_colorbar=False);
ax1.set_title('Seafloor at time t=' + str(t0));
t0 = 0.5*tfinal # time to plot deformation
dtopotools.plot_dZ_colors(X,Y,dZ_at_t(t0),axes=ax2,
cmax_dZ = dz_max, add_colorbar=False);
ax2.set_title('Seafloor at time t=' + str(t0));
t0 = tfinal # time to plot deformation
dtopotools.plot_dZ_colors(X,Y,dZ_at_t(t0),axes=ax3,
cmax_dZ = dz_max, add_colorbar=True);
ax3.set_title('Seafloor at time t=' + str(t0));
#fig.savefig('CSZ_triangular.png');
Plot the rupture time and rise time of each subfault¶
This shows where the rupture originates and how it propagates outward. Each vertical bar shows the rupture time and duration of one subfault.
plt.figure(figsize=(14,8))
plt.axes()
latitudes = [s.latitude for s in fault0.subfaults]
rise_times = [s.rise_time for s in fault0.subfaults]
rupture_times = [s.rupture_time for s in fault0.subfaults]
for j,lat in enumerate(latitudes):
plt.plot([lat,lat],[rupture_times[j],rupture_times[j]+rise_times[j]],'b')
plt.xlabel('latitude')
plt.ylabel('seconds')
plt.title('rupture time + rise time of each triangle vs. latitude')
Plot the deformation as a function of time at a few locations¶
Same longitude, increasing latitude...
plt.figure(figsize=(10,9))
for kk in range(5):
plt.subplot(5,1,kk+1)
j = 60
k = 50 + 25*kk
plt.title('dynamic deformation at fixed location (' + '{:4.2f}'.format(x[j]) \
+ ',' + '{:4.2f}'.format(y[k]) + ')' )
plt.plot(dtopo0.times,dtopo0.dZ[:,k,j]);
plt.ylabel('dZ (meters)');
plt.xlabel('time (seconds)');
plt.tight_layout()
Create a dtopo file for GeoClaw¶
ruptno = rupt_fname.split('.')[1]
fname = 'cascadia' + ruptno + '.dtt3'
dtopo0.write(fname, dtopo_type=3)
print('Created %s, with dynamic rupture of a Mw %.2f event' % (fname, fault0.Mw()))
Animate the rupture:¶
This fails in v5.5.0 due to a typo that was later fixed in clawpack/geoclaw commit 9422715, and should be fixed in v5.6.0.
dz_max = abs(dtopo0.dZ).max()
# Incorporate this function in dtopotools to replace animate_dz_colors?
def plot_subfaults_dZ(t, fig):
fig.clf()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
fault0.plot_subfaults(axes=ax1, slip_color=True, plot_box=False,
slip_time=t)
dtopo0.plot_dZ_colors(axes=ax2, t=t, cmax_dZ=dz_max)
return fig
figs = []
times = dtopo0.times[::5] # only use every 5th time for animation
if dtopo0.times[-1] not in times:
times = np.hstack((times, dtopo0.times[-1])) # include final dz
print('Animation will include %i times' % len(times))
for k,t in enumerate(times):
fig = plt.figure(figsize=(12,5))
plot_subfaults_dZ(t,fig)
figs.append(fig)
plt.close(fig)
anim = animation_tools.animate_figs(figs, figsize=(12,6))
HTML(anim.to_jshtml())