API

ILSI.ilsi.Michael1984_inversion(strikes, dips, rakes, Tarantola_kwargs={}, return_eigen=True, return_stats=False)

Linear inversion described in Michael 1984.

This method assumes:

• The tectonic stress field is uniform.

• Wallace-Bott hypothesis: The slip vector points in the same direction as shear stress on the fault.

• The resolved shear stress magnitude is constant on all faults.

The parameters we invert for are the directions of the three principal stresses and the shape ratio. Because this inversion does not aim at infering the absolute stress values, we only consider the deviatoric stress tensor, therefore Trace(sigma) = 0. Furthermore, we cannot determine the norm of the stress tensor, therefore sum sigma**2 = 1. Each iteration of this inversion scheme is a linear inversion. N.B.: This routine is written assuming outward footwall normals and slip vectors of the hanging wall w.r.t. the footwall. Therefore, the stress tensor sign convention is compression negative.

Parameters

• strikes (list or numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• dips (list or numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• rakes (list or numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• Tarantola_kwargs (Dictionary, default to {}) – If not None, should contain key word arguments for the Tarantola and Valette inversion. An empty dictionary uses the default values in Tarantola_Valette. If None, uses the Moore-Penrose inverse.

• return_eigen (boolean, default to True) – If True, returns the eigendecomposition of the inverted stress tensor in addition to returning the stress tensor.

• return_stats (boolean, default to True) – If True, the posterior data and model parameter distributions estimated from the Tarantola and Valette formula (cf. Tarantola_Valette routine).

Returns

output –

• output[“stress_tensor”]: (3, 3) numpy.ndarray

  The inverted stress tensor in the (north, west, up) coordinate system.

• output[“principal_stresses”]: (3,) numpy.ndarray, optional

  The three eigenvalues of the stress tensor, ordered from most compressive (sigma1) to least compressive (sigma3). Returned if return_eigen is True.

• output[“predicted_shear_stress”]: (n_earthquakes, 3) numpy.ndarray

  The shear tractions resolved on the fault planes described by strikes, dips and rakes computed with the inverted stress tensor.

• output[“principal_directions”]: (3, 3) numpy.ndarray, optional

  The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i]. Returned if return_eigen is True.

• output[“C_m_posterior”]: (5, 5) numpy.ndarray, optional

  Posterior covariance of the model parameter distribution estimated from the Tarantola and Valette formula. Returned if return_stats is True.

• output[“C_d_posterior”]: (3 x n_earthquakes, 3 x n_earthquakes) numpy.ndarray, optional

  Posterior covariance of the data distribution estimated from the Tarantola and Valette formula. Returned if return_stats is True.

Return type

dict {str: numpy.ndarray}

ILSI.ilsi.Tarantola_Valette(G, data, C_d=None, C_d_inv=None, C_m=None, C_m_inv=None, m_prior=None, inversion_space='model_space')

Returns Tarantola’s and Valette’s least square solution for a given linear operator G and observation vector data. If the covariance matrices of the observations and of the model parameters are not known, we assume them to be identity. The inversion can be performed either in the data space or in the model space.

Parameters

• G ((n, m) numpy.ndarray) – The linear operator projecting elements of the model space m onto the data space: d = G.m n is the dimension of the data space, m is the dimension of the model space.

• data ((3k,) or (3k, 1) or (k, 3) numpy.ndarray) – Vector of observations. k is the number of focal mechanisms. data is reshaped to (n=3k, 1) before the inversion.

• C_d ((n, n) numpy.ndarray, default to None) – Covariance matrix of the observations. It quantifies the errors in the observations and propagates them in the inversion to give more weight to the observations with low errors. If None, then C_d is filled with zeros (assume no error in data).

• C_m ((m, m) numpy.ndarray, default to None) – Covariance matrix of the model parameters. It quantifies the errors in the model parameters and propagates them in the inversion to determine the range of acceptable model parameters for a given set of observations. If None, then C_m is identity.

• m_prior ((m,) or (m, 1) numpy.ndarray, default to None) – If one already has a rough estimate of what the model parameters are, then m_prior should be filled with this estimate. If None, m_prior is set to zero.

Returns

• m_inv ((m, 1) numpy.ndarray) – The inverted model parameters.

• C_m_posterior ((5, 5) numpy.ndarray) – Posterior covariance of the model parameter distribution.

• C_d_posterior ((3 x n_earthquakes, 3 x n_earthquakes) numpy.ndarray) – Posterior covariance of the data distribution.

ILSI.ilsi.compute_instability_parameter(principal_directions, R, friction, strike_1, dip_1, rake_1, strike_2, dip_2, rake_2, return_fault_planes=False, signed_instability=False)

Compute the instability parameter as introduced by Lund and Slunga 1999, re-used by Vavrycuk 2013-2014 and modified by Beauce 2022. For a given stress field characterized by the principal stress directions and shape ratio R=(sig1-sig2)/(sig1-sig3), and for a given rock friction, this routine computes an instability parameter based on the Mohr-Coulomb failure criterion to determine which of the two nodal planes of a focal mechanism solution is more likely to be the fault plane. Beauce 2022 includes the sign of the dot product between the shear stress and the slip vector on the fault. This instability ranges from -1 to +1, instead of from 0 to +1.

Parameters

• principal_directions ((3, 3) numpy.ndarray, float) – The three eigenvectors of the reference stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i]

• R (float) – Shape ratio of the reference stress tensor.

• strikes_1 (list or numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• dips_1 (list or numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• rakes_1 (list or numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• strikes_2 (list or numpy.ndarray, float) – The strike of nodal planes 2, angle between north and the fault’s horizontal (0-360).

• dips_2 (list or numpy.ndarray, float) – The dip of nodal planes 2, angle between the horizontal plane and the fault plane (0-90).

• rakes_2 (list or numpy.ndarray, float) – The rake of nodal planes 2, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• return_fault_planes (boolean, default to False) – If True, return the strikes, dips, rakes of the selected fault planes.

• signed_instability (boolean, default to False) – If True, the instability parameter ranges from -1 to +1. Negative values mean that the predicted and observed slip have opposite directions. If False, the instability parameter is the one defined in Vavrycuk 2013, 2014.

Returns

• instability_parameter ((n_earthquakes, 2) numpy.ndarray) – The instability parameter as defined in Beauce 2022 for the two nodal planes of each focal mechanism datum.

• strikes (list or numpy.ndarray, float, optional) – Strikes of the fault planes with largest instability. Only provided if return_fault_planes=True.

• dips (list or numpy.ndarray, float, optional) – Dips of the fault planes with largest instability. Only provided if return_fault_planes=True.

• rakes (list or numpy.ndarray, float, optional) – Rakes of the fault planes with largest instability. Only provided if return_fault_planes=True.

ILSI.ilsi.find_optimal_friction(strikes_1, dips_1, rakes_1, strikes_2, dips_2, rakes_2, principal_directions, R, friction_min=0.1, friction_max=0.8, friction_step=0.05, signed_instability=False)

Find the friction that maximizes the instability parameter I based on V. Vavrycuk 2013,2014 and B. Lund and R. Slunga 1999.

Parameters

• strikes_1 (list or numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• dips_1 (list or numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• rakes_1 (list or numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• strikes_2 (list or numpy.ndarray, float) – The strike of nodal planes 2, angle between north and the fault’s horizontal (0-360).

• dips_2 (list or numpy.ndarray, float) – The dip of nodal planes 2, angle between the horizontal plane and the fault plane (0-90).

• rakes_2 (list or numpy.ndarray, float) – The rake of nodal planes 2, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• principal_directions ((3, 3) numpy.ndarray, float) – The three eigenvectors of the reference stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i].

• R (float) – Shape ratio of the reference stress tensor.

• friction_min (float, default to 0.1) – Lower bound of explored friction values.

• friction_max (float, default to 0.8) – Upper bound of explored friction values.

• friction_step (float, default to 0.05) – Step employed in the grid search of the friction value that maximizes the instability parameter.

• signed_instability (boolean, default to False) – If True, the instability parameter ranges from -1 to +1. Negative values mean that the predicted and observed slip have opposite directions. If False, the instability parameter is the one defined in Vavrycuk 2013, 2014.

Returns

optimal_friction – The friction value that maximizes the mean instability parameter.

Return type

float

ILSI.ilsi.find_optimal_friction_one_set(strikes_1, dips_1, rakes_1, principal_directions, R, friction_min=0.2, friction_max=0.8, friction_step=0.05, signed_instability=False)

Find the friction that maximizes the instability parameter I based on V. Vavrycuk 2013,2014 and B. Lund and R. Slunga 1999.

Parameters

• strikes_1 (list or numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• dips_1 (list or numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• rakes_1 (list or numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• principal_directions ((3, 3) numpy.ndarray, float) – The three eigenvectors of the reference stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i].

• R (float) – Shape ratio of the reference stress tensor.

• friction_min (float, default to 0.1) – Lower bound of explored friction values.

• friction_max (float, default to 0.8) – Upper bound of explored friction values.

• friction_step (float, default to 0.05) – Step employed in the grid search of the friction value that maximizes the instability parameter.

• signed_instability (boolean, default to False) – If True, the instability parameter ranges from -1 to +1. Negative values mean that the predicted and observed slip have opposite directions. If False, the instability parameter is the one defined in Vavrycuk 2013, 2014.

Returns

optimal_friction – The friction value that maximizes the mean instability parameter.

Return type

float

ILSI.ilsi.forward_model(n_)

Build the forward modeling matrix G given a collection of fault normals.

Parameters

n ((n_earthquakes, 3) numpy.ndarray) – The i-th row n_ are the components of the i-th fault normal in the (north, west, south) coordinate system.

Returns

G – The forward modeling matrix giving the slip (shear stress) directions on the faults characterized by n_, given the 5 elements of the deviatoric stress tensor.

Return type

(3 x n_earthquakes, 5) numpy.ndarray

ILSI.ilsi.inversion_bootstrap(strikes, dips, rakes, n_random_selections=1, n_resamplings=100, variable_shear=True, max_n_iterations=300, shear_update_atol=1e-05, Tarantola_kwargs={})

Inverts one set of focal mechanisms without seeking which nodal planes are more likely to be the fault planes. Performs bootstrap resampling of the data set to return an ensemble of solutions.

Parameters

• strikes (list or numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• dips (list or numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• rakes (list or numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• n_random_selections (integer, default to 5) – Number of random selections of subsets of nodal planes on which the stress inversion is run. The final stress tensor is averaged over the n_random_selections solutions.

• n_resamplings (integer, default to 100) – Number of times the data set is resampled following the bootstrapping method (sampling with replacement). n_resamplings stress tensors are returned, allowing to estimate uncertainties from the distribution of solutions.

• shear_update_atol (float, default to 1e-5) – Convergence criterion on the shear stress magnitude updates. Convergence is reached when the RMS difference between two estimates of shear stress magnitudes falls below that threshold.

• max_n_iterations (integer, default to 300) – The maximum number of iterations if shear stress magnitude update does not fall below shear_update_atol.

• variable_shear (boolean, default to True) – If True, use the iterative linear method described in Beauce et al. 2022, else use the classic linear method due to Michael 1984.

• Tarantola_kwargs (Dictionary, default to {}) – If not None, should contain key word arguments for the Tarantola and Valette inversion. An empty dictionary uses the default values in Tarantola_Valette. If None, uses the Moore-Penrose inverse.

Returns

output –

• output[“boot_stress_tensor”]: (n_resamplings, 3, 3) numpy.ndarray

  The inverted stress tensor in the (north, west, up) coordinate system.

• output[“boot_principal_stresses”]: (n_resamplings, 3) numpy.ndarray

  The three eigenvalues of the stress tensor, ordered from most compressive (sigma1) to least compressive (sigma3).

• output[“boot_principal_directions”]: (n_resamplings, 3, 3) numpy.ndarray

  The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i for the b-th bootstrap replica is given by: boot_principal_directions[b, :, i].

Return type

dict {str: numpy.ndarray}

ILSI.ilsi.inversion_bootstrap_instability(principal_directions, R, strikes, dips, rakes, friction_coefficient, n_resamplings=100, n_stress_iter=10, stress_tensor_update_atol=1e-05, Tarantola_kwargs={}, variable_shear=True, max_n_iterations=300, shear_update_atol=1e-05, signed_instability=False, weighted=False, n_threads=1, parallel=False)

Invert one set of focal mechanisms with the instability parameter to seek which nodal planes are more likely to be the fault planes (cf. B. Lund and R. Slunga 1999, V. Vavrycuk 2013,2014). Performs bootstrap resampling of the data set to return an ensemble of solutions.

Use a previously determined stress tensor (e.g. the output of inversion_one_set_instability) described by its principal stress directions and shape ratio as the prior model in the Tarantola and Valette formula. In general, you can keep the default parameter values, except for n_resamplings which depends on the time you can afford spending.

Parameters

• principal_directions ((3, 3) numpy.ndarray, float) – The three eigenvectors of the reference stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i].

• R (float) – Shape ratio of the reference stress tensor.

• friction_coefficient (float) – Value of the friction coefficient used in the instability parameter. This can be the value output by inversion_one_set_instability.

• strikes (list or numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• dips (list or numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• rakes (list or numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• n_stress_iter (integer, default to 10) – Maximum number of iterations to seek for the best fault planes. See Beauce et al. 2022 for explanations.

• stress_tensor_update_atol (float, default to 1.e-4) – If the RMS difference of the stress tensors between two iterations fall below this threshold, convergence has been reached.

• shear_update_atol (float, default to 1e-5) – Convergence criterion on the shear stress magnitude updates. Convergence is reached when the RMS difference between two estimates of shear stress magnitudes falls below that threshold.

• signed_instability (boolean, default to False) – If True, the instability parameter ranges from -1 to +1. Negative values mean that the predicted and observed slip have opposite directions. If False, the instability parameter is the one defined in Vavrycuk 2013, 2014.

• max_n_iterations (integer, default to 300) – The maximum number of iterations if shear stress magnitude update does not fall below shear_update_atol.

• variable_shear (boolean, default to True) – If True, use the iterative linear method described in Beauce et al. 2022, else use the classic linear method due to Michael 1984.

• Tarantola_kwargs (Dictionary, default to {}) – If not None, should contain key word arguments for the Tarantola and Valette inversion. An empty dictionary uses the default values in Tarantola_Valette. If None, uses the Moore-Penrose inverse.

• weighted (boolean, default to False) –

  This option is exploratory. If True,

  1. More weight is given to the fault planes that are clearly more unstable than their auxiliary counterpart in the stress field estimated at iteration t-1

  2. Randomly mixes the set of fault planes at iterations t-1 and t, giving larger probability to the planes belonging to the set that produced the larger instability.

  This option can be interesting for reaching convergence on data sets of bad quality.

• n_threads (scalar int, optional) – Default to n_threads=1. If different from 1, the task is parallelized across n_threads threads. If n_threads is 0, None or “all”, use all available CPUs.

Returns

output –

• output[“boot_stress_tensor”]: (n_resamplings, 3, 3) numpy.ndarray

  The inverted stress tensor in the (north, west, up) coordinate system.

• output[“boot_principal_stresses”]: (n_resamplings, 3) numpy.ndarray

  The three eigenvalues of the stress tensor, ordered from most compressive (sigma1) to least compressive (sigma3).

• output[“boot_principal_directions”]: (n_resamplings, 3, 3) numpy.ndarray

  The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i for the b-th bootstrap replica is given by: boot_principal_directions[b, :, i].

Return type

dict {str: numpy.ndarray}

ILSI.ilsi.inversion_jackknife(jack_strikes, jack_dips, jack_rakes, n_random_selections=1, n_resamplings=100, max_n_iterations=300, shear_update_atol=1e-05, variable_shear=True, Tarantola_kwargs={}, bootstrap_events=False)

This routine was tailored for one of my application, but it can be of interest to others. Each earthquake comes with an ensemble of focal mechanism solutions that were obtained by resampling the set of seismic stations used in the focal mechanism inversion. The resampling was done with the delete-k-jackknife method, hence the name of the routine. This routine randomly samples focal mechanisms from these ensembles and runs the stress inversion. This is a way of propagating the focal mechanism uncertainties into the stress inversion. In this routine we do not seek which nodal planes are more likely to be the fault planes.

Parameters

• jack_strikes ((n_earthquakes, n_jackknifes) numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• jack_dips ((n_earthquakes, n_jackknifes) numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• jack_rakes ((n_earthquakes, n_jackknifes) numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• n_random_selections (integer, default to 20) – Number of random selections of subsets of nodal planes on which the stress inversion is run. The final stress tensor is averaged over the n_random_selections solutions.

• n_resamplings (integer, default to 100) – Number of times the data set is resampled from the ensembles of focal mechanism solutions available for each earthquake. n_resamplings stress tensors are returned, allowing to estimate uncertainties from the distribution of solutions.

• bootstrap_events (boolean, default to False) – If True, the resampling is also done accross earthquakes, following the bootstrapping method.

• shear_update_atol (float, default to 1e-5) – Convergence criterion on the shear stress magnitude updates. Convergence is reached when the RMS difference between two estimates of shear stress magnitudes falls below that threshold.

• max_n_iterations (integer, default to 300) – The maximum number of iterations if shear stress magnitude update does not fall below shear_update_atol.

• variable_shear (boolean, default to True) – If True, use the iterative linear method described in Beauce et al. 2022, else use the classic linear method due to Michael 1984.

• Tarantola_kwargs (Dictionary, default to {}) – If not None, should contain key word arguments for the Tarantola and Valette inversion. An empty dictionary uses the default values in Tarantola_Valette. If None, uses the Moore-Penrose inverse.

Returns

output –

• output[“jack_stress_tensor”]: (n_resamplings, 3, 3) numpy.ndarray

  The inverted stress tensor in the (north, west, up) coordinate system.

• output[“jack_principal_stresses”]: (n_resamplings, 3) numpy.ndarray

  The three eigenvalues of the stress tensor, ordered from most compressive (sigma1) to least compressive (sigma3).

• output[“jack_principal_directions”]: (n_resamplings, 3, 3) numpy.ndarray

  The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i for the b-th jackknife replica is given by: jack_principal_directions[b, :, i].

Return type

dict {str: numpy.ndarray}

ILSI.ilsi.inversion_jackknife_instability(principal_directions, R, jack_strikes, jack_dips, jack_rakes, friction_coefficient, n_resamplings=100, n_stress_iter=10, stress_tensor_update_atol=0.0001, Tarantola_kwargs={}, bootstrap_events=False, n_earthquakes=None, variable_shear=True, max_n_iterations=300, shear_update_atol=1e-05, signed_instability=False, weighted=False, n_threads=1, parallel=False)

This routine was tailored for one of my application, but it can be of interest to others. Each earthquake comes with an ensemble of focal mechanism solutions that were obtained by resampling the set of seismic stations used in the focal mechanism inversion. The resampling was done with the delete-k-jackknife method, hence the name of the routine. This routine randomly samples focal mechanisms from these ensembles and runs the stress inversion. This is a way of propagating the focal mechanism uncertainties into the stress inversion. Use the instability parameter to seek which nodal planes are more likely to be the fault planes (cf. B. Lund and R. Slunga 1999, V. Vavrycuk 2013,2014).

Use a previously determined stress tensor (e.g. the output of inversion_one_set_instability) described by its principal stress directions and shape ratio as the prior model in the Tarantola and Valette formula. In general, you can keep the default parameter values, except for n_resamplings which depends on the time you can afford spending.

Parameters

• principal_directions ((3, 3) numpy.ndarray, float) – The three eigenvectors of the reference stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i].

• R (scalar float) – Shape ratio of the reference stress tensor.

• friction_coefficient (scalar float) – Friction value used in the instability parameter. This can be the value output by inversion_one_set_instability.

• jack_strikes ((n_earthquakes, n_jackknifes) numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• jack_dips ((n_earthquakes, n_jackknifes) numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• jack_rakes ((n_earthquakes, n_jackknifes) numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• n_stress_iter (integer, default to 10) – Maximum number of iterations to seek for the best fault planes. See Beauce et al. 2022 for explanations.

• stress_tensor_update_atol (float, default to 1.e-4) – If the RMS difference of the stress tensors between two iterations fall below this threshold, convergence has been reached.

• shear_update_atol (float, default to 1e-5) – Convergence criterion on the shear stress magnitude updates. Convergence is reached when the RMS difference between two estimates of shear stress magnitudes falls below that threshold.

• signed_instability (boolean, default to False) – If True, the instability parameter ranges from -1 to +1. Negative values mean that the predicted and observed slip have opposite directions. If False, the instability parameter is the one defined in Vavrycuk 2013, 2014.

• max_n_iterations (integer, default to 300) – The maximum number of iterations if shear stress magnitude update does not fall below shear_update_atol.

• variable_shear (boolean, default to True) – If True, use the iterative linear method described in Beauce et al. 2022, else use the classic linear method due to Michael 1984.

• Tarantola_kwargs (Dictionary, default to {}) – If not None, should contain key word arguments for the Tarantola and Valette inversion. An empty dictionary uses the default values in Tarantola_Valette. If None, uses the Moore-Penrose inverse.

• bootstrap_events (boolean, default to False) – If True, the resampling is also done accross earthquakes, following the bootstrapping method.

• weighted (boolean, default to False) –

  This option is exploratory. If True:

  1. More weight is given to the fault planes that are clearly more unstable than their auxiliary counterpart in the stress field estimated at iteration t-1

  2. Randomly mixes the set of fault planes at iterations t-1 and t, giving larger probability to the planes belonging to the set that produced the larger instability.

  This option can be interesting for reaching convergence on data sets of bad quality.

• n_threads (scalar int, optional) – Default to n_threads=1. If different from 1, the task is parallelized across n_threads threads. If n_threads is 0, None or “all”, use all available CPUs.

Returns

output –

• output[“jack_stress_tensor”]: (n_resamplings, 3, 3) numpy.ndarray

  The inverted stress tensor in the (north, west, up) coordinate system.

• output[“jack_principal_stresses”]: (n_resamplings, 3) numpy.ndarray

  The three eigenvalues of the stress tensor, ordered from most compressive (sigma1) to least compressive (sigma3).

• output[“jack_principal_directions”]: (n_resamplings, 3, 3) numpy.ndarray

  The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i for the b-th jackknife replica is given by: jack_principal_directions[b, :, i].

Return type

dict {str: numpy.ndarray}

ILSI.ilsi.inversion_one_set(strikes, dips, rakes, n_random_selections=20, max_n_iterations=300, shear_update_atol=1e-05, variable_shear=True, Tarantola_kwargs={}, return_eigen=True, return_stats=False)

Invert one set of focal mechanisms without seeking which nodal planes are more likely to be the fault planes.

Parameters

• strikes (list or numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• dips (list or numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• rakes (list or numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• n_random_selections (integer, default to 20) – Number of random selections of subsets of nodal planes on which the stress inversion is run. The final stress tensor is averaged over the n_random_selections solutions.

• shear_update_atol (float, default to 1e-5) – Convergence criterion on the shear stress magnitude updates. Convergence is reached when the RMS difference between two estimates of shear stress magnitudes falls below that threshold.

• max_n_iterations (integer, default to 300) – The maximum number of iterations if shear stress magnitude update does not fall below shear_update_atol.

• variable_shear (boolean, default to True) – If True, use the iterative linear method described in Beauce et al. 2022, else use the classic linear method due to Michael 1984.

• Tarantola_kwargs (Dictionary, default to {}) – If not None, should contain key word arguments for the Tarantola and Valette inversion. An empty dictionary uses the default values in Tarantola_Valette. If None, uses the Moore-Penrose inverse.

• return_stats (boolean, default to True) – If True, the posterior data and model parameter distributions estimated from the Tarantola and Valette formula (cf. Tarantola_Valette routine).

Returns

output –

• output[“stress_tensor”]: (3, 3) numpy.ndarray

  The inverted stress tensor in the (north, west, up) coordinate system.

• output[“principal_stresses”]: (3,) numpy.ndarray, optional

  The three eigenvalues of the stress tensor, ordered from most compressive (sigma1) to least compressive (sigma3). Returned if return_eigen is True.

• output[“principal_directions”]: (3, 3) numpy.ndarray, optional

  The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i]. Returned if return_eigen is True.

• output[“C_m_posterior”]: (5, 5) numpy.ndarray, optional

  Posterior covariance of the model parameter distribution estimated from the Tarantola and Valette formula. Returned if return_stats is True.

• output[“C_d_posterior”]: (3 x n_earthquakes, 3 x n_earthquakes) numpy.ndarray, optional

  Posterior covariance of the data distribution estimated from the Tarantola and Valette formula. Returned if return_stats is True.

Return type

dict {str: numpy.ndarray}

ILSI.ilsi.inversion_one_set_instability(strikes, dips, rakes, friction_coefficient=0.6, friction_min=0.1, friction_max=0.8, friction_step=0.05, n_stress_iter=10, n_random_selections=20, stress_tensor_update_atol=0.0001, Tarantola_kwargs={}, max_n_iterations=300, shear_update_atol=1e-05, n_averaging=1, signed_instability=False, verbose=True, variable_shear=True, return_stats=False, weighted=False, plot=False)

Invert one set of focal mechanisms with the instability parameter to seek which nodal planes are more likely to be the fault planes (cf. B. Lund and R. Slunga 1999, V. Vavrycuk 2013,2014). In general, you can keep the default parameter values.

Parameters

• strikes (list or numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• dips (list or numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• rakes (list or numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• friction_coefficient (float or None, default to 0.6) – If not None, the inversion is made assuming a friction coefficient equal to friction_coefficient. If None, the friction coefficient is taken as the one that maximizes instability based on a first approximation of the stress tensor.

• friction_min (float, default to 0.1) – Lower bound of explored friction values.

• friction_max (float, default to 0.8) – Upper bound of explored friction values.

• friction_step (float, default to 0.05) – Step employed in the grid search of the friction value that maximizes the instability parameter.

• n_stress_iter (integer, default to 10) – Maximum number of iterations to seek for the best fault planes. See Beauce et al. 2022 for explanations.

• stress_tensor_update_atol (float, default to 1.e-4) – If the RMS difference of the stress tensors between two iterations fall below this threshold, convergence has been reached.

• n_random_selections (integer, default to 20) – Number of random selections of subsets of nodal planes on which the stress inversion is run. The final stress tensor is averaged over the n_random_selections solutions.

• shear_update_atol (float, default to 1e-5) – Convergence criterion on the shear stress magnitude updates. Convergence is reached when the RMS difference between two estimates of shear stress magnitudes falls below that threshold.

• max_n_iterations (integer, default to 300) – The maximum number of iterations if shear stress magnitude update does not fall below shear_update_atol.

• variable_shear (boolean, default to True) – If True, use the iterative linear method described in Beauce et al. 2022, else use the classic linear method due to Michael 1984.

• n_averaging (integer, default to 1) – The inversion can be sensitive to initial conditions. To improve reproducibility of the results it is good to repeat the inversion several times and average the results. Set n_averaging to ~5 if you can afford the increase in run time.

• signed_instability (boolean, default to False) – If True, the instability parameter ranges from -1 to +1. Negative values mean that the predicted and observed slip have opposite directions. If False, the instability parameter is the one defined in Vavrycuk 2013, 2014.

• Tarantola_kwargs (Dictionary, default to {}) – If not None, should contain key word arguments for the Tarantola and Valette inversion. An empty dictionary uses the default values in Tarantola_Valette. If None, uses the Moore-Penrose inverse.

• return_stats (boolean, default to True) – If True, the posterior data and model parameter distributions estimated from the Tarantola and Valette formula (cf. Tarantola_Valette routine).

• weighted (boolean, default to False) –

  This option is exploratory. If True:

  1. More weight is given to the fault planes that are clearly more unstable than their auxiliary counterpart in the stress field estimated at iteration t-1

  2. Randomly mixes the set of fault planes at iterations t-1 and t, giving larger probability to the planes belonging to the set that produced the larger instability.

  This option can be interesting for reaching convergence on data sets of bad quality.

• plot (boolean, default to False) – If True, plot the set of nodal planes selected at each iteration, and the weight attributed to each of these planes. Can be used with weighted=True to see if the inversion convergences to a well defined set of planes.

• verbose (integer, default to 1) –

  Level of verbosity. 0: No print statements. 1: Print whether the algorithm converged. 2: Print the stress tensor at the end of each fault plane

  selection iteration.

Returns

output –

• output[“stress_tensor”]: (3, 3) numpy.ndarray

  The inverted stress tensor in the (north, west, up) coordinate system.

• output[“friction_coefficient”]: scalar float

  The best friction coefficient determined by the inversion or the input friction coefficient (see friction_coefficient).

• output[“principal_stresses”]: (3,) numpy.ndarray

  The three eigenvalues of the stress tensor, ordered from most compressive (sigma1) to least compressive (sigma3).

• output[“principal_directions”]: (3, 3) numpy.ndarray

  The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i].

• output[“C_m_posterior”]: (5, 5) numpy.ndarray, optional

  Posterior covariance of the model parameter distribution estimated from the Tarantola and Valette formula. Returned if return_stats is True.

• output[“C_d_posterior”]: (3 x n_earthquakes, 3 x n_earthquakes) numpy.ndarray, optional

  Posterior covariance of the data distribution estimated from the Tarantola and Valette formula. Returned if return_stats is True.

Return type

dict {str: numpy.ndarray}

ILSI.ilsi.iterative_linear_si(strikes, dips, rakes, max_n_iterations=300, shear_update_atol=1e-05, Tarantola_kwargs={}, return_eigen=True, return_stats=False)

Iterative stress inversion described in Beauce et al. 2022.

This method assumes:

• The tectonic stress field is uniform.

• Wallace-Bott hypothesis: The slip vector points in the same direction as shear stress on the fault.

The parameters we invert for are the directions of the three principal stresses and the shape ratio. Because this inversion does not aim at infering the absolute stress values, we only consider the deviatoric stress tensor, therefore Trace(sigma) = 0. Furthermore, we cannot determine the norm of the stress tensor, therefore sum sigma**2 = 1. Each iteration of this inversion scheme is a linear inversion. N.B.: This routine is written assuming outward footwall normals and slip vectors of the hanging wall w.r.t. the footwall. Therefore, the stress tensor sign convention is compression negative.

Parameters

• strikes (list or numpy.ndarray, float) – The strike of nodal planes 1, angle between north and the fault’s horizontal (0-360).

• dips (list or numpy.ndarray, float) – The dip of nodal planes 1, angle between the horizontal plane and the fault plane (0-90).

• rakes (list or numpy.ndarray, float) – The rake of nodal planes 1, angle between the fault’s horizontal and the slip direction of the hanging wall w.r.t. the foot wall (0-360 or -180-180).

• shear_update_atol (float, default to 1e-5) – Convergence criterion on the shear stress magnitude updates. Convergence is reached when the RMS difference between two estimates of shear stress magnitudes falls below that threshold.

• max_n_iterations (integer, default to 300) – The maximum number of iterations if shear stress magnitude update does not fall below shear_update_atol.

• Tarantola_kwargs (Dictionary, default to {}:) – If not None, should contain key word arguments for the Tarantola and Valette inversion. An empty dictionary uses the default values in Tarantola_Valette. If None, uses the Moore-Penrose inverse.

• return_eigen (boolean, default to True) – If True, returns the eigendecomposition of the inverted stress tensor in addition to returning the stress tensor.

• return_stats (boolean, default to True) – If True, the posterior data and model parameter distributions estimated from the Tarantola and Valette formula (cf. Tarantola_Valette routine).

Returns

output –

• output[“stress_tensor”]: (3, 3) numpy.ndarray

  The inverted stress tensor in the (north, west, up) coordinate system.

• output[“principal_stresses”]: (3,) numpy.ndarray, optional

  The three eigenvalues of the stress tensor, ordered from most compressive (sigma1) to least compressive (sigma3). Returned if return_eigen is True.

• output[“principal_directions”]: (3, 3) numpy.ndarray, optional

  The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i]. Returned if return_eigen is True.

• output[“C_m_posterior”]: (5, 5) numpy.ndarray, optional

  Posterior covariance of the model parameter distribution estimated from the Tarantola and Valette formula. Returned if return_stats is True.

• output[“C_d_posterior”]: (3 x n_earthquakes, 3 x n_earthquakes) numpy.ndarray, optional

  Posterior covariance of the data distribution estimated from the Tarantola and Valette formula. Returned if return_stats is True.

Return type

dict {str: numpy.ndarray}

ILSI.utils_stress.A_phi_(principal_stresses, principal_directions)

Compute A_phi as defined by Simpson 1997.

Parameters

• principal_stresses ((3,) numpy.ndarray) – The three eigenvalues of the stress tensor, ordered from most compressive (sigma1) to least compressive (sigma3).

• principal_directions ((3, 3) numpy.ndarray) – The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i].

Returns

A_phi

Return type

scalar, float

ILSI.utils_stress.R_(principal_stresses)

Computes the shape ratio R=(sig1-sig2)/(sig1-sig3).

Parameters

pinricpal_stresses (numpy.ndarray or list) – Contains the three eigenvalues of the stress tensor ordered such that: principal_stresses[0] < principal_stresses[1] < principal_stresses[2] with principal_stresses[0] being the most compressional stress.

Returns

shape_ratio

Return type

scalar float

ILSI.utils_stress.angular_residual(stress_tensor, strikes, dips, rakes)

Compute the angle between the direction of the resolved shear stress predicted by the stress tensor and the direction of slip given by the strike/dip/rake data.

Parameters

• stress_tensor ((3, 3) numpy.ndarray) – The Cauchy stress tensor.

• strikes ((n_earthquakes) list or numpy.ndarray) – Fault strikes.

• dips ((n_earthquakes) list or numpy.ndarray) – Fault dips.

• rakes ((n_earthquakes) list or numpy.ndarray) – Fault rakes.

Returns

angles – Angles between shear stress and slip.

Return type

(n_earthquakes) numpy.ndarray

ILSI.utils_stress.aux_plane(s1, d1, r1)

Get Strike and dip of second plane.

Adapted from MATLAB script bb.m written by Andy Michael, Chen Ji and Oliver Boyd.

Taken from <https://docs.obspy.org/_modules/obspy/imaging/beachball.html#aux_plane>. See Obspy project at <https://github.com/obspy/obspy>.

ILSI.utils_stress.check_right_handedness(basis)

Make sure the matrix of column vectors forms a right-handed basis. This is particularly important when re-ordering the principal stress directions based on their eigenvalues.

Parameters

basis ((3, 3) numpy.ndarray) – Matrix with column vectors that form the basis of interest.

Returns

rh_basis – Matrix with column vectors that form the right-handed version of the input basis. One of the unit vectors might have been reversed in the process.

Return type

(3, 3) numpy.ndarray

ILSI.utils_stress.compute_traction(stress_tensor, normal)

Compute traction and its normal and shear components on a given plane.

Parameters

• stress_tensor ((3, 3) numpy.ndarray) – Cauchy stress tensor.

• normal ((n_earthquakes, 3) numpy.ndarray) – Matrix of n_earthquakes row vectors of fault normals.

Returns

• traction ((n_earthquakes, 3) numpy.ndarray) – Tractions on the surfaces defined by normal.

• normal_traction ((n_earthquakes, 3) numpy.ndarray) – Normal component of the tractions.

• shear_traction ((n_earthquakes, 3) numpy.ndarray) – Tangential component of the tractions.

ILSI.utils_stress.errors_in_data(strike, dip, rake, jack_strikes_1, jack_dips_1, jack_rakes_1, jack_strikes_2, jack_dips_2, jack_rakes_2)

This routines was tailored for my applications. Use the multiple solutions obtained during the jackknife resampling of the focal mechanism inversion to compute the deviation of these multiple solutions from the best solution. A low deviation means a good quality focal mechanism. Because there are two possible slip vectors for each focal mechanism solution, we systematically look among the jackknife solutions 1 and 2 for the closest slip vector to the target vector, defined by (strike, dip, rake). We recommend to run this function for (strike, dip, rake)_1 and (strike, dip, rake)_2 of the best focal mechanism solution, and average the outputs.

ILSI.utils_stress.get_CI_levels(azimuths, plunges, confidence_intervals=[95.0, 90.0], nbins=200, smoothing_sig=1, return_count=False, plot=False)

Computes the 2d histogram in the stereographic space of a collection lines described by their azimuth and plunge.

Parameters

• azimuths ((n_lines) list or numpy.ndarray, float) – Azimuths of the lines.

• plunges ((n_lines) list or numpy.ndarray, float) – Plunges (angle from horizontal) of the lines.

• nbins (integer, default to 200) – Number of bins, in both axes, used to discretized the 2d space.

• smoothing_sig (float, default to 1) – If greater than 0, smooth the 2d distribution with a gaussian kernel. This is useful to derive smooth confidence intervals.

• plot (boolean, default to False) – If True, plot the 2d histogram.

Returns

• count ((nbins, nbins) numpy.ndarray, integer, optional) – 2D histogram of the lines dsecribed by azimuths and plunges. Only provided if return_count is True.

• lons_g ((nbins, nbins) numpy.ndarray, float, optional) – 2D grid of the longitudinal coordinate of each bin. Only provided if return_count is True.

• lats_g ((nbins, nbins) numpy.ndarray, float, optional) – 2D grid f the latitudinal coordinate of each bin. Only provided if return_count is True.

• confidence_intervals ((nbins, nbins) numpy.ndarray, float) – 2D distribution of the mass.

ILSI.utils_stress.get_CI_levels_joint(azimuths, plunges, confidence_intervals=[90.0, 95.0], nbins=200, smoothing_sig=1, return_count=False, plot=False)

Computes the 2d histogram in the stereographic space of a collection lines described by their azimuth and plunge. This is an EXPERIMENTAL function.

Parameters

• azimuths ((n_lines) list or numpy.ndarray, float) – Azimuths of the lines.

• plunges ((n_lines) list or numpy.ndarray, float) – Plunges (angle from horizontal) of the lines.

• nbins (integer, default to 200) – Number of bins, in both axes, used to discretized the 2d space.

• smoothing_sig (float, default to 1) – If greater than 0, smooth the 2d distribution with a gaussian kernel. This is useful to derive smooth confidence intervals.

• plot (boolean, default to False) – If True, plot the 2d histogram.

Returns

• count ((nbins, nbins) numpy.ndarray, integer, optional) – 2D histogram of the lines dsecribed by azimuths and plunges. Only provided if return_count is True.

• lons_g ((nbins, nbins) numpy.ndarray, float, optional) – 2D grid of the longitudinal coordinate of each bin. Only provided if return_count is True.

• lats_g ((nbins, nbins) numpy.ndarray, float, optional) – 2D grid f the latitudinal coordinate of each bin. Only provided if return_count is True.

• confidence_intervals ((nbins, nbins) numpy.ndarray, float) – 2D distribution of the mass.

ILSI.utils_stress.get_bearing_plunge(u, degrees=True, hemisphere='lower')

The vectors are in the coordinate system (x1, x2, x3): x1: north x2: west x3: upward

Parameters

• u – Vector for which we want the bearing (azimuth) and plunge.

• degrees (boolean, default to True) – If True, returns bearing and plunge in degrees. In radians otherwise.

• hemisphere (string, default to 'lower') – Consider the intersection of the line defined by u with the lower hemisphere if hemisphere is ‘lower’, or with the upper hemisphere if hemisphere is ‘upper’.

Returns

• bearing (float) – Angle between the north and the line.

• plunge (float) – Angle between the horizontal plane and the line.

ILSI.utils_stress.hist2d(azimuths, plunges, nbins=200, smoothing_sig=0, plot=False)

Computes the 2d histogram in the stereographic space of a collection lines described by their azimuth and plunge.

Parameters

• azimuths ((n_lines) list or numpy.ndarray, float) – Azimuths of the lines.

• plunges ((n_lines) list or numpy.ndarray, float) – Plunges (angle from horizontal) of the lines.

• nbins (integer, default to 200) – Number of bins, in both axes, used to discretized the 2d space.

• smoothing_sih (float, default to 0) – If greater than 0, smooth the 2d distribution with a gaussian kernel. This is useful to derive smooth confidence intervals.

• plot (boolean, default to False) – If True, plot the 2d histogram.

ILSI.utils_stress.joint_CDF(count)

ILSI.utils_stress.kagan_angle(tensor1, tensor2)

Compute the minimum rotation about some axis required to match the two tensors. This angle is a measure of their difference.

Parameters

• tensor1 ((3, 3) numpy.ndarray) – First tensor, e.g. moment or stress tensor.

• tensor2 ((3, 3) numpy.ndarray) – Second tensor, e.g. moment of stress tensor.

Returns

rotation_angle – Smallest angle, in degrees, required to superimpose the two tensors.

Return type

scalar float

ILSI.utils_stress.mean_angular_residual(stress_tensor, strikes, dips, rakes)

Mean of the absolute value of the angles returned by angular_residual. See angular_residual for more info.

ILSI.utils_stress.mean_kagan_angle(strikes, dips, rakes, strike0=None, dip0=None, rake0=None)

Computes the mean kagan angle as a measure of dispersion.

The mean kagan angle within a population of focal mechanisms described by strikes/dips/rakes. If strike0, dip0, and rake0 are specified, then the mean kagan angle is computed not from all pairs of focal mechanisms, but only between all focal mechanisms and the reference focal mechanism described by strike0/dip0/rake0. The mean kagan angle can be interpreted as a measure of dispersion within the population.

Parameters

• strikes (numpy.ndarray or list, float) – Strikes of the moment tensors.

• dips (numpy.ndarray or list, float) – Dips of the moment tensors.

• rakes (numpy.ndarray or list, float) – Rakes of the moment tensors.

• strike0 (scalar, float, default to None) – Strike of the reference moment tensor.

• dip0 (scalar, float, default to None) – Dip of the reference moment tensor.

• rake0 (scalar, float, default to None) – Rake of the reference moment tensor.

Returns

mean_angle – Mean kagan angle between the moment tensors given as input.

Return type

scalar, float

ILSI.utils_stress.normal_slip_vectors(strike, dip, rake, direction='inward')

Determine the normal and the slip vectors of the focal mechanism defined by (strike, dip, rake). From Stein and Wysession 2002.

N.B.: This is the normal of the FOOT WALL and the slip of the HANGING WALL w.r.t the foot wall. It means that the normal is an inward-pointing normal for the hanging wall, and an outward pointing-normal for the foot wall.

The vectors are in the coordinate system (x1, x2, x3): x1: north x2: west x3: upward

Parameters

• strike (float) – Strike of the fault.

• dip (float) – Dip of the fault.

• rake (float) – Rake of the fault.

• direction (string, default to 'inward') – If ‘inward’, returns the inward normal of the HANGING wall, which is the formula given in Stein and Wysession. Equivalently, this is the outward normal of the foot wall. If ‘outward’, returns the outward normal of the HANGING wall, or, equivalently, the inward normal of the hanging wall.

Returns

• n ((3) numpy.ndarray) – The fault normal.

• d ((3) numpy.ndarray) – The slip vector given as the direction of motion of the hanging wall w.r.t. the foot wall.

ILSI.utils_stress.p_t_b_axes(normal, slip)

Determine the P (most compressive), T (least compressive) and B (intermediate, or neutral axis) axes from the normal and the slip vectors, following Stein and Wysession 2002, Section 4.5.2. (P, T, B) forms an orthogonal basis.

The vectors are in the coordinate system (x1, x2, x3): x1: north x2: west x3: upward

ILSI.utils_stress.principal_faults(stress_tensor, friction_coefficient)

Compute the orientation of the most unstable fault planes given a stress tensor and a coefficient of friction. These faults are called the principal faults.

Parameters

• stress_tensor ((3, 3) numpy.ndarray) – Cauchy stress tensor.

• friction_coefficient (scalar float) – Coefficient of friction used for the Mohr-Coulomb failure criterion.

Returns

• n1 ((3, 1) numpy.ndarray) – Normal of the first principal faults.

• n2 ((3, 1) numpy.ndarray) – Normal of the second principal faults. The two faults form a pair of conjugate faults.

ILSI.utils_stress.quaternion(t, p, b)

Formula of quaternion of rotation matrix with t (least compressive), p (most compressive), b (neutral) components expressed in the (north, east, down) frame of reference. t, p, b can equivalently be the sigma_3, sigma_1, sigma_2 components. Make sure (t, p, b) form a right-handed basis. This routine was copied from the _tpb2q routine of the Pyrocko Python project (see at https://pyrocko.org/docs/current/_modules/pyrocko/moment_tensor.html#kagan_angle).

Parameters

• t ((3,) numpy.ndarray or list) –

• p ((3,) numpy.ndarray or list) –

• b ((3,) numpy.ndarray or list) –

Returns

quaternion – The quaternion that represents the rotation represented by the matrix (t, p, b), where t, p, b are column vectors.

Return type

(4,) numpy.ndarray

ILSI.utils_stress.random_rotation(max_angle=360.0)

Generate a random rotation matrix.

Generate a random rotation matrix by:

1. Generate a random unit vector in 3D.

2. Generate a random rotation angle between 0 and max_angle (degrees)

Parameters

max_angle (scalar float, default to 360) – Upper bound of the uniform distribution from which the rotation angle is randomly drawn.

Returns

R – Rotation matrix.

Return type

(3, 3) numpy.ndarray

ILSI.utils_stress.reduced_stress_tensor(principal_directions, R)

Computes a normalized stress tensor where the most and least compressive principal stresses are set to -1 and +1, respectively, and the intermediate stress is determined by the shape ratio.

Parameters

• principal_directions ((3, 3) numpy.ndarray.) – The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i].

• R (float) – The shape ratio (sig1 - sig2)/(sig1 - sig3).

Returns

stress_tensor – The stress tensor built from the principal directions and the shape ratio.

Return type

(3, 3) numpy.ndarray

ILSI.utils_stress.rotation(axis, angle)

Compute the rotation matrix about axis with angle angle.

Parameters

• axis – Axis about which the rotation is computed.

• angle (scalar, float) – Angle, in degrees, of the rotation.

Returns

R – Rotation matrix of angle angle degrees about axis.

Return type

(3, 3) numpy.ndarray, float

ILSI.utils_stress.round_cos(x)

Clip x so that it fits with the [-1,1] interval.

If x is slightly outside the [-1,1] because of numerical imprecision, x is rounded, and can then be safely passed to arccos or arcsin. If x is truly outside of [-1,1], x is returned unchanged.

Parameters

x (scalar, float) – Float variable that represents a cos or sin that is supposed to be within the [-1,1] interval.

Returns

x_r – A rounded version of x, if necessary.

Return type

scalar, float

ILSI.utils_stress.shear_slip_angle_difference(stress_tensor, strike, dip, rake)

Return the angle difference between the slip vector from the focal mechanism solution and the shear traction on the fault determined from the inverted stress tensor. Given that the stress inversion is made under the Wallace-Bott assumption, shear stress on the fault is parallel to slip, then this angle difference is a measure of misfit.

Parameters

• stress_tensor ((3, 3) numpy.ndarray) – The Cauchy stress tensor.

• strike (float) – Strike of the fault.

• dip (float) – Dip of the fault.

• rake (float) – Rake of the fault.

Returns

angle – The angle between shear stress and slip, in degrees.

Return type

float

ILSI.utils_stress.stress_tensor_eigendecomposition(stress_tensor)

Compute the eigendecomposition of stress tensor.

Parameters

stress_tensor ((3, 3) numpy.ndarray.) – The stress tensor for which to solve the eigenvalue problem.

Returns

• principal_stresses ((3,) numpy.ndarray.) – The three eigenvalues of the stress tensor, ordered from most compressive (sigma1) to least compressive (sigma3).

• principal_directions ((3, 3) numpy.ndarray.) – The three eigenvectors of the stress tensor, stored in a matrix as column vectors and ordered from most compressive (sigma1) to least compressive (sigma3). The direction of sigma_i is given by: principal_directions[:, i].

ILSI.utils_stress.strike_dip_rake(n, d)

Invert the relationships between strike/dip/rake and normal (n) and slip (d) vectors found in Stein. n and d are required to be given as the default format returned by normal_slip_vectors.

Parameters

• n – The outward pointing normal of the FOOT wall.

• d – The slip direction of the hanging wall w.r.t. the foot wall.

Returns

• strike (float) – Strike of the fault, in degress.

• dip (float) – Dip of the fault, in degrees.

• rake (float) – Rake of the fault, in degrees.

ILSI.utils_stress.strike_dip_rake_to_mt(strike, dip, rake)

Compute the normalized moment tensor described by strike/dip/rake.

Parameters

• strike (scalar, float) – Strike of the input focal mechanism.

• dip (scalar, float) – Dip of the input focal mechanism.

• rake (scalar, float) – Rake of the input focal mechanism.

Returns

mt – Normalized moment tensor. Its columns are the (p, b, t) axes. Note: we return (p, b, t) to be consistent with our stress tensor convention (sig1, sig2, sig3).

Return type

(3, 3) numpy.ndarray, float