Rocks and Water

• About
• Bio
• Research
• Publications
• Global Block Model
• Archive

Low-angle Normal Fault Likelihood IPython Notebook

Richard Styron

Sat 18 January 2014

This is an IPython Notebook version of a paper I have just submitted to Geophysical Research Letters. This is the complete paper, but with all the code embedded, and is executable. However, it does need a Python module called eq_stats.py as well as some of the typical Python scientific stack. The notebook can be downloaded from github here. I will also probably post the updated version once it's through review.

Estimated likelihood of observing a large earthquake on a continental low-angle normal fault, and implications for low-angle normal fault activity

Richard Styron* and Eric Hetland

Dept. Earth and Environmental Sciences, University of Michigan

*[email protected]

Abstract

The lack of observed continental earthquakes that clearly occurred on low-angle normal faults (LANFs) may indicate that these structures are not seismically active, or that these earthquakes are simply rare events. To address this, we compile all potentially active continental LANFs (twenty in total) and calculate the likelihood of observing a significant earthquake on them over periods of 1-100 years. This probability depends on several factors including the frequency-magnitude distribution. For either a characteristic or Gutenberg-Richter distribution, we calculate a probability of about 0.5 that an earthquake greater than $M6.5$ (large enough to avoid ambiguity in dip angle) will be observed on any LANF in a period of 35 years, which is the current length of the Global CMT catalog. We then use Bayes' Theorem to illustrate how the absence of observed significant LANF seismicity over the catalog period moderately decreases the likelihood that the structures generate large earthquakes.

1 Introduction

Low-angle normal faults (LANFs), with dips less than 30$^\circ$, are well described in the geologic record. They are thought to play an important role in accommodating large-magnitude continental extension [Howard and John, 1987] and crustal thinning [Lister et al., 1986], and their recognition has been a major development in continental tectonics [Wernicke, 2009]. However, despite widespread field observations of inactive LANFs and their central role in extensional tectonic theory, they remain enigmatic and contentious structures, and it is not clear if they are seismically active at low dip angles in the upper crust. This is for two reasons: because brittle faulting on LANFs is in apparent conflict with standard Andersonian rock mechanical theory as typically applied to the upper crust [Axen, 2004], and because observations of active faulting on LANFs are sparse and at times ambiguous [Wernicke, 1995]. A considerable amount of research has been performed to address the former concern, reconciling LANF slip with rock mechanics [e.g., Axen and Bartley, 1997; Collettini, 2011]. The latter issue is highlighted by studies that have searched the focal mechanism catalogs and found no normal faulting earthquakes with focal mechanisms and surface ruptures clearly indicating slip on planes $\le30^\circ$ [Jackson, 1987; Collettini and Sibson, 2001], which is taken as conclusive evidence that LANFs are inactive or aseismic. However, the lack of observed seismic slip on continental LANFs may be simply because they are rare structures with long recurrence intervals, so earthquakes on them are very infrequent. Without knowing the likelihood of observing an LANF rupture in a time window of a few decades, it is not clear if an empty search result is strong evidence against LANF seismicity. If this likelihood is known, though, Bayesian probability theory provides a framework for quantifying how the negative search results impact the probability that LANFs are seismogenic.

In this work, we estimate the maximum likelihood of a significant LANF event occurring in time windows from 1 to 100 years, and then we interpret the lack of observed LANF seismicity in a quantified, probabilistic context using Bayes' Theorem. We estimate the maximum observation likelihood by treating all potentially active LANFs described in the literature as seismically active at their surface dip angles throughout the upper crust. Under these assumptions, we create synthetic earthquake catalogs with both Gutenberg-Richter and `characteristic' frequency--magnitude distributions, using each fault's geometry and slip rate. We then calculate the probability of observing earthquakes on at least one LANF over different observation periods. Finally, we use Bayes' Theorem to incorporate the negative catalog search results and the observance likelihood to show how the negative results reduce the probability that LANFs are seismically active, but do not bring the final probability to zero.

1.1 LANF Slip, Mohr-Coulomb Failure Theory, and Earthquakes

Areas of the crust undergoing active extension are generally assumed to have a subvertical maximum compressive stress. Mohr-Coulomb theory, as applied to the crust, predicts that a fault with a typical coefficient of friction for rocks (0.6-0.8) should lock up if it is oriented at an angle greater than 60$^\circ$ to the maximum compressive stress (i.e., fault dips less than 30$^\circ$), and new, optimally oriented faults should form [Sibson, 1985]. Therefore, for normal faults with dips less than 30$^\circ$, either much lower fault friction or elevated pore fluid pressure is required for fault slip.

Evidence for seismic slip on LANFs is sparse. This is partly due to the ambiguity of the rupture plane in earthquake focal mechanisms, as a focal mechanism with a low angle nodal plane will also by definition have a high angle nodal plane. Without ancillary information indicating which nodal plane corresponds to the slip surface, searches of earthquake catalogs cannot yield unique results as to whether they contain LANF events. Several collections of normal fault earthquakes with known surface breaks Jackson, 1987; Collettini and Sibson, 2001], thereby resolving dip ambiguity, contain no low-angle events, although we note the total number of events in these collections are small ($\le$ 25 events). Some candidate LANF events exist, but they are undersea [e.g., Abers, 2001] or difficult to verify [e.g., Doser, 1987].

2. Potentially Active LANFs

Over the past decade or so, many field studies have found evidence for LANF activity in orogens throughout the world. These studies typically find arrays of Quaternary normal fault scarps on the fault traces and/or in the hanging walls of mapped or inferred low-angle detachment faults [e.g., Axen et al., 1999] . Some studies also have bedrock thermochronology data from the exhumed detachment footwalls that are suggestive of ongoing rapid exhumation [e.g., Sundell et al. 2013], although this data does not preclude a recent cessation of faulting. In some cases, additional evidence for LANF activity comes from geophysical data such as GPS geodesy [e.g., Hreinsdóttir and Bennett, 2009] and seismic waves [e.g., Doser, 1987].

from IPython.core.display import Image 
Image(filename='active_lanfs_map_insets.png')

Figure 1: Map of known, potentially active continental LANFs (blue lines), with insets showing the physiographic context of the faults. DXV=Dixie Valley fault. PV=Panamint Valley fault. DV=Death Valley fault. CD=Cañada David detachment. SD=Sevier Desert detachment. CB=Cordillera Blanca detachment. AT=Alto-Tiberina fault. KZ=Kuzey detachment. GN=Guney detachment. KS=Kongur Shan fault. LP=Leo Pargil detachment. GM=Gurla Mandhata detachment. NLR=North Lunggar detachment. SLR=South Lunggar detachment. PXN=Pum Qu--Xainza north fault. PXQ=Pum Qu--Xainza Qingdu fault. NQTL=Nyainqentanglha detachment. PP=Pompangeo detachment. TK=Tokorondo detachment. DD=Dayman Dome.

