Generating dtopo files for CSZ fakequakes¶
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 at https://zenodo.org/record/59943#.WgHuahNSxE4.
This sample uses one rupture scenario extracted from data/cascadia.001297.
In [1]:
%matplotlib inline
In [2]:
from clawpack.geoclaw import dtopotools
import numpy as np
from copy import copy
import matplotlib.pyplot as pl
import os
Set up CSZ geometry¶
In [3]:
fault_geometry_file = './cascadia30.mshout'
print('Reading fault geometry from %s' % fault_geometry_file)
print('\nHeader:\n')
print(open(fault_geometry_file).readline())
In [4]:
# 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
In [5]:
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.
In [6]:
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:
In [7]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy
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 = numpy.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 = numpy.array(c)
ax.plot(c[:,0],c[:,1], 'b')
ax.set_aspect(1./numpy.cos(45*numpy.pi/180.))
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_title('Plan view')
Out[7]:
Rupture scenario¶
We now read in rupture scenario, using data from [https://zenodo.org/record/59943#.WgHuahNSxE4].
In [8]:
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.
In [9]:
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.
In [10]:
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);
In [11]:
x.shape,y.shape, dtopo0.dZ.shape
Out[11]:
In [12]:
fig,(ax0,ax1,ax2, ax3) = pl.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.
In [13]:
pl.figure(figsize=(14,8))
pl.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):
pl.plot([lat,lat],[rupture_times[j],rupture_times[j]+rise_times[j]],'b')
pl.xlabel('latitude')
pl.ylabel('seconds')
pl.title('rupture time + rise time of each triangle vs. latitude')
Out[13]:
Plot the deformation as a function of time at a few locations¶
Same longitude, increasing latitude...
In [14]:
pl.figure(figsize=(10,9))
for kk in range(5):
pl.subplot(5,1,kk+1)
j = 60
k = 50 + 25*kk
pl.title('dynamic deformation at fixed location (' + '{:4.2f}'.format(x[j]) \
+ ',' + '{:4.2f}'.format(y[k]) + ')' )
pl.plot(dtopo0.times,dtopo0.dZ[:,k,j]);
pl.ylabel('dZ (meters)');
pl.xlabel('time (seconds)');
pl.tight_layout()
Create a dtopo file for GeoClaw¶
In [15]:
ruptno = rupt_fname.split('.')[1]
fname = 'cascadia' + ruptno + '.dtt3'
dtopo0.write(fname, dtopo_type=3)
In [16]:
print('Created %s, with dynamic rupture of a Mw %.2f event' % (fname, fault0.Mw()))
Animate the rupture:¶
In [17]:
from clawpack.visclaw.JSAnimation import IPython_display
import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J
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
plotdir = '_plots'
J.make_plotdir(plotdir, clobber=True)
fig = pl.figure(figsize=(12,5))
for k,t in enumerate(dtopo0.times[::10]):
plot_subfaults_dZ(t,fig)
J.save_frame(k, verbose=False)
pl.close(fig)
In [18]:
anim = J.make_anim(plotdir)
anim
Out[18]: