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())
Reading fault geometry from ./cascadia30.mshout

Header:

#   fault No. , centroid(lon,lat,z[km]) , node1(lon,lat,z[km]) , node2(lon,lat,z[km]) , node3(lon,lat,z[km]) , mean vertex length(km) , area(km^2) , strike(deg) , dip(deg)

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])
Loaded geometry for 963 triangular subfaults

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])
[-125.36073    46.302805 8795.382   ]
[ -125.311684    46.526803 10094.313   ]
[ -125.038743    46.426986 13526.373   ]

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]:
<matplotlib.text.Text at 0x1a1ddd1438>

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)
Reading earthquake data from _cascadia.001297.rupt

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);
Will create dtopo on arrays of shape 106 by 181
Making Okada dz for each of 963 subfaults
0..1..2..3..4..5..6..7..8..9..10..11..12..13..14..15..16..17..18..19..20..21..22..23..24..25..26..27..28..29..30..31..32..33..34..35..36..37..38..39..40..41..42..43..44..45..46..47..48..49..50..51..52..53..54..55..56..57..58..59..60..61..62..63..64..65..66..67..68..69..70..71..72..73..74..75..76..77..78..79..80..81..82..83..84..85..86..87..88..89..90..91..92..93..94..95..96..97..98..99..100..101..102..103..104..105..106..107..108..109..110..111..112..113..114..115..116..117..118..119..120..121..122..123..124..125..126..127..128..129..130..131..132..133..134..135..136..137..138..139..140..141..142..143..144..145..146..147..148..149..150..151..152..153..154..155..156..157..158..159..160..161..162..163..164..165..166..167..168..169..170..171..172..173..174..175..176..177..178..179..180..181..182..183..184..185..186..187..188..189..190..191..192..193..194..195..196..197..198..199..200..201..202..203..204..205..206..207..208..209..210..211..212..213..214..215..216..217..218..219..220..221..222..223..224..225..226..227..228..229..230..231..232..233..234..235..236..237..238..239..240..241..242..243..244..245..246..247..248..249..250..251..252..253..254..255..256..257..258..259..260..261..262..263..264..265..266..267..268..269..270..271..272..273..274..275..276..277..278..279..280..281..282..283..284..285..286..287..288..289..290..291..292..293..294..295..296..297..298..299..300..301..302..303..304..305..306..307..308..309..310..311..312..313..314..315..316..317..318..319..320..321..322..323..324..325..326..327..328..329..330..331..332..333..334..335..336..337..338..339..340..341..342..343..344..345..346..347..348..349..350..351..352..353..354..355..356..357..358..359..360..361..362..363..364..365..366..367..368..369..370..371..372..373..374..375..376..377..378..379..380..381..382..383..384..385..386..387..388..389..390..391..392..393..394..395..396..397..398..399..400..401..402..403..404..405..406..407..408..409..410..411..412..413..414..415..416..417..418..419..420..421..422..423..424..425..426..427..428..429..430..431..432..433..434..435..436..437..438..439..440..441..442..443..444..445..446..447..448..449..450..451..452..453..454..455..456..457..458..459..460..461..462..463..464..465..466..467..468..469..470..471..472..473..474..475..476..477..478..479..480..481..482..483..484..485..486..487..488..489..490..491..492..493..494..495..496..497..498..499..500..501..502..503..504..505..506..507..508..509..510..511..512..513..514..515..516..517..518..519..520..521..522..523..524..525..526..527..528..529..530..531..532..533..534..535..536..537..538..539..540..541..542..543..544..545..546..547..548..549..550..551..552..553..554..555..556..557..558..559..560..561..562..563..564..565..566..567..568..569..570..571..572..573..574..575..576..577..578..579..580..581..582..583..584..585..586..587..588..589..590..591..592..593..594..595..596..597..598..599..600..601..602..603..604..605..606..607..608..609..610..611..612..613..614..615..616..617..618..619..620..621..622..623..624..625..626..627..628..629..630..631..632..633..634..635..636..637..638..639..640..641..642..643..644..645..646..647..648..649..650..651..652..653..654..655..656..657..658..659..660..661..662..663..664..665..666..667..668..669..670..671..672..673..674..675..676..677..678..679..680..681..682..683..684..685..686..687..688..689..690..691..692..693..694..695..696..697..698..699..700..701..702..703..704..705..706..707..708..709..710..711..712..713..714..715..716..717..718..719..720..721..722..723..724..725..726..727..728..729..730..731..732..733..734..735..736..737..738..739..740..741..742..743..744..745..746..747..748..749..750..751..752..753..754..755..756..757..758..759..760..761..762..763..764..765..766..767..768..769..770..771..772..773..774..775..776..777..778..779..780..781..782..783..784..785..786..787..788..789..790..791..792..793..794..795..796..797..798..799..800..801..802..803..804..805..806..807..808..809..810..811..812..813..814..815..816..817..818..819..820..821..822..823..824..825..826..827..828..829..830..831..832..833..834..835..836..837..838..839..840..841..842..843..844..845..846..847..848..849..850..851..852..853..854..855..856..857..858..859..860..861..862..863..864..865..866..867..868..869..870..871..872..873..874..875..876..877..878..879..880..881..882..883..884..885..886..887..888..889..890..891..892..893..894..895..896..897..898..899..900..901..902..903..904..905..906..907..908..909..910..911..912..913..914..915..916..917..918..919..920..921..922..923..924..925..926..927..928..929..930..931..932..933..934..935..936..937..938..939..940..941..942..943..944..945..946..947..948..949..950..951..952..953..954..955..956..957..958..959..960..961..962..
Done
In [11]:
x.shape,y.shape, dtopo0.dZ.shape
Out[11]:
((106,), (181,), (100, 181, 106))
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]:
<matplotlib.text.Text at 0x1a4d7c5390>

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()))
Created cascadia001297.dtt3, with dynamic rupture of a Mw 9.09 event

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)
Figure files for each frame will be stored in  _plots
In [18]:
anim = J.make_anim(plotdir)
anim
Out[18]:


Once Loop Reflect