Tatic stress change from a kinematic fault is ~10–100× larger than expected — what could cause this?

General
1

Hello,

I’m using PyLith 2.2.2 (MPI, tet4 mesh from Cubit) to compute the static stress change (stress drop) on a subduction megathrust, by imposing an interseismic slip deficit as kinematic slip on the fault with FaultCohesiveKin. The computed fault traction change is roughly 10–100× larger than expected, and I’ve been unable to find the cause. I’d appreciate any pointers.

Model setup

  • 3D, unstructured tet4 mesh (Cubit/Exodus).
  • Domain ≈ 170 × 180 × 80 km (x = strike-normal, y = along-strike, z = depth).
  • Fault: smoothly curved megathrust, ≈ 130 km down-dip × 180 km along-strike, dip ~11°→48°.
  • Material: linear elastic; Vs ≈ 3.0–4.5 km/s, ρ ≈ 2800 kg/m³ → μ ≈ 25–40 GPa (I also tested a homogeneous half-space — same result).
  • Boundary conditions: Dirichlet on all faces except the free surface (top). I have tried both “roller” (one component per face) and fully fixed (all 3 components) — essentially the same result.
  • Fault: FaultCohesiveKin, label = fault, with and without edge = fault_edge (no difference).
  • Imposed slip: pure reverse (dip) slip, no strike-slip/opening. I tested both the spatially variable slip deficit (max ≈ 5.08 m) and a UNIFORM 1 m slip (to isolate the effect).

Observation (uniform 1 m slip)

  • Shear traction change, dip component (≈ −traction_change[:,1]): ≈ 29 MPa at fault center; 20–47 MPa over the interior.
  • Simple estimate Δτ ~ μ·s/L ~ (30 GPa × 1 m)/(130 km) ~ 0.2 MPa → observed is ~100× larger.
  • For the real slip deficit (max 5.08 m), Δτ reaches ~200 MPa (expected ~10 MPa).
  • Domain total_strain reaches ~1% (expected ~0.004%).

What I’ve already checked (no improvement)

  1. Material properties (μ ≈ 25–40 GPa — looks correct).
  2. Slip is genuinely pure dip-slip (slip components = [0, s, 0]).
  3. Boundary conditions: roller vs. fully fixed — no change.
  4. Fault edge nodeset: with / without — no change.
  5. Extending the fault sheet beyond the domain boundary — no change (or slightly worse).
  6. up_dir / fault orientation.

One thing I suspect: the fault currently spans the full along-strike extent of the domain (fault y = 0–180 km = domain y = 0–180 km) and its down-dip end reaches the domain face, so the fault’s edges lie directly on the fixed (Dirichlet) boundaries. Could this “pinning” of the fault edges to the fixed boundaries produce such amplification? What is the recommended way to avoid it (e.g., making the domain much larger than the fault so the fault is fully embedded, or excluding the fault edges from the BC nodesets)?

Question: What common mistakes could cause the static stress change to be ~10–100× too large? Any guidance on geometry, up_dir, boundary conditions, or units would be very helpful.

Thank you.

image
image2023×1093 344 KB

2

For the geometry you describe, the fault slip is definitely interacting with the boundaries of the domain.

If you are trying to simulate an infinitely-long fault along strike, then this may be an appropriate model, but you should not impose Dirichlet BC on the vertices that are on the fault.

If you are trying to simulate a portion of a subduction zone, then the fault should be in the middle of the domain and you should extend the domain several fault lengths away from the edges of the fault.

If the traction changes are spatially heterogeneous, then the fault may not be smooth relative to the discretization size. It is difficult to tell from your figure, what the discretization size is over the fault surface, although it looks quite small at the edges. If the discretization size is as fine in the middle of the fault as it is at the edges, it looks like the discretization size may be too small rather than too coarse, which could be leading to problems with the nondimensionalization.

You may want to start with a planar fault and a rather coarse discretization to test the basic model setup (material properties, boundary conditions, fault slip, etc). For example, with a similar model setup to what you have a planar fault with prescribed slip (without slip pinned at the edges of the domain) should produce rigid body motion and zero change in tractions.