Benchmark Description 

Summary

This benchmark problem computes the viscoelastic (Maxwell) relaxation of stresses from repeated infinite, strike-slip earthquakes in 3D without gravity. The analytical solution for this problem is presented by Savage and Prescott (JGR, 1978). Although the analytical solution is for an elastic layer overlying a viscoelastic half-space assuming an infinitely long strike-slip fault, the problem can be very well approximated using a 3D finite element model. The problem is driven by a combination of constant velocities applied along the x-boundaries of the mesh, along with periodic (coseismic) fault slip applied along the upper (locked) portion of the fault and steady fault slip applied in the portion of the elastic layer below the locked section of the fault. Since the finite element solution is an approximation of a half-space solution, the mesh dimensions are somewhat arbitrary. The dimensions listed below, however, have been shown to provide a good fit to the analytical solution.

Problem Specification

Model size

-1000 km ≤ x ≤ 1000 km; -500 km ≤ y ≤ 500 km; -400 km ≤ z ≤ 0 km

Top layer (elastic)

-40 km ≤ z ≤ 0 km

Bottom layer (viscoelastic)

-400 km ≤ z ≤ -40 km

Material properties

The top layer is elastic whereas the bottom layer is viscoelastic.

Elastic

Poisson solid, G = 30 GPa

Viscoelastic

Poisson solid, G = 30 GPa, η = 4.73364e19 Pa-s (yielding a relaxation time of 2η/G = 100 years)

Fault specifications

Type

Vertical left-lateral strike-slip fault.

Location

Strike parallel to y-direction at center of model (x = 0 km), extending along the entire y-dimension of the model (-500 km ≤ y ≤ 500 km), with a depth extent of -40 km ≤ z ≤ 0 km.

Slip distribution

The locked portion of the fault (red section in the Figure) extends from -20 km ≤ z ≤ 0 km, while the creeping section (blue) extends from -40 km ≤ z ≤ -20 km. Along the line where the two sections coincide (z = -20 km), half of the coseismic displacement and half of the steady creep is applied. On the locked portion of the fault, 4 m of left-lateral slip is applied every 200 years, representing a complete release of the accumulated slip deficit over that time. On the creeping portion of the fault, steady left-lateral slip of 2 cm/year is applied for the entire simulation.

Boundary conditions

Along the bottom of the mesh, the z-displacements are set to zero. On the x = -1000 km plane and the x = 1000 km plane, the y-velocities are set to -1 cm/year and +1 cm/year, respectively, giving a relative plate velocity of 2 cm/year. The x-displacements are set to zero on these planes, and they are also set to zero on the y = -500 km plane and the y = 500 km plane.

Discretization

The model should first be tried using a uniform resolution of 20 km (coarse resolution). A variable resolution mesh should also be tried. Options that have been used so far include a mesh with a nominal resolution of 20 km, with an inner refined region with extents of -240 km ≤ x ≤ 240 km; -120 km ≤ y ≤ 120 km; -100 km ≤ z ≤ 0 km. When using Cubit, this yields a resolution near the fault of 6.67 km for hexahedral meshes and 10 km for tetrahedral meshes.

Element types

Linear and/or quadratic and tetrahedral and/or hexahedral.

Requested Output

Solution on z = 0 km

Displacements at all nodes on the upper (free) surface at times of 10, 50, 100, 150, and 190 years as well as the mesh topology (i.e., element connectivity arrays and coordinates of vertices) and basis functions. Note that it is necessary to run the problem through multiple earthquake cycles to obtain a solution approximating the analytical solution. This is due to a difference in how the two models are spun up. Therefore, run the model through 9 earthquake cycles and then output the solution for the tenth cycle.

May 19, 2010

Use ASCII output for now. In the future we will switch to using HDF5 files.

Performance

• CPU time
• Wallclock time
• Memory usage

"Truth"

The analytical solution for displacements on the free surface has been given by Savage and Prescott (JGR, 1978). A Python utility to compute this solution is available in the geodynamics.org subversion repository. Also contained in the repository are PyLith input files, along with additional utilities and figures comparing analytical results with results using PyLith. The files may be obtained by entering the following SubVersion command:

svn checkout http://geodynamics.org/svn/cig/short/3D/PyLith/benchmarks/trunk/quasistatic/sceccrustdeform/savageprescott