-
Notifications
You must be signed in to change notification settings - Fork 63
Pisco earthquake 2: Velocity structure and fault geometry files
Now that the folder structure has been created we need to setup the relevant files for the inversion. Here we will setup the velocity structure file, which contains information on the 1D layered Earth model assumed in the inversion and the fault geometry files.
Mudpy uses velocity model files denoted by the extension .mod
a file template.mod
is copied into the structure
folder when you initalize the folder structure so that you may see it's format. You can either edit this file or make your own. The columns are layer thickness in km, vs in km/s, vp in km/s, density in g/cm^3, s-wave quality and, p-wave quality. Make sure you do not leave an empty line at the end of the file, this will produce empty Green's functions.
Now we must decide on a velocity structure for the Pisco region. A quick search of the literature reveals that for northern Chile and southern Peru a good 1D model can be found in Husen et al. (1999) which gives P and S wave velocities, that velocity model does not include density which is necessary for slip inversion. To compute density we use the empirical relation ρ=0.77+0.32*Vp. We also need P and S wave attenuation (Qa and Qb), we use values of 1500 and 600 as recommended in Pritchard et al. (2007). We arrive at the following velocity model:
Thickness (km) | Vs (km/s) | Vp (km/s) | Density (g/cm^3) | Qb | Qa |
---|---|---|---|---|---|
3.0 | 2.990 | 5.210 | 2.437 | 600 | 1500 |
2.0 | 3.090 | 5.370 | 2.488 | 600 | 1500 |
2.0 | 3.190 | 5.550 | 2.546 | 600 | 1500 |
2.0 | 3.290 | 5.720 | 2.600 | 600 | 1500 |
2.0 | 3.390 | 5.890 | 2.655 | 600 | 1500 |
4.5 | 3.440 | 5.980 | 2.684 | 600 | 1500 |
5.0 | 3.750 | 6.800 | 2.946 | 600 | 1500 |
5.0 | 3.880 | 6.810 | 2.949 | 600 | 1500 |
5.0 | 3.940 | 6.950 | 2.994 | 600 | 1500 |
5.0 | 4.050 | 6.980 | 3.004 | 600 | 1500 |
5.0 | 4.110 | 7.110 | 3.045 | 600 | 1500 |
5.0 | 4.180 | 7.410 | 3.141 | 600 | 1500 |
5.0 | 4.300 | 7.690 | 3.231 | 600 | 1500 |
10.0 | 4.390 | 8.050 | 3.346 | 600 | 1500 |
10.0 | 4.730 | 8.480 | 3.484 | 600 | 1500 |
0.0 | 4.780 | 8.480 | 3.484 | 600 | 1500 |
I have already created this .mod
file it should be in examples/pisco/pisco.mod
you need to copy it into home/project_name/structure
. Mudpy will look for .mod
files in this folder, you cannot put .mod
files anywhere else.
There is a function forward.makefault
that you can use to make a fault geometry (.fault) file in the format needed by mudpy. We will use the GlobalCMT moment tensor fault geometry and make a fault that has 10x10km subfaults. But first specify where you are going to save the fault geometry file:
>>> fout='/Users/dmelgar/Slip_inv/Pisco_insar/data/model_info/pisco_10km.fault'
The fault geometry file has to be in this folder mudpy will look for it here. Now make the file by typing
>>> from mudpy import forward
>>> forward.makefault(fout,strike=321,dip=28,nstrike=40,dx_dip=10,dx_strike=10,epicenter=[-76.509,-13.354,39],num_updip=7,num_downdip=4,rise_time=5)
You can make a quick plot of the location of the subfault centers to check that it covers the area youd esire by doing:
>>> from numpy import genfromtxt
>>> from matplotlib import pyplot as plt
>>> lonlat=genfromtxt(u'/Users/dmelgar/Slip_inv/Pisco_insar/data/model_info/pisco_10km.fault',usecols=[1,2])
>>> plt.scatter(lonlat[:,0],lonlat[:,1])
>>> plt.scatter(-76.509,-13.354,c='r',s=100)
>>> plt.axis('equal')
>>> plt.show()
This should make the following plot (the red circle is the hypocenter):
This fault model looks good, it has 40x12 subfaults (480 total) and it covers what we now know was the rupture area, but the earthquake was unilateral to the south so we might want to trim the northern faults, this also makes the inversion faster by reducing the number of parameters. Let's say we want to trim 18 columns of faults to the north, to do this we can do:
>>> from numpy import genfromtxt,savetxt,arange,r_
>>> fault=genfromtxt(u'/Users/dmelgar/Slip_inv/Pisco_insar/data/model_info/pisco_10km.fault')
>>> ndip=12
>>> nstrike=40
>>> good_faults=arange(0,nstrike-18) #We want to throw away the first 18 columns of subfaults
>>> keep_faults=good_faults.copy() #Initialize
>>> for k in range(ndip-1):
>>> keep_faults=r_[keep_faults,good_faults+nstrike*(k+1)]
>>> #Ok now throw away faults you don't want
>>> fault=fault[keep_faults,:]
>>> #Renumber the faults so the numbers are consecutive
>>> fault[:,0]=arange(1,len(fault)+1)
>>> #Save to file
>>> savetxt('/Users/dmelgar/Slip_inv/Pisco_insar/data/model_info/pisco_10km.fault',fault,fmt='%i\t%.6f\t%.6f\t%.3f\t%i\t%i\t%.1f\t%.1f\t%.2f\t%.2f')
Now the fault model is unilateral and has only 264 subfaults. You can plot it to check the geometry, again, the red circle is the hypocenter.
>>> from matplotlib import pyplot as plt
>>> fault=genfromtxt(u'/Users/dmelgar/Slip_inv/Pisco_insar/data/model_info/pisco_10km.fault')
>>> plt.scatter(fault[:,1],fault[:,2])
>>> plt.scatter(-76.509,-13.354,c='r',s=100)
>>> plt.axis('equal')
Again, you should confirm that the file is in $MUD/Pisco_insar/data/model_info/pisco_10km.fault
. Next we'll look at preparing the InSAR data.
Husen, S., Kissling, E., Flueh, E., & Asch, G. (1999). Accurate hypocentre determination in the seismogenic zone of the subducting Nazca Plate in northern Chile using a combined on-/offshore network. Geophysical Journal International, 138(3), 687-701.
Pritchard, M. E., E. O. Norabuena, C. Ji, R. Boroschek, D. Comte, M. Simons, T. Dixon, and P. A. Rosen (2007), Teleseismic, geodetic, and strong motion constraints on slip from recent southern Peru subduction zone earthquakes, J. Geophys. Res., 112, B03307