About Bstrain

This code has been developed based on previous work done by Thomas Bodin and coauthors for tomography. Find the source code https://forge.univ-lyon1.fr/marianne.metois/bayesianstrainrate/

It has been modified to apply the transdimensional bayesian inversion framework to the problem of deriving strain rates from discrete 2D GNSS velocity fields or displacements.

Its primary version has been developed at LGL-TPE by Colin Pagani during his PhD (obtained in 2021) and applied to California as a proof of concept (see Pagani et al. 2021).

Warning

If you use this code, please cite :

Pagani, C., Bodin, T., Métois, M., & Lasserre, C. (2021). Bayesian estimation of surface strain rates from Global Navigation Satellite System Measurements: Application to the Southwestern United States. Journal of Geophysical Research: Solid Earth, 126(6), e2021JB021905.https://doi.org/10.1029/2021JB021905

Tip

If you enconter problems, please post an issue in the gitlab, or contact marianne.metois-at-univ-lyon1.fr

The code is still under developpement. In particular, we plan to adapt the code for 3D GNSS and to adapt it to InSAR data sets (LOS directed deformation). If you are interested in collaborating, get in touch, our team (Colin, Thomas, Cécile and Marianne) will be pleased to discuss.

Compute infinitesimal vs Lagrangian strain rate tensor

Two different formalism exist to compute the strain rate or strain tensor from velocity (resp. displacement) fields. Bstrain has been firstly designed to deal with small deformation, and therefore to compute the infinitesimal strain rate tensor \(\dot\epsilon_{ij}\). We set \(\mathbf{\nabla v}\) that can be decomposed into a symmetric strain rate tensor \( \mathbf{\dot\epsilon}\) and an antisymmetric rigid body rotation matrix \(W\) so that :

\[\dot\epsilon_{ij} = \frac{1}{2}(\partial_jv_i + \partial_iv_j), \text{ and}\]

\[\begin{split}W = \begin{pmatrix} 0 & \omega \\ -\omega & 0 \\ \end{pmatrix} \text{ with } \omega=\frac{1}{2}(\partial_yv_x - \partial_xv_y).\end{split}\]

However, when dealing with large deformations such as coseismic motions, or creeping plate boundaries, one may prefer to compute the Lagrangian version of the strain rate tensor. Bstrain has an option to compute this tensor rather than the infinitesimal one if needed.

Let’s consider two points \(\overrightarrow{X} \) and \( \overrightarrow{x}\) as described in the following figure, moving by \(\overrightarrow{dX}\) and \(\overrightarrow{dx}\) respectively.

We can write:

\[\overrightarrow{x} + \overrightarrow{dx} = \overrightarrow{X} + \overrightarrow{dX} + \overrightarrow{u}(\overrightarrow{X} + \overrightarrow{dX})\]

\[\overrightarrow{dx} = \overrightarrow{X} - \overrightarrow{x} + \overrightarrow{dX} + \overrightarrow{u}(\overrightarrow{X} + \overrightarrow{dX})\]

\[\overrightarrow{dx} = \overrightarrow{dX} + \underbrace{\overrightarrow{u}(\overrightarrow{X} + \overrightarrow{dX}) - \overrightarrow{u}(\overrightarrow{X})}_{\frac{\partial{\overrightarrow{u}}}{\partial{\overrightarrow{X}}}\cdot\overrightarrow{dX}}\]

Now,

\[\frac{\partial{\overrightarrow{u}}}{\partial{\overrightarrow{X}}}\cdot\overrightarrow{dX} = \nabla\overrightarrow{u}\cdot\overrightarrow{dX}\]

hence

\[\overrightarrow{dx} = (\mathbf{I}+\nabla\overrightarrow{u})\cdot\overrightarrow{dX}\]

By definition, we denote \(F\) as the deformation matrix given by \(F=(\mathbf{I}+\nabla\overrightarrow{u})\), which can also be written as \(F_{ij}=\frac{\partial{x_i}}{\partial{X_j}}\).

The Cauchy tensor \(C\) is given by

\[C=F^T\cdot F\]

\[C=(\mathbf{I}+\nabla\overrightarrow{u})^T\cdot(\mathbf{I}+\nabla\overrightarrow{u})\]

\[C= \mathbf{I} + \underbrace{\nabla\overrightarrow{u} + \nabla\overrightarrow{u}^T + \nabla\overrightarrow{u}^T\cdot\nabla\overrightarrow{u}}\]

\[C=\mathbf{I} + 2E\]

where \(E\) is by definition the Lagrangian strain tensor given by:

\[E=0.5(\nabla\overrightarrow{u} + \nabla\overrightarrow{u}^T + \nabla\overrightarrow{u}^T\cdot\nabla\overrightarrow{u})\]

that is, in index notation

\[E_{ij} = \frac{1}{2}(\partial_j u_i + \partial_i u_j + \underbrace{\partial_i u_j\cdot\partial_j u_i}_{A})\]

\(A\) is generally considered negligible if \(\nabla\overrightarrow{u}<<0\), i.e., for small deformations, and we then have \(E_{ij}=\epsilon_{ij}=\frac{1}{2}(\partial_j u_i + \partial_i u_j)\). This is not the case for a coseismic displacement, for example, but potentially also not for a very fast plate boundary like the San Andreas.

Specifically, in the 2D tensor \(\mathbf{E}\), we thus need to add:

\[E_{xx}=\partial_x u_x + \boxed{\frac{1}{2}(\partial_xu_x)^2}\]

and

\[E_{xy}=\frac{1}{2}(\partial_x u_y + \partial_y u_x + \boxed{\partial_x u_y\cdot\partial_y u_x})\]