Slip on a Fault

In this example we illustrate how to calculate the stress distribution around a fault in the Earth's crust caused by a slip through an earthquake.

To simplify the presentation we assume a simple domain $\Omega=[0,1]^2$ with a vertical fault in its center. We assume that the slip distribution $s\hackscore{i}$ on the fault is known. We want to calculate the distribution of the displacements $u\hackscore{i}$ and stress $\sigma_{ij}$ in the domain. We assume an isotropic, linear elastic material model of the form

$$\sigma\hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})$$ (64)

where $\lambda$ and $\mu$ are the Lame coefficients and $\delta\hackscore{ij}$ denotes the Kronecker symbol. On the boundary the normal stress is given by

$$\sigma\hackscore{ij}n\hackscore{j}=0$$ (65)

and normal displacements are set to zero:

$$u\hackscore{i}n\hackscore{i} =0$$ (66)

The stress needs to fulfill the momentum equation

$$- \sigma\hackscore{ij,j}=0$$ (67)

This problem is very similar to the elastic deformation problem presented in Section .

However, we need to address an additional challenge: The displacement $u\hackscore{i}$ is in fact discontinuous across the fault but we are in the lucky situation that we know the jump of the displacements across the fault. This is in fact the the given slip $s\hackscore{i}$. So we can split the total distribution $u\hackscore{i}$ into a component $v\hackscore{i}$ which is continuous across the fault and the known slip $s\hackscore{i}$

$$u\hackscore{i} = v\hackscore{i} + \frac{1}{2} s^{\pm}\hackscore{i}$$ (68)

where $s^{\pm}=s$ when right of the fault and $s^{\pm}=-s$ when left of the fault. We assume that $s^{\pm}=0$ sufficiently away from the fault.

We insert this into the stress definition 

$$\sigma\hackscore{ij} = \sigma^c\hackscore{ij} + \frac{1}{2} \sigma^s\hackscore{ij}$$ (69)

with

$$\sigma^c\hackscore{ij} = \lambda v\hackscore{k,k} \delta\hackscore{ij} + \mu ( v\hackscore{i,j} + v\hackscore{j,i})$$ (70)

and

$$\sigma^s\hackscore{ij} = \lambda s^{\pm}\hackscore{k,k} \delta\hackscore{ij} + \mu ( s^{\pm}\hackscore{i,j} + s^{\pm}\hackscore{j,i})$$ (71)

In fact $\sigma^s\hackscore{ij}$ defines a stress jump across the fault. In easy way to construct this function is to us a function $\chi$ which is $1$ on the right and $-1$ on the left side from the fault. One can then set

$$\sigma^s\hackscore{ij} = \chi \cdot ( \lambda s\hackscore{k,k} \delta\hackscore{ij} + \mu ( s\hackscore{i,j} + s\hackscore{j,i}) )$$ (72)

assuming that $s$ is extended by zero away from the fault. After inserting  into  we get the differential equation

$$- \sigma^c\hackscore{ij,j}=\frac{1}{2} \sigma^s\hackscore{ij,j}$$ (73)

Together with the definition  we have a differential equation for the continuous function $v_i$. Notice that the boundary condition  and  transfer to $v_i$ and $\sigma^c\hackscore{ij}$ as $s$ is zero away from the fault. In Section  we have discussed how this problem is solved using the LinearPDE class. We refer to this section for further details.

To define the fault we use the FaultSystem class introduced in Section . The following statements define a fault system fs and add the fault 1 to the system.

The fault added starts at point $(0.5,0.25)$ has length $0.5$ and is pointing north. The main purpose of the FaultSystem class is to define a parameterization of the fault using a local coordinate system. One can inquire the class to get the range used to parameterize a fault.

Typically p0 is equal to zero while p1 is equal to the length of the fault. The parameterization is given as a mapping from a set of local coordinates onto parameter range (in our case the range p0 to p1). For instance to map the entire domain mydomain onto the fault one can use

Of course there is the problem that not all location are on the fault. For those locations which are on the fault m is set to $1$ otherwise 0 is used. So on return the values of p are defining the value of the fault parameterization (typically the distance from the starting point of the fault along the fault) where m is positive. On all other locations the value if p is undefined. Now p can be used to define a slip distribution on the fault via

Notice that the factor m which makes sure that s is zero away from the fault. It is important that the slip is zero at the ends of the faults.

Figure: Total Displacement after slip event.

We can now put all components together to get the script:

The script creates the slip.vtu giving the total displacements and stress. These values are stored as cell centered data. See Figure  for the result.