Calculating displacement due to faulting#
This notebook calculates displacement at specified points due to slip on a fault. The opening code cells load necessary files that define the problem: the fault geometry (in this case, the Cascadia subduction interface, discretized into triangular dislocation elements) and the location of GNSS stations where displacement can be measured.
The next code cells carry out a forward calculation of displacement due to fault slip. A radial Gaussian function, with variable position and gradient, defines the slip magnitude distribution, which is multiplied by the matrix of partial derivatives relating unit slip to station displacement.
The final code cells carry out an inversion of a displacement field (a slow slip event recorded at Cascadia GNSS stations). The estimation routine includes Laplacian regularization (smoothing), with a user-specified weight. Additional constraints imposed include a priori specification of slip magnitude and a bounded (non-negative) estimation.
# Install additional libraries
%pip install meshio
%pip install cutde
# Import libraries
import pandas as pd
import numpy as np # Numerical analysis
import scipy # Scientific analysis
import meshio # Interaction between fault model files and Python
import matplotlib.pyplot as plt # Visualize
from pyproj import Proj # Geographic coordinate transformation
import pyproj
import h5py # Efficient IO of arrays
import cutde.halfspace as cutde # Elastic dislocation code
# Make figures interactive
%matplotlib ipympl
0. Basic spherical geometry functions#
These functions convert between spherical and geocentric Cartesian coordinates, which we use to calculate some geometric parameters of the subduction interface geometry.
# Define some basic coordinate transformation functions
GEOID = pyproj.Geod(ellps="WGS84") # Define the geoid
KM2M = 1.0e3 # Conversion factor, km to m
RADIUS_EARTH = np.float64((GEOID.a + GEOID.b) / 2) # Average radius of Earth
# Convert spherical to geocentric Cartesian coordinates
def sph2cart(lon, lat, radius):
lon_rad = np.deg2rad(lon)
lat_rad = np.deg2rad(lat)
x = radius * np.cos(lat_rad) * np.cos(lon_rad)
y = radius * np.cos(lat_rad) * np.sin(lon_rad)
z = radius * np.sin(lat_rad)
return x, y, z
# Convert geocentric Cartesian to spherical coordinates
def cart2sph(x, y, z):
azimuth = np.arctan2(y, x)
elevation = np.arctan2(z, np.sqrt(x ** 2 + y ** 2))
r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
return azimuth, elevation, r
# Convert negative longitudes to 0–360º convention
def wrap2360(lon):
lon[np.where(lon < 0.0)] += 360.0
return lon
1. Reading in geometry files#
The problem geometry is defined by the fault geometry and the GNSS station locations. We are using triangular elements to represent the fault geometry, and in this case, we represent the geometry of the Cascadia subduction interface using a mesh of elements generated by the open-source tool Gmsh. The .msh file can be read using the library meshio. The GNSS information is stored in a .csv file containing station locations and displacements, and we can read that using pandas.
1.1. Read fault geometry file and parse triangular elements#
The fault geometry is defined using triangular dislocation elements (TDEs), which allow for approximation of geometrically complex surfaces.
The nodes of the mesh are defined as 3-dimensional x, y, z coordinates and can be specified for geographic meshes as longitude (º), latitude (º), depth (km). These nodes are stored as the m-by-3 field points of the mesh dictionary:
mesh["points"] = np.array([[1, 1, 0], [0.5, 0, 0], [0, 1, 0]])
The nodes are connected to form elements, with the entries of rows of the n-by-3 field verts of the mesh dictionary specifying the row indices of mesh["points"] that are connected to circulate the element in a counterclockwise manner. This circulation allows the normal vector (cross product of the TDE legs) to point upward for a dipping element.
mesh["verts"] = np.array([0, 2, 1])
In the drawing below, element 0 is defined by nodes 0, 2, and 1, and the nodal coordinates are as specified above.
2 ------- 0
\ /
\ 0 /
\ /
\ /
1
With this basic geometry, we can expand some geometric properties of the elements, identifying the coordinates of each element’s 3 nodes and, after converting the nodal coordinates to Cartesian coordinates, calculating the orientation (strike and dip).
# Read in source fault
filename = "../shared_data/cascadia.msh"
mesh = dict()
meshobj = meshio.read(filename)
mesh["file_name"] = filename
# The points field contains longitude, latitude, depth coordinates of the nodes of the triangles (m-by-3)
mesh["points"] = meshobj.points
# The verts field is a set of row indices into the points field, indicating which nodes comprise each triangle (n-by-3)
mesh["verts"] = meshio.CellBlock("triangle", meshobj.get_cells_type("triangle")).data
# Count of triangular elements
ntri = len(mesh["verts"])
# Expand mesh coordinates
mesh["lon1"] = mesh["points"][mesh["verts"][:, 0], 0]
mesh["lon2"] = mesh["points"][mesh["verts"][:, 1], 0]
mesh["lon3"] = mesh["points"][mesh["verts"][:, 2], 0]
mesh["lat1"] = mesh["points"][mesh["verts"][:, 0], 1]
mesh["lat2"] = mesh["points"][mesh["verts"][:, 1], 1]
mesh["lat3"] = mesh["points"][mesh["verts"][:, 2], 1]
mesh["dep1"] = mesh["points"][mesh["verts"][:, 0], 2]
mesh["dep2"] = mesh["points"][mesh["verts"][:, 1], 2]
mesh["dep3"] = mesh["points"][mesh["verts"][:, 2], 2]
# Calculate centroids
mesh["centroids"] = np.mean(mesh["points"][mesh["verts"], :], axis=1)
# Convert nodes to Cartesian coordinates in meters
mesh["x1"], mesh["y1"], mesh["z1"] = sph2cart(
mesh["lon1"],
mesh["lat1"],
RADIUS_EARTH + KM2M * mesh["dep1"],
)
mesh["x2"], mesh["y2"], mesh["z2"] = sph2cart(
mesh["lon2"],
mesh["lat2"],
RADIUS_EARTH + KM2M * mesh["dep2"],
)
mesh["x3"], mesh["y3"], mesh["z3"] = sph2cart(
mesh["lon3"],
mesh["lat3"],
RADIUS_EARTH + KM2M * mesh["dep3"],
)
# Cartesian triangle centroids
mesh["x_centroid"] = (mesh["x1"] + mesh["x2"] + mesh["x3"]) / 3.0
mesh["y_centroid"] = (mesh["y1"] + mesh["y2"] + mesh["y3"]) / 3.0
mesh["z_centroid"] = (mesh["z1"] + mesh["z2"] + mesh["z3"]) / 3.0
# Calculate cross products, which give element orientations
tri_leg1 = np.transpose([np.deg2rad(mesh["lon2"] - mesh["lon1"]), np.deg2rad(mesh["lat2"] - mesh["lat1"]), (1 + KM2M * mesh["dep2"] / RADIUS_EARTH) - (1 + KM2M * mesh["dep1"] / RADIUS_EARTH)])
tri_leg2 = np.transpose([np.deg2rad(mesh["lon3"] - mesh["lon1"]), np.deg2rad(mesh["lat3"] - mesh["lat1"]), (1 + KM2M * mesh["dep3"] / RADIUS_EARTH) - (1 + KM2M * mesh["dep1"] / RADIUS_EARTH)])
# Normal vector for each element
mesh["nv"] = np.cross(tri_leg1, tri_leg2)
# Convert Cartesian normal vector to angular units
azimuth, elevation, r = cart2sph(mesh["nv"][:, 0], mesh["nv"][:, 1], mesh["nv"][:, 2])
# Extract geometric quantities for each element
mesh["strike"] = wrap2360(-np.rad2deg(azimuth))
mesh["dip"] = 90 - np.rad2deg(elevation)
mesh["dip_flag"] = mesh["dip"] != 90
1.2. Read GNSS data#
We’re reading in an example displacement field from a slow slip event in Cascadia, using pandas to generate a DataFrame. The columns of the .csv file and resulting DataFrame are:
name, lon, lat, east_disp, north_disp, up_disp, east_sig, north_sig, up_sig
gnss = pd.read_csv('../shared_data/casc_sse_07-2011.csv')
nsta = len(gnss)
lon_corr = 1
# Check longitude convention of mesh
# If longitudes are already 0–360º, we don't need to convert
if np.max(gnss.lon) > 180:
lon_corr = 0
# Reshape grid into a 3 column array
pts = np.array([gnss.lon, gnss.lat, 0 * gnss.lat]).reshape((3, -1)).T.copy()
1.3. Plot fault and stations#
To visualize the problem geometry, we can plot both the fault elements and GNSS station locations, along with a coastline for geographic reference.
# Setting up axis limits around fault
lonmin = np.mean(mesh["points"][:, 0]) - 3*np.std(mesh["points"][:, 0])
lonmax = np.mean(mesh["points"][:, 0]) + 3*np.std(mesh["points"][:, 0])
latmin = np.mean(mesh["points"][:, 1]) - 3*np.std(mesh["points"][:, 1])
latmax = np.mean(mesh["points"][:, 1]) + 3*np.std(mesh["points"][:, 1])
# Read in coastline file
coast = pd.read_csv("../shared_data/coastline.csv")
# Define the figure and its axis
fig, ax = plt.subplots()
# Draw the fault
ax.triplot(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"], linewidth=0.5)
# Add grid points
ax.plot(pts[:, 0], pts[:, 1], '.r')
# Add coastline
ax.plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
ax.set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
plt.show()
1.4. Coordinate transformation#
Calculating the relationship between fault slip and displacement is done in Cartesian coordinates, so we’ll need to convert our longitude, latitude coordinates for both the fault and the GNSS stations. The approach below uses a UTM projection, which is generally a good choice for regions that are aligned in a north-south direction, such as Cascadia. Other map projections would work, too.
# Convert the fault and station coordinates to UTM coordinates
# Find UTM zone. UTM zones are 6º wide, and this formula finds the EPSG code for the zone
# WGS 84 zones' EPSG codes start at 32600, and we add the zone to that. Add 100 if in the southern hemisphere.
utmzone=int(32600 + np.floor((-180 + np.mean(gnss.lon))/6))
if np.mean(gnss.lat) < 0:
utmzone += 100
# Define target UTM coordinate system
target_crs = 'epsg:'+str(utmzone)
# Define original source geographic coordinate system
source_crs = 'epsg:4326' # Global lat-lon coordinate system
# Define the coordinate transformation using pyproj
latlon_to_utm = pyproj.Transformer.from_crs(source_crs, target_crs)
# Convert coordinates of the fault and the stations
# We're just converting the longitude, latitude coordinates of the fault, as the depths are already in linear units
# Dividing by 1e3 expresses UTM coordinates in km
faultxy = np.array(latlon_to_utm.transform(mesh["points"][:, 1], mesh["points"][:, 0])).T/1e3
gnssxy = np.array(latlon_to_utm.transform(pts[:, 1], pts[:, 0])).T/1e3
# Place UTM coordinates into 3-column arrays (x, y, z)
# Allocate space for fault coordinates
cart_fault_pts = np.zeros_like(mesh["points"])
# Insert x, y
cart_fault_pts[:, 0:2] = faultxy
# Insert depths as the third column
cart_fault_pts[:, 2] = mesh["points"][:, 2]
# Allocate space for grid coordinates
cart_gnss_pts = np.zeros_like(pts)
# Insert x, y
cart_gnss_pts[:, 0:2] = gnssxy
# Insert z (all zero, at the surface of the half-space)
cart_gnss_pts[:, 2] = pts[:, 2]
# Make a plot of the fault and stations, both in geographic and UTM coordinates
fig, ax = plt.subplots(1, 2)
# Draw the fault
ax[0].triplot(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"], linewidth=0.5)
# Add grid points
ax[0].plot(pts[:, 0], pts[:, 1], '.r')
# Add coastline
ax[0].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
ax[0].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[0].set_xlabel("Longitude")
ax[0].set_ylabel("Latitude")
ax[0].set_title("Geographic")
ax[1].triplot(cart_fault_pts[:,0], cart_fault_pts[:, 1], mesh["verts"], linewidth=0.5)
ax[1].plot(cart_gnss_pts[:, 0], cart_gnss_pts[:, 1], '.r')
ax[1].set_aspect("equal")
ax[1].set_xlabel("Eastings (km)")
ax[1].set_ylabel("Northings (km)")
ax[1].set_title("UTM")
plt.show()
1.5. Elastic dislocation partial derivatives#
The Python library cutde calculates the partial deriviatives relating unit slip on the TDEs to displacement at GNSS station locations.
Each TDE can have 3 slip (dislocation) components, expressed in an element-local coordinate system. The first component is parallel to the element strike, the second is parallel to the element dip, and the third is normal to the element (tensile slip). For fault slip problems, we typically ignore the third component, focusing on the element-parallel slip components.
Unit slip in each component direction yields displacement at all GNSS stations in 3 dimensions, expressed in a global x, y, z coordinate system; we make an assumption here that the UTM (Cartesian) coordinate system describes displacement in the east, north, and up directions.
Upon calculating the partial derivative matrix, we reshape it to shape 3 n-by-3 m, where n is the number of stations and m is the number of TDEs. The structure of this matrix is:
where \(g^{S_0}_{E_0}\) is the Green’s function, \(g\), that gives displacement in the east (\(E\)) direction at station 0 due to unit slip in the strike (\(S\)) direction on TDE 0. This matrix can then be multiplied by a slip vector of size 3 m-by-1 to give the total displacement at all stations due to slip on all TDEs.
# Calculate the mathematical relationship between unit fault slip and displacement
# This can take a bit of time, but it only needs to be done once if your fault geometry and observation coordinates don't change
disp_mat = cutde.disp_matrix(obs_pts=cart_gnss_pts, tris=cart_fault_pts[mesh["verts"]], nu=0.25).reshape((-1, 3*ntri))
# If running on your own computer, uncomment the following lines to write to an HDF5 file
# This file could be reloaded from disk (see next code cell; all lines intentionall commented) rather than recalculating the matrix, if your problem geometry remains the same
# hf = h5py.File("disp_mat.h5", "w")
# hf["disp_mat"] = disp_mat
# hf.close()
# Optional: If you've already calculated the slip-to-displaceent matrix, read disp_mat from file
# hf = h5py.File('disp_mat.h5', 'r')
# disp_mat = np.array(hf["disp_mat"])
# hf.close()
2. Relating a specific fault slip distribution to a surface displacement field#
Using the partial deriviates relating unit slip on each element to displacement at the surface observation coordinates, we can now determine the total displacement at GNSS stations that results from a defined slip distribution on the fault. We can do this by integrating the partial derivatives across the domain of the fault, and in practice we do this through matrix multiplication: matrix multiplication of the partial derivatives times the slip vector yields the total displacement arising from the contributions of all elements’ slip magnitudes.
To get a sense for how a given slip distribution impacts the surface displacement field, we can prescribe a slip distribution and then predict the displacement at GNSS stations. We can propose a slip distribution defined by a 2-dimensional Gaussian distribution, smoothly decaying to zero away from a peak value. We define the center of the distribution, the “steepness” of the decay in the x and y directions, and a rotation of the y-axis. These specifications allow us to specify the distribution in a way that may vary the dimensions along the strike and dip directions of the source fault.
2.1. Forward model slip function#
This function produces a slip value for each TDE on the fault geometry, based on a 2-D Gaussian distribution. A true radial Gaussian function would produce a “bullseye” distribution of slip magnitudes, but this function can have a different gradient in the x and y directions, as well as a rotation applied to the axes, allowing for an elliptical shape.
The function only calculates the slip magnitude. Partitioning of this magnitude into slip in the element strike and dip directions would require some trigonometric modification.
def radgauss(x0, y0, a, x, y, c):
#
# RADGAUSS calculates a radial Gaussian distribution.
#
# S = RADGAUSS(X0, Y0, A, X, Y, C) uses the formula:
#
# S(Xi, Yi) = A exp(-Ri^2/2C^2)
#
# to calculate the value S at point (Xi,Yi), which is located a distance
# Ri from point (X0,Y0), where the value is A. C represents a Gaussian
# curvature parameter: larger values of C produce more abrupt changes in
# value with distance.
#
# X0, Y0 represent the x,y coordinates of the point of maximum magnitude A.
# X and Y are n x 1 arrays containing the x,y coordinates of all points at
# which the value S should be calculated. C is a single value that applies
# to the entire calculation.
#
# If C is a 2-element vector, the pattern will not be radial but elliptical,
# with the first element of C giving the curvature in the X direction and the
# second giving the curvature in Y:
#
# S(Xi, Yi) = A exp( [ -(Xi-X0)^2 / 2C[0]^2 ] + [ -(Yi-Y0)^2 / 2C[1]^2] )
#
# If C is a 3-element vector, an elliptical pattern will result, and the
# ellipse will be rotated C[2] degrees from the y-axis, where a positive
# number gives a clockwise rotation.
c = np.array(c)
if np.size(c) == 1:
r = np.sqrt((x-x0)**2+(y-y0)**2) # calculate distances
s = a*np.exp(-r**2/(2*c**2)) # calculate Gaussian value
elif np.size(c) == 2:
s = a*np.exp((-(x-x0)**2./(2*c[0]**2)) + (-(y-y0)**2./(2*c[1]**2)))
elif np.size(c) == 3:
ang = np.deg2rad(c[2])
rot = np.array([[np.cos(ang), -np.sin(ang)], [np.sin(ang), np.cos(ang)]])
coords = np.vstack((np.hstack((x.T, x0)), np.hstack((y.T, y0))))
rcoords = rot.dot(coords)
s = a*np.exp((-(rcoords[0, :-1]-rcoords[0, -1])**2/(2*c[0]**2)) + (-(rcoords[1, :-1]-rcoords[1, -1])**2/(2*c[1]**2)))
return s
2.2. Specify slip distribution and predict displacements#
This code block can be run repeatedly to generate a trial-and-error prediction of the observed GNSS displacement field. The parameters of the radial Gaussian slip distribution can be edited in the opening lines, and then the rest of the block will calculate the predicted displacement arising from that specified distribution. This is commonly expressed as the algebraic equation \(Gm = d\), where \(G\) is the matrix of partial derivatives relating unit slip to displacement; \(m\) is a set of model parameters, which in this case is a vector of slip component magnitudes; and \(d\) is a data vector, which here is a set of displacements at GNSS stations. Finally, the code generates a figure that shows the spatial distribution of dip slip and the predicted and observed displacement field.
# Call radial slip distribution function
# Center coordinates of radial slip distribution
lon_center = 236 #<-- MODIFY
lat_center = 45 #<-- MODIFY
# 3 parameters: slip gradient in x, y directions plus rotation angle
# Larger gradient value = more smeared out
slip_shape = np.array([1, 0.5, 0]) #<-- MODIFY
# Calculate Gaussian slip distribution with selected parameters
s = radgauss(lon_center, lat_center, 1, mesh["centroids"][:, 0], mesh["centroids"][:, 1], slip_shape)
# Create full slip distribution, adding radial slip as dip-slip
slip = np.vstack((np.zeros_like(s), 10*s, np.zeros_like(s))).T
# Calculate predicted displacement at grid points that arises from fault slip
pred_disp = disp_mat.dot(slip.reshape(3*ntri, 1))
plt.close('all')
# Visualize the slip and displacement pattern
fig, ax = plt.subplots(1, 2)
fso = ax[0].tripcolor(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"], facecolors=slip[:, 1])
# fig.colorbar(fso, ax=ax[0], orientation='horizontal')
ax[0].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
ax[0].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[0].title.set_text("Fault slip")
# Observed displacements
ax[1].quiver(gnss.lon, gnss.lat, gnss.east_disp, gnss.north_disp, color='k', scale=50, label="Observed")
# Predicted displacements
ax[1].quiver(gnss.lon, gnss.lat, pred_disp[0::3], pred_disp[1::3], color='r', scale=50, label="Modeled")
# Add coastline
ax[1].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
ax[1].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[1].legend()
ax[1].yaxis.tick_right()
ax[1].title.set_text("Displacement")
# Add colorbar
cbar_ax = fig.add_axes([0.45, 0.15, 0.025, 0.7])
fig.colorbar(fso, cax=cbar_ax, label="Slip (mm)")
# Add vector scale
ax[1].quiver(235, 37.3, -5, 0, color='k', scale=50)
ax[1].text(233.7, 37.8, '5 mm')
plt.show()
3. Inversion of displacement data for slip#
Rather than specifying a slip distribution a priori, we can estimate a slip distribution, using linear least-squares fitting.
3.1. Naive estimate of slip#
In simple algebraic terms, if \(Gm = d\), then \(m = G^{-1}d\). In other words, we could invert the matrix of partial derivatives and multiply by the displacement vector to estimate the model parameters (slip values). \(G\) is not a square matrix; there are more TDE slip components than there are GNSS displacement components, and so there are more columns than there are rows. We can instead estimate a set of slip values as:
\(m = (G^T G)^{-1}G^T d\),
where \((G^T G)^{-1}\) is the model covariance.
If we consider a set of data weights, \(W_d\), given as a diagonal matrix with entries as the inverse square of data uncertainties, \(1/\sigma^2\), we could carry out a weighted estimation of slip values as:
\(m = (G^T W_d G)^{-1}G^T W_d d\).
Before carrying out this estimation, we will trim components from \(G\), removing every row coresponding to vertical displacement (the vertical displacements are very noisy) and every column corresponding to tensile dislocation, assuming that the observed displacements arise from shear dislocation on TDEs.
# Simple inversion
# From the matrix of partial derivatives,
# Delete ROWS corresponding to vertical displacement
disp_mat_no_vert = np.delete(disp_mat, np.arange(2, 3*nsta, 3), axis=0)
# Delete COLUMNS corresponding to tensile slip
disp_mat_final = np.delete(disp_mat_no_vert, np.arange(2, 3*ntri, 3), axis=1)
# Weights on GNSS displacements are 1/sigma^2
data_uncertainties = np.array([gnss.east_sig, gnss.north_sig]).T.reshape((2*nsta, -1)).copy()
data_weights = 1./data_uncertainties**2
# Displacement vector
disp_vector = np.array([gnss.east_disp, gnss.north_disp]).T.reshape((2*nsta, -1)).copy()
# Calculate model covariance
cov = np.linalg.inv(disp_mat_final.T @ disp_mat_final)
# Estimate slip using pre-calculated covariance
est_slip = cov @ disp_mat_final.T @ disp_vector
# Predict displacement at stations
pred_sse_disp = disp_mat_final.dot(est_slip)
# Visualize the slip and displacement pattern
fig, ax = plt.subplots(1, 2)
fso = ax[0].tripcolor(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"], facecolors=est_slip[1::2, 0])
ax[0].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
ax[0].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[0].title.set_text("Fault slip")
# Observed displacements (columns of the GNSS DataFrame)
ax[1].quiver(gnss.lon, gnss.lat, gnss.east_disp, gnss.north_disp, color='k', scale=50, label="Observed")
# Predicted displacements (alternating rows of the pred_sse_disp array)
ax[1].quiver(gnss.lon, gnss.lat, pred_sse_disp[0::2], pred_sse_disp[1::2], color='r', scale=50, label="Modeled")
# Add coastline
ax[1].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
# Axis limits
ax[1].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[1].legend()
ax[1].yaxis.tick_right()
ax[1].title.set_text("Displacement")
# Add colorbar
cbar_ax = fig.add_axes([0.45, 0.15, 0.025, 0.7])
fig.colorbar(fso, cax=cbar_ax, label="Slip (mm)")
# Add vector scale
ax[1].quiver(235, 37.3, -5, 0, color='k', scale=50)
ax[1].text(233.7, 37.8, '5 mm')
plt.show()
3.2. Regularized estimation of slip#
We can see above that the slip distribution estimated from our simple psuedoinverse is physically unreasonable and yields a very poor fit to the observed SSE displacements. One reason for this is that the problem is underdetermined, that is, there are more model parameters (slip components) to estimate than there are constraining data (GNSS displacement components). To generate a more physically reasonable slip distribution, we can regularize the inversion. There are many approaches to regularization (as well as many approaches to model parameter estimation without regularization), and we will use a simple method of Laplacian smoothing. In carrying out this regularized estimation, we impose a new constraint on the problem. Not only should the estimated slip distribution predict displacements that match the observations, but the slip distribution should also vary smoothly in space. A potential physical justification for this is that, if a slip event features a uniform stress drop across the slipped area of the fault, the corresponding slip distribution will smoothly vary in magnitude.
The Laplacian smoothing operator, \(L\) serves to take a discrete second derivative of the slip distribution. By setting the product \(Lm = 0\), we enforce minimization of the curvature in the spatial distribution of slip magnitude. The smoothest possible slip distribution is one with homogeneous magnitudes, as it has zero curvature.
Constructing the square matrix \(L\) involves setting the diagonal entries equal to 3 and, in a given row corresponding to a slip component on a target element \(i\), setting the columns corresponding to slip components of the elements neighboring \(i\) to -1. These values can be visualized for a given element in the code cell below.
How well the estimated slip distribution fits the true data (GNSS displacements) versus the pseudodata (smoothing criterion) can be controlled through weighting, similar to what we did above with the simple inversion. We can consider the relative balance in the weights on the data and pseudodata by setting the data weights, \(W_d\), to be the inverse squared uncertainties as above while varying a weighting factor applied to the pseudodata, \(W_L\). The higher the chosen pseudodata weight, the more the smoothing criterion influences the estimated slip.
3.2.1. Constructing the smoothing operator#
The smoothing operator is constructed first by finding the indices of the (up to) 3 elements that share a side with each target element. These indices can be used to construct the smoothing operator.
Below, we find the shared sides and construct the smoothing operator, then show that operator for a single element to highlight the construction.
# We first need to determine the neighbors of each element, defined as those that share a side
# From https://github.com/brendanjmeade/celeri/blob/main/celeri/spatial.py
def get_shared_sides(vertices):
"""Determine the indices of the triangular elements sharing
one side with a particular element.
Inputs:
vertices: n x 3 array containing the 3 vertex indices of the n elements,
assumes that values increase monotonically from 1:n.
Outputs:
share: n x 3 array containing the indices of the m elements sharing a
side with each of the n elements. "-1" values in the array
indicate elements with fewer than m neighbors (i.e., on
the edge of the geometry).
In general, elements will have 1 (mesh corners), 2 (mesh edges), or 3
(mesh interiors) neighbors, but in the case of branching faults that
have been adjusted with mergepatches, it's for edges and corners to
also up to 3 neighbors.
"""
# Make side arrays containing vertex indices of sides
side_1 = np.sort(np.vstack((vertices[:, 0], vertices[:, 1])).T, 1)
side_2 = np.sort(np.vstack((vertices[:, 1], vertices[:, 2])).T, 1)
side_3 = np.sort(np.vstack((vertices[:, 0], vertices[:, 2])).T, 1)
sides_all = np.vstack((side_1, side_2, side_3))
# Find the unique sides - each side can part of at most 2 elements
_, first_occurence_idx = np.unique(sides_all, return_index=True, axis=0)
_, last_occurence_idx = np.unique(np.flipud(sides_all), return_index=True, axis=0)
last_occurence_idx = sides_all.shape[0] - last_occurence_idx - 1
# Shared sides are those whose first and last indices are not equal
shared = np.where((last_occurence_idx - first_occurence_idx) != 0)[0]
# These are the indices of the shared sides
sside1 = first_occurence_idx[shared]
sside2 = last_occurence_idx[shared]
el1, sh1 = np.unravel_index(
sside1, vertices.shape, order="F"
) # "F" is for fortran ordering.
el2, sh2 = np.unravel_index(sside2, vertices.shape, order="F")
share = -1 * np.ones((vertices.shape[0], 3))
for i in range(el1.size):
share[el1[i], sh1[i]] = el2[i]
share[el2[i], sh2[i]] = el1[i]
share = share.astype(int)
return share
# Using the indices of the shared sides, we can construct the smoothing matrix
# From https://github.com/brendanjmeade/celeri/blob/main/celeri/spatial.py
def get_tri_smoothing_matrix_simple(share, n_dim):
"""Produces a smoothing matrix based without scale-dependent
weighting.
Inputs:
share: n x 3 array of indices of the up to 3 elements sharing a side with each of the n elements
n_dim: Number of dimensions (slip components) for each element
Outputs:
smoothing matrix: n_dim * n x n_dim * n smoothing matrix
"""
# Allocate space for contructing smoothing matrix
# For large problems, this can be made as a sparse matrix
n_shared_tri = share.shape[0]
smoothing_matrix = np.zeros((n_dim * n_shared_tri, n_dim * n_shared_tri))
for j in range(n_dim):
for i in range(n_shared_tri):
smoothing_matrix[n_dim * i + j, n_dim * i + j] = 3
if share[i, j] != -1:
k = n_dim * i + np.arange(n_dim)
m = n_dim * share[i, j] + np.arange(n_dim)
smoothing_matrix[k, m] = -1
return smoothing_matrix
# Call functions to create smoothing matrix
share = get_shared_sides(mesh["verts"])
smoothing_mat = get_tri_smoothing_matrix_simple(share, n_dim=3)
# Trim rows and columns to remove tensile slip components
smoothing_mat = np.delete(smoothing_mat, np.arange(2, 3*ntri, 3), axis=0)
smoothing_mat = np.delete(smoothing_mat, np.arange(2, 3*ntri, 3), axis=1)
# Figure showing a single row of the smoothing operator, corresponding to a single TDE
fig, ax = plt.subplots(1, 2)
this_el = 102 #<-- MODIFY
fso = ax[0].tripcolor(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"], facecolors=smoothing_mat[2*this_el, 0::2])
ax[0].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
ax[0].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[0].title.set_text("Smoothing operator")
# Zoom in on target element
zoom_inc = 0.5
lonmin_el = mesh["centroids"][this_el, 0] - zoom_inc
lonmax_el = mesh["centroids"][this_el, 0] + zoom_inc
latmin_el = mesh["centroids"][this_el, 1] - zoom_inc
latmax_el = mesh["centroids"][this_el, 1] + zoom_inc
ax[1].tripcolor(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"], color="white", linewidth=0.1, facecolors=smoothing_mat[2*this_el, 0::2])
ax[1].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
# Axis limits
ax[1].set(xlim=(lonmin_el, lonmax_el), ylim=(latmin_el, latmax_el), aspect='equal')
ax[1].yaxis.tick_right()
ax[1].title.set_text(f"Element {this_el}")
ax[0].plot([lonmin_el, lonmax_el, lonmax_el, lonmin_el, lonmin_el], [latmin_el, latmin_el, latmax_el, latmax_el, latmin_el], 'r')
# Add colorbar
cbar_ax = fig.add_axes([0.575, 0.15, 0.3, 0.025])
fig.colorbar(fso, cax=cbar_ax, orientation="horizontal", label="Operator value")
plt.show()
3.3. Assemble matrices, weighting matrix, and data vector#
To estimate slip constrained by the displacement data and Laplacian smoothing, we can use a modified version of the equation in Section 3.1:
\(m = \left( \begin{bmatrix}G \\ L\end{bmatrix}^T \begin{bmatrix}W_d \\ W_L \end{bmatrix} \begin{bmatrix}G \\ L\end{bmatrix}\right)^{-1}\begin{bmatrix}G \\ L\end{bmatrix}^T \begin{bmatrix}W_d \\ W_L \end{bmatrix} \begin{bmatrix}d \\ \mathbf{0} \end{bmatrix}\).
The top row of the block matrices is an estimate of slip aiming to satisfy \(Gm=d\), while the bottom row aims to meet the smoothing criterion, \(Lm=0\). By solving the simultaneous system of equations, we find a slip distribution that is a compromise fit to both constraints, with the goodness-of-fit to each constraint governed by the chosen weighting values, \(W_d\) and \(W_L\). In practice, we can just vary \(W_L\) and assess how it impacts the resulting slip distribution and fit to the displacement data.
In the code cell below, we construct the block matrices by vertically stacking their components.
# Assemble matrices:
assembled_mat = np.vstack([disp_mat_final, smoothing_mat])
# List of smoothing weights to test
smoothing_coefficients = [1e-2, 1e0, 1e2] # <-- MODIFY
# Select a smoothing weight and apply to meshes
smoothing_coefficient = smoothing_coefficients[1] # <-- MODIFY
# Assemble weights:
# Weights on pseudodata controlling regularization are applied to all slip rate rows
smoothing_weights = smoothing_coefficient*np.ones((len(smoothing_mat), 1))
# Complete weighting vector
all_weights = np.vstack([data_weights, smoothing_weights])
# Assemble data:
# Displacement vector
disp_vector = np.array([gnss.east_disp, gnss.north_disp]).T.reshape((2*nsta, -1)).copy()
# Psuedodata vector
zero_vector = np.zeros_like(smoothing_weights)
# Complete data vector
all_data = np.vstack([disp_vector, zero_vector])
3.4. Carry out smoothed inversion#
With the matrices assembled, we calculate the model covariance, then use it to estimate the slip. As before, we can then multiply the estimated slip vector by the Green’s functions to yield predicted displacements. We then can plot the slip distribution and compare the predicted vs. observed displacements to gauge the goodness-of-fit.
# Calculate model covariance
cov = np.linalg.inv(assembled_mat.T * all_weights.T @ assembled_mat)
# Estimate slip using pre-calculated covariance
est_slip = cov @ assembled_mat.T *all_weights.T @ all_data
# Predict displacement at stations
pred_sse_disp = disp_mat_final.dot(est_slip)
# Visualize the slip and displacement pattern
fig, ax = plt.subplots(1, 2)
fso = ax[0].tripcolor(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"], facecolors=est_slip[1::2, 0])
ax[0].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
ax[0].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[0].title.set_text("Fault slip")
# Observed displacements
ax[1].quiver(gnss.lon, gnss.lat, gnss.east_disp, gnss.north_disp, color='k', scale=50, label="Observed")
# Predicted displacements
ax[1].quiver(gnss.lon, gnss.lat, pred_sse_disp[0::2], pred_sse_disp[1::2], color='r', scale=50, label="Modeled")
# Add coastline
ax[1].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
# Axis limits
ax[1].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[1].legend()
ax[1].yaxis.tick_right()
ax[1].title.set_text("Displacement")
# Add colorbar
cbar_ax = fig.add_axes([0.45, 0.15, 0.025, 0.7])
fig.colorbar(fso, cax=cbar_ax, label="Slip (mm)")
# Add vector scale
ax[1].quiver(235, 37.3, -5, 0, color='k', scale=50)
ax[1].text(233.7, 37.8, '5 mm')
plt.show()
3.5. Independent exercise: Assess models with a variety of smoothing weights#
To get a sense for how the chosen smoothing weight impacts the overall model, there are several metrics that might be worth exploring. We want the estimated slip distribution to predict a good fit to the GNSS displacements, so evaluating the residual displacements is important. We might also want to characterize how physically reasonable the fault slip distribution is, gauging features such as the seismic moment of the estimated distribution, the minimum/maximum magnitudes, and rake distribution.
Add an additional code cell that allows you to run several models with varying smoothing weights, assessing some metrics of model quality for each test.
3.6. Additional slip distribution constraints#
In addition to enforcing a spatially smooth slip distribution, we can add more constraints. For example, we may want to force the slip to go to zero around the edges of the fault geometry, which would prevent a significant discontinuity in slip (large magnitude to zero slip, which is implicit when we are off the fault). We may also want to force the sign of slip (dip-slip in particular) to be positive (reverse), avoiding estimation of normal slip that is more consistent with coupling and strain accumulation on the fault, rather than the strain release that is achieved by slow slip events.
3.6.1. A priori setting of slip values#
To force slip on a subset of TDEs to be a particular value, we can augment our system of equations from Section 3.3 with some additional, simple rows:
\(m = \left( \begin{bmatrix}G \\ L \\ \mathbf{I} \end{bmatrix}^T \begin{bmatrix}W_d \\ W_L \\ W_C \end{bmatrix} \begin{bmatrix}G \\ L \\ \mathbf{I}\end{bmatrix}\right)^{-1}\begin{bmatrix}G \\ L \\ \mathbf{I}\end{bmatrix}^T \begin{bmatrix}W_d \\ W_L \\ W_C \end{bmatrix} \begin{bmatrix}d \\ \mathbf{0} \\ \mathbf{s} \end{bmatrix}\),
where \(\mathbf{I}\) represents a sparse matrix of size c-by-3 m with values of 1 in the columns corresponding to the slip components whose magnitudes are being constrained (with c as the total number of constraints), \(W_C\) is a set of weights applied to the a priori slip constraints, and \(\mathbf{s}\) are the a priori slip rates themselves, usually zeros for enforicing a no-slip constraint.
In the code cell below, we will set a zero-slip constraint on the bottom edge of fault, encouraging the slip to taper as the edge of the geometry is approached.
# Find indices of bottom elements
# Maximum depth of any node used by an element
max_depth = np.min(mesh["points"][mesh["verts"], 2])
# True/false array indicating element nodes at max depth
check_node_depth = np.column_stack((mesh["dep1"] == max_depth, mesh["dep2"] == max_depth, mesh["dep3"] == max_depth))
# Bottom edge eleements are those with two of their nodes at the max depth
bottom_el_index = np.where(np.sum(check_node_depth, axis=1) == 2)[0]
# Construct constraint array
# We're setting both strike and dip slip components to zero, so we need two rows for every bottom edge element
slip_constraint_mat = np.zeros((2*len(bottom_el_index), 2*ntri))
for i in range(len(bottom_el_index)):
# Set strike slip entry to 1
slip_constraint_mat[2*i, 2*bottom_el_index[i]] = 1
# Set dip slip entry to 1
slip_constraint_mat[2*i+1, 2*bottom_el_index[i]+1] = 1
# Add this array to the bottom of the combined array from before
# Assemble matrices:
assembled_mat_constrained = np.vstack([assembled_mat, slip_constraint_mat])
# Assemble weights:
# Weights on slip constraints
slip_constraint_weight = 1e15 # A big number to strongly enforce the constraint
slip_constraint_weights = slip_constraint_weight*np.ones((len(slip_constraint_mat), 1))
# Complete weighting vector
all_weights_constrained = np.vstack([all_weights, slip_constraint_weights])
# Assemble data:
# Slip constraint vector
slip_constraint_vector = np.zeros_like(slip_constraint_weights)
# Complete data vector
all_data_constrained = np.vstack([all_data, slip_constraint_vector])
# Calculate model covariance
cov_constrained = np.linalg.inv(assembled_mat_constrained.T * all_weights_constrained.T @ assembled_mat_constrained)
# Estimate slip using pre-calculated covariance
est_slip_constrained = cov_constrained @ assembled_mat_constrained.T *all_weights_constrained.T @ all_data_constrained
# Predict displacement at stations
pred_sse_disp_constrained = disp_mat_final.dot(est_slip_constrained)
# Visualize the slip and displacement pattern
fig, ax = plt.subplots(1, 2)
fso = ax[0].tripcolor(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"], facecolors=est_slip_constrained[1::2, 0])
# Highlight constrained elements
ax[0].triplot(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"][bottom_el_index, :], color='r', linewidth=0.5)
ax[0].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
ax[0].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[0].title.set_text("Fault slip")
# Observed displacements
ax[1].quiver(gnss.lon, gnss.lat, gnss.east_disp, gnss.north_disp, color='k', scale=50, label="Observed")
# Predicted displacements
ax[1].quiver(gnss.lon, gnss.lat, pred_sse_disp_constrained[0::2], pred_sse_disp_constrained[1::2], color='r', scale=50, label="Modeled")
# Add coastline
ax[1].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
# Axis limits
ax[1].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[1].legend()
ax[1].yaxis.tick_right()
ax[1].title.set_text("Displacement")
# Add colorbar
cbar_ax = fig.add_axes([0.45, 0.15, 0.025, 0.7])
fig.colorbar(fso, cax=cbar_ax, label="Slip (mm)")
# Add vector scale
ax[1].quiver(235, 37.3, -5, 0, color='k', scale=50)
ax[1].text(233.7, 37.8, '5 mm')
plt.show()
3.6.2. Bounded slip values#
Depending on the smoothing weight set, the above estimated slip distributions feature some patches of negative (normal-sense) dip-slip. If our expectation for a slow slip event is that all slip is positive (reverse-sense), we can impose additional constraints on the estimated slip by way of a bounded optimization. The scipy package’s optimization routines allow for straightforward specification of bounds to be placed on the estimated model parameters.
The slow slip event we’re using as an example is located near the prominent bend in the Cascadia subduction zone. Therefore, we may expect a nominally reverse-sense slow slip event to feature some components of strike-slip that vary in sigh. Therefore, we may wish to place wide bounds on the sign of the estimated strike-slip while forcing the dip slip values to be positive.
# Initialize slip bounds array
# Two columns (min and max), with one row for each slip component
min_bounds = np.zeros((len(est_slip), ))
max_bounds = np.zeros((len(est_slip), ))
min_bounds[0::2] = -100
max_bounds[0::2] = 100
min_bounds[1::2] = 0
max_bounds[1::2] = 100
# Using the constrained problem above, use scipy.optimize.lsq_linear
# Create weighted versions of the matrix and data vector
weighted_assembled_mat_constrained = np.sqrt(all_weights_constrained)*assembled_mat_constrained
weighted_all_data_constrained = all_data_constrained*np.sqrt(all_weights_constrained)
est_slip_bounded = scipy.optimize.lsq_linear(weighted_assembled_mat_constrained, weighted_all_data_constrained[:, 0], bounds=(min_bounds, max_bounds), lsmr_tol='auto')
# Predict displacement at stations
pred_sse_disp_bounded = disp_mat_final.dot(est_slip_bounded.x)
# Visualize the slip and displacement pattern
fig, ax = plt.subplots(1, 2)
fso = ax[0].tripcolor(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"], facecolors=est_slip_bounded.x[1::2])
# Highlight constrained elements
ax[0].triplot(mesh["points"][:, 0], mesh["points"][:, 1], mesh["verts"][bottom_el_index, :], color='r', linewidth=0.5)
ax[0].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
ax[0].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[0].title.set_text("Fault slip")
# Observed displacements
ax[1].quiver(gnss.lon, gnss.lat, gnss.east_disp, gnss.north_disp, color='k', scale=50, label="Observed")
# Predicted displacements
ax[1].quiver(gnss.lon, gnss.lat, pred_sse_disp_bounded[0::2], pred_sse_disp_bounded[1::2], color='r', scale=50, label="Modeled")
# Add coastline
ax[1].plot(coast.lon, coast.lat, color="gray", linewidth=0.5)
# Axis limits
ax[1].set(xlim=(lonmin, lonmax), ylim=(latmin, latmax), aspect='equal')
ax[1].legend()
ax[1].yaxis.tick_right()
ax[1].title.set_text("Displacement")
# Add colorbar
cbar_ax = fig.add_axes([0.45, 0.15, 0.025, 0.7])
fig.colorbar(fso, cax=cbar_ax, label="Slip (mm)")
# Add vector scale
ax[1].quiver(235, 37.3, -5, 0, color='k', scale=50)
ax[1].text(233.7, 37.8, '5 mm')
plt.show()