Topographic and tectonic stress inversions with halfspace: Calculating topographic stresses on a fault
Richard Styron
In order to complete the tutorial, please first clone/download the github repo which comes with the DEM, coseismic slip model, and some auxiliary files, all with the correct pathnames.
(Or just download this notebook)
This part is pretty easy. It's just an interpolation of the topographic stress field \(M(x,y,z)\) onto the fault, as represented by a coseismic slip model.
In [1]:
import numpy as np
import pandas as pd
import halfspace.load as hs
import halfspace.projections as hsp
import halfspace.sandbox as hbx
import scipy.ndimage as nd
import h5py
import json
In [2]:
data_dir = '../test_data/'
dem_dir = data_dir + 'dem/'
stress_dir = '../stress_arrays/'
fault_slip_file = data_dir + 'baloch_fault_model.csv'
stress_file = stress_dir + 'baloch_topo_stress.h5'
dem_metadata_file = dem_dir + 'baloch_dem_metadata.json'
Load files¶
In [3]:
d_meta = json.load(open(dem_metadata_file, 'r')) # DEM metadata dictionary
fault_df = pd.read_csv(fault_slip_file, index_col=0) # fault slip DataFrame
stress_db = h5py.File(stress_file, 'r') # open wrapper for stress files
stress_shape = stress_db['xx_MPa'][:,:,0].shape # size of the stress arrays
In [4]:
x_res = int(d_meta['x_res_m']) # get info from metadata dict
y_res = int(d_meta['y_res_m'])
z_res = 1000
z_min = int(d_meta['x_res_m'])
z_max = z_min + 26000 # make sure this matches topo stress calcs
z_len = (z_max - z_min) / z_res + 1
z_vec = np.linspace(z_min, z_max, num=z_len)
topo_res = x_res
Calculate some more important numbers¶
In [5]:
d_meta
Out[5]:
{'north_max': 3329810.0066404096,
'upper_left_y': 3329810.0066404096,
'y_res_m': 445.9049868175547,
'lat_min': 23.993050290002728,
'lon_max': 68.27979759187075,
'east_max': 1009621.1791108495,
'x_rotation': 0.0,
'lat_max': 29.993372267246755,
'x_res_m': 445.9049868175545,
'lon_min': 62.000053658151806,
'east_min': 398285.44218398223,
'y_rotation': 0.0,
'n_cols': 1371,
'north_min': 2653818.0466249967,
'upper_left_x': 398285.44218398223,
'n_rows': 1516}
In [6]:
e_range_dem = np.linspace(d_meta['east_min'], d_meta['east_max'], # orig dem
num=d_meta['n_cols'])
n_range_dem = np.linspace(d_meta['north_min'], d_meta['north_max'],
num=d_meta['n_rows'])
e_range = hs._centered(e_range_dem, stress_shape[1]) # cut down to final size
n_range = hs._centered(n_range_dem, stress_shape[0])
x0 = e_range[0] # first index
y0 = n_range[0]
Interpolate stress tensor components onto fault¶
First we have to calculate the coordinates of the fault patches in 'pixel coordinates' of the DEM for interpolation
In [7]:
fault_xyz = hbx.coord_map_inverse_3d([ fault_df.east_utm.values,
fault_df.north_utm.values,
fault_df.z.values],
x_step=x_res, x_shift=x0,
y_step=y_res, y_shift=y0,
z_step=1, z_shift=fault_df.z.min())
In [8]:
# quick check to see if the conversion is set up right.
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
plt.figure(figsize=(14,12))
plt.pcolormesh(stress_db['zz_MPa'][:,:,0], cmap=plt.cm.gist_earth)
plt.scatter(fault_xyz[0], fault_xyz[1], c=fault_xyz[2],
s=5, lw=0)
plt.colorbar(label='depth (km)')
plt.axis('equal')
plt.show()
Good, now we can do the interpolation.
In [9]:
comp_list = ['zz', 'xy', 'xz', 'yz', 'xx', 'yy']
for comp in comp_list:
print('interpolating {} stresses'.format(comp) )
fault_df['{}_stress'.format(comp)] = nd.map_coordinates(
stress_db['{}_MPa'.format(comp)][:,:,:],
[fault_xyz[1], fault_xyz[0], fault_xyz[2]],
order=1)
stress_db.close()
interpolating zz stresses interpolating xy stresses interpolating xz stresses interpolating yz stresses interpolating xx stresses interpolating yy stresses
Project stress tensor components onto fault¶
In [10]:
fault_df['tau_dd'] = 0.
fault_df['tau_ss'] = 0.
fault_df['sig_nn'] = 0.
stress_tens_xyz = {}
for i, fault_patch in enumerate(fault_df.index):
stress_tens_xyz[i] = hsp.make_xyz_stress_tensor(
sig_xx = fault_df.xx_stress.iloc[i],
sig_yy = fault_df.yy_stress.iloc[i],
sig_zz = fault_df.zz_stress.iloc[i],
sig_xy = fault_df.xy_stress.iloc[i],
sig_xz = fault_df.xz_stress.iloc[i],
sig_yz = fault_df.yz_stress.iloc[i])
fault_df['tau_dd'].iloc[i] = hsp.dip_shear_stress_from_xyz(
fault_df.strike.iloc[i], fault_df.dip.iloc[i],
stress_tens_xyz[i] )
fault_df['tau_ss'].iloc[i] = hsp.strike_shear_stress_from_xyz(
fault_df.strike.iloc[i], fault_df.dip.iloc[i],
stress_tens_xyz[i] )
fault_df['sig_nn'].iloc[i] = hsp.normal_stress_from_xyz(
fault_df.strike.iloc[i], fault_df.dip.iloc[i],
stress_tens_xyz[i] )
In [11]:
fault_df.head()
Out[11]:
| lat | lon | rake | rise | slip_m | rupture_time | x | y | z | dip | ... | north_utm | zz_stress | xy_stress | xz_stress | yz_stress | xx_stress | yy_stress | tau_dd | tau_ss | sig_nn | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 27.1152 | 65.5036 | -2 | 11.5 | 1.2632 | 11.494 | 17.7058 | 27.2426 | 0.2171 | 60 | ... | 3001667.271426 | 34.409794 | -3.208107 | 4.297898 | -4.019554 | 25.450635 | 28.030340 | -5.473018 | -3.743854 | 34.502948 |
| 1 | 27.0489 | 65.4721 | 7 | 7.0 | 1.2632 | 8.880 | 14.5800 | 19.8786 | 0.2171 | 60 | ... | 2994257.958404 | 30.553797 | -1.898560 | 1.930903 | -1.001144 | 23.029290 | 24.616245 | -3.646218 | -1.720015 | 27.994459 |
| 2 | 26.9826 | 65.4406 | 10 | 4.0 | 4.0000 | 6.387 | 11.4541 | 12.5145 | 0.2171 | 60 | ... | 2986849.691310 | 28.060829 | -1.194816 | 1.097168 | -1.762246 | 20.601227 | 22.731160 | -3.566387 | -1.978959 | 24.825576 |
| 3 | 26.9164 | 65.4090 | 10 | 4.0 | 5.0000 | 4.239 | 8.3283 | 5.1505 | 0.2171 | 60 | ... | 2979453.359426 | 26.534736 | -2.712429 | 1.860519 | -1.841293 | 20.933961 | 22.053082 | -2.722381 | -2.464341 | 26.031898 |
| 4 | 26.8501 | 65.3775 | 13 | 11.5 | 6.0000 | 3.218 | 5.2024 | -2.2136 | 0.2171 | 60 | ... | 2972047.177332 | 27.472174 | -1.913355 | 3.559284 | -1.770538 | 20.819765 | 20.779945 | -4.271303 | -1.258188 | 26.947084 |
5 rows × 23 columns
Done with calculations. Now save and we're through.¶
In [12]:
fault_df.to_csv(fault_slip_file)
- Installation and configuration
- Preparing a Digital Elevation Model file
- Preparing the coseismic slip model
- Calculating the topographic stress field
- Calculating the topographic stresses on the fault
- Bayesian inversion of tectonic stress pt. 1: Rake inversion
- Bayesian inversion of tectonic stress pt. 2: Conditions at failure