We have compiled all potentially active LANFs with known subareal fault traces from a thorough review of the literature, finding twenty total (Figure 1). We have then mapped the approximate fault traces into a GIS file (available at https://github.com/cossatot/LANF_gis), with metadata such as slip rate and source. Though the fault traces of many LANFs considered here are obscured by vegetation, others display large fault scarps in Quaternary sediments, particularly those in Tibet [e.g., Styron et al., 2013; Kapp et al., 2005] and the western US [e.g., Axen et al., 1999, Hayman et al., 2003], which are commonly interpreted as evidence for past seismic slip. About half are in Tibet, consistent with hypotheses that LANFs and metamorphic core complexes form in areas of hot, thick crust [e.g., Buck, 1991]. The rest are distributed through other areas of active continental extension: the North American Basin and Range, the Malay Archipelago, western Turkey, Italy, and Peru.

Several of the most-commonly cited candidates for seismically active LANFs were not included because they do not have a clearly-defined, mappable fault trace, which is necessary for our earthquake likelihood calculations. These include the 1995 Aigion, Greece earthquake fault [Bernard et al., 1997] and other potential LANFs underneath the Gulf of Corinth, and the 1952 Ancash, Peru earthquake fault [Doser, 1987]. Furthermore, though submarine core complexes with superficially low-angle detachments are well-described in the literature and some of these structures may have produced recent earthquakes [Abers, 2001], we do not include these in our calculations for several reasons: because mid-ocean ridges have not been structurally mapped with the completeness or resolution of subareal extensional provinces, it is not currently possible to come up with a reasonably complete inventory of ocean LANFs; without high-resolution structural mapping and geodesy of oceanic LANFs, it is not possible to determine which structures in a mid-ocean ridge segment are currently active (seismically or not), and it is difficult to confidently associate particular earthquakes with a specific fault, given the high spatial density of normal faults at mid-ocean ridges.

To estimate the likelihood of observing a significant earthquake on an individual LANF over some contiguous time window of length $t$ (in years), we perform a Monte Carlo simulation in which we create 4000 synthetic time series of earthquakes, with unique values for fault geometry and slip rate for each time series. Then, for each time series we calculate the fraction of unique time windows of length $t$ in which an earthquake as large or larger than a given magnitude occurs. We take this value as the probability of observing an earthquake greater than or equal to moment magnitude \emph{M} over time period $t$, which we will refer to in general as $P(M,t)$. All calculations are performed with Python, using the Numpy [Oliphant, 2007], IPython [Pérez and Granger, 2007], Pandas [McKinney, 2010], and Joblib Parallel [Varoquaux and Grisel, 2009] packages. All code and data for this project is available at https://github.com/cossatot/lanf%5C_earthquake%5C_likelihood/.

The geometry for each fault is estimated based on the length of the fault trace, the dip of the fault, and the estimated fault locking depth in the area. The fault is treated as planar for simplicity of calculations, even though the exposed footwalls of many detachment faults are nonplanar. We determine the fault length by measuring the approximate length of the mapped fault trace perpendicular to the assumed extension direction; for faults that change dip significantly along strike, we only consider the low-angle segments of the fault. Values for the dip are taken from the literature in most cases, and measurements of the dip of footwall triangular facets (interpreted as the exhumed fault plane) from SRTM data otherwise. In all cases, ranges of fault geometries are considered, encompassing the degree to which the values are known. The fault locking depth is assumed to be 10 km in the absence of other evidence (such as a geodetic study, [e.g., Hreinsdóttir and Bennett, 2009]).

Slip rates of the 20 LANFs are gathered from the literature if possible, or given broad ranges if not (e.g., 1--10 mm yr$^{-1}$). In the Monte Carlo simulation, samples for slip rate and dip are drawn from uniform distributions defined by the maximum and minimum values. Based on field observations, some faults have dip ranges that go above 30$^\circ$, although for these faults dip values are sampled from the minimum to 30$^\circ$, as here we only consider slip on faults shallower than 30$^\circ$. The resulting probabilities on these faults are then multiplied by the fraction of the dip range that is $\le30^\circ$.

Each synthetic earthquake sequence is generated by randomly sampling either 50,000 events from a tapered Gutenberg-Richter (GR) distribution with corner magnitude $M_c = 7.64$ and $\beta = 0.65$ (from values estimated by Bird and Kagan [2004] for continental rifts), or a 25,000 events from `characteristic' distribution. It is not certain which distribution more appropriately describes seismicity on a single LANF, though studies of many individual fault rupture histories suggests that the characteristic distribution is more accurate [e.g., Hecker et al., 2003]. The smaller number of samples drawn from the characteristic distribution is due to the increased computation time associated with a higher proportion of large events, leading to much longer time series for a given number of events. The samples are taken from an interval $M = [5.0, \, M_{max}]$, where $M_{max}$ is the moment magnitude associated with 15 m of slip over the given fault plane. We use the standard relations between fault slip, $D$, and moment magnitude, $M$, given by

\begin{equation} M_o = \mu L z D \,/ \, \sin \delta \end{equation}

and

\begin{equation} M = 2/3 \; \log_{10} (M_o) - 6 \end{equation}

where $L$ is the fault length, $z$ is the seismogenic thickness, $\delta$ is the fault dip, $\mu = 30$ GPa is the shear modulus, and $M_o$ is the seismic moment in N m [Aki and Richards, 2002; Kagan, 2003]. The characteristic distribution has a large-magnitude mode corresponding to $D$ = 1.5 m on the fault, a typical slip distance for normal fault events [e.g. Wesnousky, 2008]. The distributions are shown in Figure 2.

from IPython.core.display import Image
Image(filename='F-Ms.jpg')

Figure 2: Gutenberg-Richter and characteristic frequency-magnitude distributions for the South Lunggar detachment

These calculations rely on two important assumptions that warrant some discussion. The first is that each earthquake ruptures the entire fault patch uniformly. Though this is unlikely fault behavior, the long-term statistical distribution of earthquake recurrence is insensitive to assumptions about slip distribution in individual events as long as earthquakes are unclustered in time (the second assumption discussed below). Specifically, if $n$ different, equal fault patches rupture independently, each requires $n$ times the interseismic strain accumulation time to rupture with an earthquake of magnitude $M$ compared to the accumulation time for a single fault rupturing uniformly with much lower coseismic slip in each earthquake. Thus, magnitude $M$ events would happen with the same long-term frequency. The next assumption is that earthquakes are ordered randomly and separated by the time necessary for sufficient strain to accumulate for each earthquake to occur. This means that foreshock and aftershock sequences and other types of event clustering are not taken into account. However, the modal inter-event times for earthquakes $\ge M \,6$ or so are greater than a hundred years for most LANFs, so the ordering of events does not impact the results, as this is longer than our maximum observation window. Furthermore, any clustering resulting in event spacing less than the observation window would decrease $P(M,t)$, and here we choose to calculate the maximum $P(M,t)$ using the simplest assumptions, rather that choose the model assumptions such that the calculated probabilities are the minimum.

Now for some coding...

Setting up the problem

Import necessary modules

%pylab inline
%config InlineBackend.figure_format = 'retina'

Populating the interactive namespace from numpy and matplotlib

import sys
sys.path.append('../eq_stats')

import numpy as np
import pandas as pd
import eq_stats as eqs  # Custom module for this project
import time
from joblib import Parallel, delayed
from itertools import chain

Read in fault data table

Makes a Pandas dataframe of fault data (length, slip rates, etc.)

f = pd.read_csv('../data/lanf_stats.csv', index_col=0)

f

	region	L_km	z_km	slip_rate_mm_a	sr_err_mm_a	dip_deg	dip_err_deg
fault
kongur_shan	tibet	147	10.0	6.5	2.0	30.0	10.0
leo_pargil	tibet	56	10.0	1.0	0.8	25.0	10.0
gurla_mandhata	tibet	44	10.0	4.5	2.0	33.0	10.0
s_lunggar	tibet	20	10.0	2.5	0.5	22.0	5.0
n_lunggar	tibet	32	10.0	2.0	1.0	25.0	5.0
pqx_qingdu	tibet	15	10.0	1.0	0.5	30.0	5.0
pqx_n	tibet	11	10.0	1.0	0.5	25.0	7.5
nqtl	tibet	115	10.0	2.5	1.5	30.0	8.0
alto_tiberina	italy	60	4.0	2.4	0.3	25.0	5.0
death_valley	w_us	38	10.0	2.0	1.0	27.5	8.5
panamint_valley	w_us	58	10.0	3.5	1.5	20.0	10.0
canada_david	w_us	20	10.0	3.0	2.0	20.0	10.0
sevier_desert	w_us	98	8.5	2.5	1.0	13.5	1.5
dixie_valley	w_us	21	10.0	0.6	0.4	20.0	10.0
dayman_dome	se_asia	80	10.0	5.0	4.8	20.0	10.0
papangeo	se_asia	27	10.0	5.0	4.8	20.0	10.0
tokorondo	se_asia	72	10.0	5.0	4.8	10.0	8.0
cordillera_blanca	andes	100	10.0	6.0	4.0	30.0	10.0
kuzey	turkey	130	10.0	1.0	0.9	17.0	2.0
guney	turkey	124	10.0	1.0	0.9	7.0	7.0

20 rows × 7 columns

Define some variables to be used later

n_cores = -3
n_eq_samp = 3e4 # number of earthquakes in time series
time_window = np.hstack( (1, np.arange(5, 105, step=5) ) ) # observation times
mc_iters = 2e3 # number of Monte Carlo iterations [50k in published calcs]
mc_index = np.arange(mc_iters, dtype='int')
mc_cols = ['dip', 'Ddot'] + [t for t in time_window]
max_eq_slip = 15 #m
char_eq_slip= 1.5 #m
Mc = 7.64 # Hypothetical corner magnitude for Continental Rift Boundaries

Make list of minimum search magnitude $M_{min}$, and then make MultiIndex for Pandas dataframes

min_M_list = [5, 5.5, 6, 6.5, 7, 7.5]

df_ind_tuples = [[i, M] for i in mc_index for M in min_M_list]
df_multi_ind = pd.MultiIndex.from_tuples(df_ind_tuples, names=['mc_iter','M'])

rec_int_bins = np.logspace(1, 5)

Calculate $P(M,t)$ for Gutenberg-Richter frequency-magnitude distributions

Define a function for Joblib Parallel to calculate probabilities for each iteration.

Function is defined here so it can access all variables generated by script, not just passed variables. This makes the code cleaner even if it's not very abstracted.

Here is what this function does:

• Get the dip, Ddot and maximum earthquake magnitude for each iteration.
• Take this info and make the earthquake sequence:
  • Take the max earthquake magnitude and make a frequency-magnitude distribution based on a Gutenburg-Richter exponential model.
  • Take 50k samples from this distribution,
• Make an earthquake time series form the EQ sequence
  • Calculate the interseismic strain accumulation time for each event
  • Separate each earthquake in the sequence with the appropriate number of years with no events.
• Calculate the probability of observation
  • Run a rolling maximum for each $t$ in [1, 5, 10, 15, ..., 95, 100]
  • Calculate the observation probability above $M_{min}$ in [5, 5.5, 6, 6.5, 7, 7.5]
• Calculate inter-event times for EQs $\ge \, M$

def calc_iter_probs(iter):
    df_iter = fdf.loc[iter].copy()
    df_iter['dip'] = mc_d['dip_samp'][iter]
    df_iter['Ddot'] = mc_d['Ddot_samp'][iter]

    # Generate EQ sample/sequence from F(M) dist.
    m_vec = np.linspace(5, mc_d['max_M'][iter], num=1000)
    fm_vec = eqs.F(m_vec, Mc=Mc)
    M_samp = eqs.sample_from_pdf(m_vec, fm_vec, n_eq_samp)
    Mo_samp = eqs.calc_Mo_from_M(M_samp)
    
    # Make time series of earthquakes, including no eq years
    recur_int = eqs.calc_recurrence_interval(Mo=Mo_samp, 
                                             dip=mc_d['dip_samp'][iter],
                                             slip_rate=mc_d['Ddot_samp'][iter],
                                             L=params['L_km'],
                                             z=params['z_km'])

    cum_yrs = eqs.calc_cumulative_yrs(recur_int)
    eq_series = eqs.make_eq_time_series(M_samp, cum_yrs)
    
    # calculate probability of observing EQ in time_window
    for t in time_window:
        roll_max = pd.rolling_max(eq_series, t)
        df_iter[t] = (eqs.get_probability_above_value(roll_max, min_M_list)
                      * mc_d['dip_frac'])
        
    # calculate histgrams of recurrence intervals
    rec_int_counts_df = rec_int_df.loc[iter].copy()
    for mm in np.array(min_M_list):
        ints = np.diff( np.where(eq_series >= mm) )
        rec_int_counts_df.loc[mm] = np.histogram(ints, bins=rec_int_bins)[0]

Iterate through the faults in the fault database, doing all the calculations for each.

The setup of this for loop is basically this:

• Make DataFrame for each fault.

  • Columns are dip, Ddot, and observation time windows.
  • Rows are values for each Monte Carlo iteration. Values for time windows are calculated probabilities.

• Calculate maximum earthquake magnitude for each MC iteration.

• Run the above 'calc_iter_probs' function (parallelized over the MC iterations) and concatenate the results

Note! Running the calcs here is a looong process!

for fault in list(f.index):
    fdf = pd.DataFrame(index=df_multi_ind, columns=mc_cols, dtype='float')
    params = f.loc[fault]
    mc_d = {}
    mc_d['dip_samp'], mc_d['dip_frac'] = eqs.dip_rand_samp( params['dip_deg'], 
                                                         params['dip_err_deg'], 
                                                         mc_iters)

    mc_d['Ddot_samp'] = eqs.Ddot_rand_samp(params['slip_rate_mm_a'],
                                           params['sr_err_mm_a'], mc_iters)

    mc_d['max_Mo'] = eqs.calc_Mo_from_fault_params(L=params['L_km'], 
                                                   z=params['z_km'], 
                                                   dip=mc_d['dip_samp'], 
                                                   D=max_eq_slip)

    mc_d['max_M'] = eqs.calc_M_from_Mo(mc_d['max_Mo'])
    
    rec_int_df = pd.DataFrame(columns = rec_int_bins[1:],
                              index=df_multi_ind, dtype='float')
    t0 = time.time()
    out_list = Parallel(n_jobs=-3)( delayed( calc_iter_probs)(ii) 
                                    for ii in mc_index)
    print 'done with', fault, 'parallel calcs in {} s'.format((time.time()-t0))
    for ii in mc_index:
        fdf.loc[ii][:] = out_list[ii][0]
        rec_int_df.loc[ii][:] = out_list[ii][1]
        
    fdf.to_csv('../results/{}_test.csv'.format(fault))

Calculate $P(M,t)$ for faults with characteristic frequency-magnitude distributions

def calc_iter_probs(ii):
    df_iter = fdf.loc[ii].copy()
    df_iter['dip'] = mc_d['dip_samp'][ii]
    df_iter['Ddot'] = mc_d['Ddot_samp'][ii]

    # Generate EQ sample/sequence from F(M) dist.
    m_vec = np.linspace(5, mc_d['max_M'][ii], num=1000)
    fm_vec = eqs.F_char(m_vec, Mc=Mc, char_M=mc_d['char_M'][ii])
    M_samp = eqs.sample_from_pdf(m_vec, fm_vec, n_eq_samp)
    Mo_samp = eqs.calc_Mo_from_M(M_samp)
    
    # Make time series of earthquakes, including no eq years
    recur_int = eqs.calc_recurrence_interval(Mo=Mo_samp, 
                                             dip=mc_d['dip_samp'][ii],
                                             slip_rate=mc_d['Ddot_samp'][ii],
                                             L=params['L_km'],
                                             z=params['z_km'])

    cum_yrs = eqs.calc_cumulative_yrs(recur_int)
    eq_series = eqs.make_eq_time_series(M_samp, cum_yrs)
    
    # calculate probability of observing EQ in time_window
    for t in time_window:
        roll_max = pd.rolling_max(eq_series, t)
        df_iter[t] = (eqs.get_probability_above_value(roll_max, min_M_list)
                      * mc_d['dip_frac'])
        
    # calculate histgrams of recurrence intervals
    rec_int_counts_df = rec_int_df.loc[ii].copy()
    for mm in np.array(min_M_list):
        ints = np.diff( np.where(eq_series >= mm) )

for fault in list(f.index):
    fdf = pd.DataFrame(index=df_multi_ind, columns=mc_cols, dtype='float')
    params = f.loc[fault]
    mc_d = {}
    mc_d['dip_samp'], mc_d['dip_frac'] = eqs.dip_rand_samp( params['dip_deg'], 
                                                         params['dip_err_deg'], 
                                                         mc_iters)

    mc_d['Ddot_samp'] = eqs.Ddot_rand_samp(params['slip_rate_mm_a'],
                                           params['sr_err_mm_a'], mc_iters)

    mc_d['max_Mo'] = eqs.calc_Mo_from_fault_params(L=params['L_km'], 
                                                   z=params['z_km'], 
                                                   dip=mc_d['dip_samp'], 
                                                   D=max_eq_slip)

    mc_d['max_M'] = eqs.calc_M_from_Mo(mc_d['max_Mo'])
    
    mc_d['char_Mo'] = eqs.calc_Mo_from_fault_params(L=params['L_km'], 
                                                    z=params['z_km'], 
                                                    dip=mc_d['dip_samp'], 
                                                    D=char_eq_slip)
    
    mc_d['char_M'] = eqs.calc_M_from_Mo(mc_d['char_Mo'])
    
    rec_int_df = pd.DataFrame(columns = rec_int_bins[1:],
                              index=df_multi_ind, dtype='float')
    t0 = time.time()
    out_list = Parallel(n_jobs=n_cores)( delayed( calc_iter_probs)(ii) 
                                    for ii in mc_index)
    print 'done with', fault, 'parallel calcs in {} s'.format((time.time()-t0))
    for ii in mc_index:
        fdf.loc[ii][:] = out_list[ii][0]
        rec_int_df.loc[ii][:] = out_list[ii][1]
        
    fdf.to_csv('../results/{}_char.csv'.format(fault))

Examining individual fault results

Load datasets into Pandas dataframes

test with one:

pv = pd.read_csv('../results/panamint_valley_gr.csv', index_col=[0])
pv_c = pd.read_csv('../results/panamint_valley_char.csv', index_col=[0])

tw_xarray = np.tile(time_window, (mc_iters,1))

tw_cols = list(time_window.astype('str'))

make some plotting functions

def plot_probs(df, lw=0.5, alpha=0.2):

    df_5 = df[df.M==5]
    df_55 = df[df.M==5.5]
    df_6 = df[df.M==6]
    df_65 = df[df.M==6.5]
    df_7 = df[df.M==7]
    df_75 = df[df.M==7.5]

    aa = plot(tw_xarray.T, df_5[tw_cols].T, 'b', 
              lw=lw, alpha=alpha)
    bb = plot(tw_xarray.T, df_55[tw_cols].T, 'c', 
              lw=lw, alpha=alpha)
    cc = plot(tw_xarray.T, df_6[tw_cols].T, 'g', 
              lw=lw, alpha=alpha)
    dd = plot(tw_xarray.T, df_65[tw_cols].T, 'gold', 
              lw=lw, alpha=alpha)
    ee = plot(tw_xarray.T, df_7[tw_cols].T, 'orange', 
              lw=lw, alpha=alpha)
    ff = plot(tw_xarray.T, df_75[tw_cols].T, 'r', 
              lw=lw, alpha=alpha)

    return aa, bb, cc, dd, ee, ff

def pmf_x_secs(df, tw, n_bins=101, hist_type='stepfilled', alpha=0.5):
    hist_bins = np.linspace(0, 1, num=n_bins)
    tw = str(tw)
    
    aa = df[df.M==5][tw].hist(bins=hist_bins, histtype=hist_type, 
                              color='b', alpha=alpha)
    bb = df[df.M==5.5][tw].hist(bins=hist_bins, histtype=hist_type, 
                                color='c', alpha=alpha)
    cc = df[df.M==6][tw].hist(bins=hist_bins, histtype=hist_type, 
                              color='g', alpha=alpha)
    dd = df[df.M==6.5][tw].hist(bins=hist_bins, histtype=hist_type, 
                                color='gold', alpha=alpha)
    ee = df[df.M==7][tw].hist(bins=hist_bins, histtype=hist_type, 
                              color='orange', alpha=alpha)
    ff = df[df.M==7.5][tw].hist(bins=hist_bins, histtype=hist_type, 
                                color='red', alpha=alpha)
    
    return aa, bb, cc, dd, ee, ff

tw=35
nbins=1001
x_sec_lims = [0,0.4,0,400]

figure()

subplot2grid((3,2), (0,0), colspan=2)
axvline(tw, lw=0.5, color='grey')
plot_probs(pv, lw=0.2, alpha=0.2)

subplot2grid((3,2), (1,0), colspan=2)
axvline(tw, lw=0.5, color='grey')
plot_probs(pv_c, lw=0.2, alpha=0.2)

subplot2grid((3,2), (2,0))
pmf_x_secs(pv, tw, n_bins=nbins)
axis(x_sec_lims)

subplot2grid((3,2), (2,1))
pmf_x_secs(pv_c, tw, n_bins=nbins)
axis([0,0.1,0,400])

show()

Figure 3. a: Probabilities of observing an earthquake greater than or equal to a given moment magnitude $M$ over a given observation window on the Panamint Valley fault, for the Gutenberg-Richter distribution. b: Probabilities of observing an earthquake greater than or equal to a given moment magnitude $M$ over a given observation window on the Panamint Valley fault, for the characteristic distribution. Note the change in the scale of the y axis. c: Cross-section through a at $t$ = 35 years, showing the distributions of $P(M)$. d: Cross-section through b at $t$ = 35 years, showing the distributions of $P(M)$.

...Back to the article

The results for faults with a GR frequency-magnitude distribution indicate that it is unlikely that any individual fault would have an earthquake greater than $M \, 5$ in any observation time window up to 100 years. As an example, the results for the Panamint Valley fault are shown in Figure 3a; this fault has the highest $P(M,t)$ of any of the well-studied LANFs. The probability of observing a $\ge M \, 6.0$ event on the Panamint Valley fault is about 0.5 for $t$ = 100 years, and about 0.15 for $t$ = 35 years, which is the current length of the Global CMT catalog. As expected given the GR distribution, $P(M,t)$ is much higher for smaller, more frequent events than for larger events.

The results for faults with a characteristic frequency-magnitude distribution yield much lower $P(M,t)$ for small to moderate events, but $P(M,t)$ is higher for large events(Figure 3b,d); this is because the earthquake sequences are dominated by large, infrequent events, so the inter-event times for moderate events are several times greater. For the Panamint Valley fault, $P(M\ge5,t=35)$ is about 0.07 (versus 0.25 for the GR distribution), but $P(M\ge 7, t=35)$ is around 0.025 (versus essentially zero for the GR distribution). As the characteristic distribution likely better represents earthquakes on an individual large fault, these results suggest that is very unlikely that we would expect to capture any significant seismicity on an single LANF in the focal mechanism catalogs. A similar conclusion was found by Wernicke [1995] based on a simple calculation, assuming perfectly repeating large earthquakes on an idealized fault.

3.2 Earthquake likelihood on all LANFs

To calculate the probability of observing at least one earthquake on any of these LANFs during a given time period, we first assume that seismicity on each fault is independent and uncorrelated with seismicity on all other faults. This assumption is likely true for most faults, but may not be true for the few proximal faults, though it is unclear how these faults may interact such that an appropriate joint probability may be calculated. We determine the probability for each time window and minimum magnitude with the equation

$P_{AT \, or \, LP\, or\, \ldots \, or \, DV} = 1 - (Q_{AT} \cdot Q_{LP} \cdot \ldots \, \cdot Q_{DV})$

where $P_{AT}$ is the probability of observing an earthquake on a single LANF (e.g., the Alto-Tiberina fault), and $Q_{AT} = 1 - P_{AT}$. Equation 3 is the union of probabilities for non-mutually exclusive random events.

More coding...

Load the Gutenberg-Richter results

ks = pd.read_csv('../results/kongur_shan_gr.csv', index_col=0)
gm = pd.read_csv('../results/gurla_mandhata_gr.csv', index_col=0)
slr = pd.read_csv('../results/s_lunggar_gr.csv', index_col=0)
nlr = pd.read_csv('../results/n_lunggar_gr.csv', index_col=0)
pqxn = pd.read_csv('../results/pqx_n_gr.csv', index_col=0)
pqxq = pd.read_csv('../results/pqx_qingdu_gr.csv', index_col=0)
lp = pd.read_csv('../results/leo_pargil_gr.csv', index_col=0)
nqtl = pd.read_csv('../results/nqtl_gr.csv', index_col=0)
at = pd.read_csv('../results/alto_tiberina_gr.csv', index_col=0)
dv = pd.read_csv('../results/death_valley_gr.csv', index_col=0)
cd = pd.read_csv('../results/canada_david_gr.csv', index_col=0)
sd = pd.read_csv('../results/sevier_desert_gr.csv', index_col=0)
dxv = pd.read_csv('../results/dixie_valley_gr.csv', index_col=0)
dd = pd.read_csv('../results/dayman_dome_gr.csv', index_col=0)
pp = pd.read_csv('../results/papangeo_gr.csv', index_col=0)
tk = pd.read_csv('../results/tokorondo_gr.csv', index_col=0)
cb = pd.read_csv('../results/cordillera_blanca_gr.csv', index_col=0)
kz = pd.read_csv('../results/kuzey_gr.csv', index_col=0)
gn = pd.read_csv('../results/guney_gr.csv', index_col=0)

Load the characteristic results

ks_c = pd.read_csv('../results/kongur_shan_char.csv', index_col=0)
gm_c = pd.read_csv('../results/gurla_mandhata_char.csv', index_col=0)
slr_c = pd.read_csv('../results/s_lunggar_char.csv', index_col=0)
nlr_c = pd.read_csv('../results/n_lunggar_char.csv', index_col=0)
pqxn_c = pd.read_csv('../results/pqx_n_char.csv', index_col=0)
pqxq_c = pd.read_csv('../results/pqx_qingdu_char.csv', index_col=0)
lp_c = pd.read_csv('../results/leo_pargil_char.csv', index_col=0)
nqtl_c = pd.read_csv('../results/nqtl_char.csv', index_col=0)
at_c = pd.read_csv('../results/alto_tiberina_char.csv', index_col=0)
dv_c = pd.read_csv('../results/death_valley_char.csv', index_col=0)
cd_c = pd.read_csv('../results/canada_david_char.csv', index_col=0)
sd_c = pd.read_csv('../results/sevier_desert_char.csv', index_col=0)
dxv_c = pd.read_csv('../results/dixie_valley_char.csv', index_col=0)
dd_c = pd.read_csv('../results/dayman_dome_char.csv', index_col=0)
pp_c = pd.read_csv('../results/papangeo_char.csv', index_col=0)
tk_c = pd.read_csv('../results/tokorondo_char.csv', index_col=0)
cb_c = pd.read_csv('../results/cordillera_blanca_char.csv', index_col=0)
kz_c = pd.read_csv('../results/kuzey_char.csv', index_col=0)
gn_c = pd.read_csv('../results/guney_char.csv', index_col=0)

make list of faults

f_list = ['ks', 'lp', 'gm', 'slr', 'nlr', 'pqxn', 'pqxq', 'nqtl', 'at',
          'dv', 'pv', 'cd', 'sd', 'dxv', 'dd', 'pp', 'tk', 'cb', 'kz', 'gn']

fault_list = [ks, lp, gm, slr, nlr, pqxn, pqxq, nqtl, at,
              dv, pv, cd, sd,  dxv, dd, pp, tk, cb, kz, gn]

fault_c_list = [ks_c, lp_c, gm_c, slr_c, nlr_c, pqxn_c, pqxq_c, nqtl_c, at_c,
              dv_c, pv_c, cd_c, sd_c, dxv_c, dd_c, pp_c, tk_c, cb_c, kz_c, gn_c]

fc_list = ['ks_c', 'lp_c', 'gm_c', 'slr_c', 'nlr_c', 'pqxn_c', 
                'pqxq_c', 'nqtl_c', 'at_c', 'dv_c', 'pv_c', 'cd_c', 
                'sd_c', 'dxv_c', 'dd_c', 'pp_c', 'tk_c', 'cb_c', 'kz_c', 'gn_c']

Calculate $Q$, probability of not observing an earthquake, for each fault

q = {}

for i, f_ in enumerate(f_list):
    q[f_] = fault_list[i].copy()
    q[f_][tw_cols] = 1 - q[f_][tw_cols]

qc = {}

for i, f_ in enumerate(fc_list):
    qc[f_] = fault_c_list[i].copy()
    qc[f_][tw_cols] = 1 - qc[f_][tw_cols]

Now estimate mutual probabilities

Make list of columns to retain in final dataframe

prob_cols = list( chain( *['M', list(tw_cols)] ) )

calculate $P_{mutual}$ as $1-(Q \cdot Q \cdot Q...)$ and make final dataframe

all_prob_gr = 1 - np.product([qq[tw_cols] for qq in  q.values()])

all_prob_gr['M'] = ks['M']

all_prob_gr = all_prob_gr.reindex_axis(prob_cols, axis=1)

all_c_prob = 1- np.product([qq[tw_cols] for qq in qc.values()])

all_c_prob['M'] = ks['M']

all_c_prob = all_c_prob.reindex_axis(prob_cols, axis=1)

tw=35
nbins=501
x_sec_lims = [0,0.4,0,400]

figure()

subplot2grid((2,2), (0,0))
axvline(tw, lw=0.5, color='grey')
plot_probs(all_prob_gr, lw=0.1, alpha=0.1)

subplot2grid((2,2), (0,1))
axvline(tw, lw=0.5, color='grey')
plot_probs(all_c_prob, lw=0.1, alpha=0.1)

subplot2grid((2,2), (1,0))
pmf_x_secs(all_prob_gr, tw, n_bins=nbins)
#axis(x_sec_lims)

subplot2grid((2,2), (1,1))
pmf_x_secs(all_c_prob, tw, n_bins=nbins)
#axis([0,0.1,0,400])
xlim([0,0.8])
show()

Figure 4. a: Probabilities of observing an earthquake greater than or equal to a given moment magnitude $M$ over a given observation window on any LANF, given a Gutenberg-Richter distribution. b: Probabilities of observing an earthquake greater than or equal to a given moment magnitude M over a given observation window on any LANF, given a characteristic distribution. c: Cross-section through a at $t$ = 35 years showing probability distributions. d: Cross-section through b at $t$ = 35 years showing probability distributions.

...And we're back.

The results of this calculation are shown in Figure 4a and b. For the Gutenberg-Richter distribution, the likelihood of observing an LANF earthquake on any fault over a given observation period is quite high. For example, $P(M,t)$ for $M \, \ge 6$ and $t$ = 35 years is about 0.85, and for the smaller events is quite close to 1. This high likelihood suggests that given the model assumptions, we should expect to find such an earthquake in the focal mechanism catalogs, though because many $M6$ events are not surface-breaking [Hecker et al., 2013], it might be difficult to unambiguously determine whether the high- or low-angle nodal plane slipped. For $M \, \ge 6.5$, the probability of observing an LANF earthquake is about 0.5, and the nodal plane ambiguity should be much less ({\it e.g.}, due to surface scarps or directivity effects). The results for the characteristic distribution are lower than the results for the GR distribution for smaller events and higher for larger events, similar to the patterns seen in results for individual faults. $P(M\ge5.5,t=35)$ through $P(M\ge6.5,t=35)$ are all close, about 0.4--0.5.

3.3 Bayesian adjustments of LANF earthquake likelihood

Because the earthquake focal mechanism catalog is much shorter than the repeat time for moderate to large earthquakes on typical normal faults with mm yr$^{-1}$ slip rates, catalog searches yielding no results for a particular class of events cannot be definitive evidence that they do not occur. Nevertheless, the absence of observations does provide some evidence against their existence. Through Bayes' Theorem, we can use the probability of observing an event (i.e., $P(M,t)$) to calculate the likelihood that LANFs are active given the negative outcome of catalog searches. In this manner, Bayes' Theorem gives an adjusted, posterior likelihood for a given prior likelihood that LANFs are capable of generating large earthquakes. Different priors may result from different evidence or assumptions, and are not likely to be constant through time or among all researchers. We do not choose a specific prior for LANF activity; rather, we calculate the posteriors over the range of prior probabilities from 0 (meaning no probability that LANFs are seismically active) to 1 (meaning LANFs are absolutely seismically active). Here $P(A)$ represents the prior probability for LANF seismic activity, and $P(O)$ is the probability of a positive test result (observation of an LANF earthquake in a catalog search). The symbol '$\sim$' indicates not, so $P(\sim A)$ is the probability that LANFs are inactive; $ P(\sim A) = 1 - P(A)$. The results of this study give the us the probability of observing or not observing an LANF event given LANF seismic activity, $P(O|A)$ and $P(\sim O | A) = 1 - P(O|A)$, respectively. $P(O | \sim A)$ is the probability of observing a 'false positive', the incorrect identification of an LANF event, when in fact LANFs are not active. The posterior $P(A |\sim O)$ is the likelihood that LANFs can generate large earthquakes given that no LANF events have been observed, through Bayes' Theorem

let's code this up too

'''priors for lanf activity'''
p_A = np.linspace(0, 1)

def no_obs(prior, false_pos, accur):
    p_A = prior
    p_nA = 1 - p_A
    
    p_OgnA = false_pos
    p_nOgnA = 1 - p_OgnA
    
    p_OgA = accur
    p_nOgA = 1 - p_OgA
    
    return p_nOgA * p_A / ( p_nOgA * p_A + p_nOgnA * p_nA)

plot(p_A, no_obs(p_A, 0.01, 0.02), 'g', label='P(O|A) = 0.1')
plot(p_A, no_obs(p_A, 0.01, 0.5), 'b', label='P(O|A) = 0.5')
plot(p_A, no_obs(p_A, 0.01, 0.8), 'r', label='P(O|A) = 0.8')

legend(loc='upper left')


xlabel('prior likelihood for LANF seismicity, P(A)')
ylabel('posterior likelihood, P(A|~O)')

show()

Figure 5: Prior likelihood for LANF seismicity $P(A)$ and posterior likelihood $P(A|\sim O)$ given no observed earthquakes. $P(O|A)$ is the likelihood of observing an earthquake given activity on all LANFs.

Figure 5 shows $P(A|\sim O)$ for $P(A) \in [0,1]$, using values for $P(O|A)$ of 0.1, 0.5, and 0.8, and a likelihood of false positives $P(O|\sim A)= 0.01$. The likelihood of LANF seismicity decreases appreciably given a moderate $P(O|A)$, but does not decrease to zero. Low values of $P(O|A)$ yield posteriors that are almost unchanged from the priors; in other words, the fact that no LANF events have occurred does not change the prior assumptions that LANF events are not expected to occur. Additionally, for strong priors with values very close to 0 or 1, the posteriors are much closer to the priors, which is to sat that it takes much more evidence to change strongly-held positions. For less strongly-held prior assumptions, the posterior probability that the LANFs are active is reduced compared to what the prior assumptions are. For example, in the case of a prior of 0.5 (meaning that there is a 50% chance that LANFs can generate large earthquakes), then if $P(O|A)=0.1$, then the posterior likelihood for LANF seismicity drops to $\approx 0.48$ (Figure 5). In this case, the fact that a catalog search results in no identified LANF earthquake is non-informative. If $P(O|A)=0.5$, the posterior drops to $\approx 0.34$, a moderate reduction. On the other hand, if $P(O|A)=0.8$, the posterior is $\approx 0.17$.

Discussion and Conclusions

Our compilation of all known potentially active LANFs shows that they are fairly uncommon structures, yet they still may be found in areas currently undergoing extension. Almost all major continental extensional regions are represented; notably, narrow, linear continental rifts, such as the East African and Rio Grande rifts, do not seem to contain active LANFs. This compilation may serve as a point of comparison for different characteristics of active normal faults or LANF geometry, or as a reference for any further study of these structures.

Any value for $P(M,t)$ we calculate here is a maximum estimate. As discussed above, part of the reason that $P(M,t)$ is the maximum is due to the fact that we assume declustered earthquake catalogs. Additionally, $P(M,t)$ is a maximum estimate as we assume that all LANFs in this study are seismically active throughout the upper crust at surface dip angles. It is quite possible that some of these faults are not tectonically active at all. It is also possible that some or all of these detachments may be seismically active but at dip angles $\ge30^\circ$. For example, the Ca~nada David detachment in Mexico may dip very steeply at seismogenic depths [Fletcher and Spelz, 2009]. Some of these may also be aseismic; the Alto-Tiberina fault appears to be creeping for much of its down-dip extent [Hreinsdóttir and Bennett, 2009], and the neighboring Zuccale inactive LANF has fault gouge suggestive of creep [Collettini and Holdsworth, 2004]. If any individual fault is not seismically active at low angles, this reduces the total $P(M, t)$ for all events.

$P(M \ge 6.5, t = 35) \approx 0.5$ for either frequency-–magnitude distribution, and is a good reference value as events in this range are likely to be surface–-breaking, and which would resolve the slip plane ambiguity inherent in earthquake focal mechanisms [Hecker et al., 2013]. Given the fact that no significant LANF earthquakes have been definitively documented, this probability of an LANF earthquake occurring during a 35 year time window results in a lowering of any prior assumption that LANFs are active, as long as that prior assumption is not a strongly-held position. The magnitude of the decrease depends on the prior likelihood, and the decrease is at most $\sim 15 \%$ (from 0.5 to 0.35). This means that the current catalog length is much too short to be used as strong evidence against LANF seismicity. $P(M\ge6.5, t=100)$ is near 0.8 for both GR and characteristic distributions; this value more strongly reduces the likelihood of LANF seismicity, yet still does not yield a definitive negative conclusion. Therefore, results of studies analyzing the dip distribution of earthquakes on continental normal faults [Jackson, 1987; Collettini and Sibson, 2001] should be interpreted as informative but not conclusive. Furthermore, alternative mechanisms for LANF occurrence or slip such as aseismic creep [e.g., Collettini, 2011; Hreinsdóttir and Bennett, 2009], isostatic flexure [e.g., Wernicke and Axen, 1988], or [Wernicke, 1995] need not be invoked to explain the lack of observed seismicity, though these mechanisms may indeed be valid or well supported by other observations.

Acknowledgements

We thank Jon Spencer for a stimulating discussion that became the impetus for this study. Mike Taylor and Kurt Sundell provided valuable comments on a draft of the manuscript.

References

Abers, G. A. (2001), Evidence for seismogenic normal faults at shallow dips in continental rifts, Geological Society, London, Special Publications, 187(1), 305–318.

Aki, K., and P. G. Richards (2002), Quantitative seismology: Theory and methods, University Science Books.

Axen, G. J. (2004), Mechanics of low-angle normal faults, in Rheology and Deformation of the Lithosphere at Continental Margins, edited by G. D. Karner et al., pp. 46–91, Columbia Univ. Press, New York.

Axen, G. J., and J. M. Bartley (1997), Field tests of rolling hinges: Existence, mechanical types, and implications for extensional tectonics, Journal of Geophysical Research, 102(B9), 20,515–20.

Axen, G. J., J. M. Fletcher, E. Cowgill, M. Murphy, P. Kapp, I. MacMillan, E. Ramos-Vela ́zquez, and J. Aranda-Gomez (1999), Range-front fault scarps of the sierra el mayor, baja california: Formed above an active low-angle normal fault?, Geology, 27(3), 247–250.

Bernard, P., P. Briole, B. Meyer, H. Lyon-Caen, J.-M. Gomez, C. Tiberi, C. Berge, R. Cattin, D. Hatzfeld, C. Lachet, et al. (1997), The ms= 6.2, june 15, 1995 aigion earthquake (greece): evidence for low angle normal faulting in the corinth rift, Journal of Seismology, 1(2), 131–150.

Bird, P., and Y. Y. Kagan (2004), Plate-tectonic analysis of shallow seismicity: Apparent boundary width, beta, corner magnitude, coupled litho- sphere thickness, and coupling in seven tectonic settings, Bulletin of the Seismological Society of America, 94(6), 2380–2399.

Buck, W. R. (1991), Modes of continental lithospheric extension, Journal of Geophysical Research: Solid Earth (1978–2012), 96(B12), 20,161– 20,178.

Collettini, C. (2011), The mechanical paradox of low-angle normal faults: Current understanding and open questions, Tectonophysics, 510(3), 253–268.

Collettini, C., and R. Holdsworth (2004), Fault zone weakening and character of slip along low-angle normal faults: insights from the Zuccale fault, Elba, Italy, Journal of the Geological Society, 161(6), 1039–1051.

Collettini, C., and R. H. Sibson (2001), Normal faults, normal friction?, Geology, 29(10), 927–930.

Doser, D. I. (1987), The Ancash, Peru, earthquake of 1946 November 10: Evidence for low-angle normal faulting in the high Andes of northern Peru, Geophysical Journal of the Royal Astronomical Society, 91(1), 57–71.

Fletcher, J. M., and R. M. Spelz (2009), Patterns of Quaternary deformation and rupture propagation associated with an active low-angle normal fault, Laguna Salada, Mexico: Evidence of a rolling hinge?, Geosphere, 5(4), 385–407.

Hayman, N. W., J. R. Knott, D. S. Cowan, E. Nemser, and A. M. Sarna- Wojcicki (2003), Quaternary low-angle slip on detachment faults in Death Valley, California, Geology, 31(4), 343–346.

Hecker, S., N. Abrahamson, and K. Wooddell (2013), Variability of displacement at a point: Implications for earthquake-size distribution and rupture hazard on faults, Bulletin of the Seismological Society of America, 103(2A), 651–674.

Howard, K. A., and B. E. John (1987), Crustal extension along a rooted system of imbricate low-angle faults: Colorado River extensional corridor, California and Arizona, Geological Society, London, Special Publications, 28(1), 299–311.

Hreinsdottir, S., and R. A. Bennett (2009), Active aseismic creep on the alto tiberina low-angle normal fault, italy, Geology, 37(8), 683–686.

Jackson, J. (1987), Active normal faulting and crustal extension, Geological Society, London, Special Publications, 28(1), 3–17.

Kagan, Y. Y. (2003), Accuracy of modern global earthquake catalogs, Physics of the Earth and Planetary Interiors, 135(2), 173–209.

Kapp, J. L., T. M. Harrison, P. Kapp, M. Grove, O. M. Lovera, and D. Lin (2005), Nyainqentanglha Shan: a window into the tectonic, thermal, and geochemical evolution of the Lhasa block, southern Tibet, Journal of Geophysical Research: Solid Earth (1978–2012), 110(B8).

Lister, G., M. Etheridge, and P. Symonds (1986), Detachment faulting and the evolution of passive continental margins, Geology, 14(3), 246–250.

McKinney, W. (2010), Data structures for statistical computing in python, in Proceedings of the 9th Python in Science Conference, edited by S. van der Walt and J. Millman, pp. 51 – 56.

Oliphant, T. E. (2007), Python for scientific computing, Computing in Science & Engineering, 9(3), 10–20.

Perez, F., and B. E. Granger (2007), IPython: a System for Interactive Scientific Computing, Comput. Sci. Eng., 9(3), 21–29.

Sibson, R. H. (1985), A note on fault reactivation, Journal of Structural Geology, 7(6), 751–754.

Styron, R. H., M. H. Taylor, K. E. Sundell, D. F. Stockli, J. A. Oalmann, A. Moller, A. T. McCallister, D. Liu, and L. Ding (2013), Miocene initiation and acceleration of extension in the South Lunggar rift, western Tibet: Evolution of an active detachment system from structural map- ping and (U-Th)/He thermochronology, Tectonics, 32(4), 880–907, doi: 10.1002/tect.20053.

Sundell, K. E., M. H. Taylor, R. H. Styron, D. F. Stockli, P. Kapp, C. Hager, D. Liu, and L. Ding (2013), Evidence for constriction and Pliocene acceleration of east-west extension in the North Lunggar rift region of west central Tibet, Tectonics, 32(5), 1454–1479, doi:10.1002/tect.20086.

Varoquaux, G., and O. Grisel (2009), Joblib: Running Python function as pipeline jobs, http://pythonhosted.org/joblib/.

Wernicke, B. (1995), Low-angle normal faults and seismicity: A review, Journal of Geophysical Research, 100(B10), 20,159–20.

Wernicke, B. (2009), The detachment era (1977–1982) and its role in revolutionizing continental tectonics, Geological Society, London, Special Publications, 321(1), 1–8.

Wernicke, B., and G. J. Axen (1988), On the role of isostasy in the evolution of normal fault systems, Geology, 16(9), 848–851.

Wesnousky, S. G. (2008), Displacement and geometrical characteristics of earthquake surface ruptures: Issues and implications for seismic-hazard analysis and the process of earthquake rupture, Bul- letin of the Seismological Society of America, 98(4), 1609–1632, doi: 10.1785/0120070111.

research active faults GIS Tibet geology programming python low-angle normal faults earthquakes