[cig-commits] commit: Added draft of Savage & Prescott benchmark section.

[Mercurial]

changeset:   115:31302b5f6e5c
parent:      113:d22fcebab00c
user:        Charles Williams <C.Williams at gns.cri.nz>
date:        Tue May 01 16:42:18 2012 +1200
files:       faultRup.tex
description:
Added draft of Savage & Prescott benchmark section.


diff -r d22fcebab00c -r 31302b5f6e5c faultRup.tex
--- a/faultRup.tex	Thu Apr 26 11:42:57 2012 -0700
+++ b/faultRup.tex	Tue May 01 16:42:18 2012 +1200
@@ -1292,6 +1292,106 @@ benchmarks.
 \subsection{Quasi-static}
 \label{sec:verification:quasi-static}
 
+As a test of our quasi-static solution, we compare our numerical results
+against the analytical solution of \citep{Savage:Prescott:1978}. This
+problem consists of a strike-slip fault in an elastic layer 
+overlying a Maxwell viscoelastic half-space. The geometry of the problem
+is shown in Figure~\ref{fig:savage:prescott::solution}, which shows an
+exaggerated view of the deformation during the tenth earthquake in a
+sequence. Between earthquakes, the upper portion of the fault is locked,
+while the lower portion slips at plate velocity. At regular intervals
+(earthquake recurrence time), the upper portion of the fault ruptures
+such that the total slip on the locked portion of the fault matches the
+slip on the creeping portion.
+
+This problem provides a good test of the kinematic fault implementation
+in PyLith, including the ability to specify complex, time-varying
+boundary conditions on the fault. It also provides a good test of the
+Maxwell viscoelastic bulk rheology. The analytical solution for this
+problem provides surface displacements as a function of distance from
+the fault and time since the last earthquake. An infinite strike-slip
+fault is assumed, so there is no geometry along the strike of the
+fault. The solution is controlled by the ratio of the fault locking
+depth ($D$) to the thickness of the elastic layer ($H$), and by the
+ratio of the earthquake recurrence time ($T$) to the viscoelastic
+relaxation time:  $tau_0 = \mu T/2\eta$, where $\mu$ is the shear
+modulus and $\eta$ is the viscosity.
+
+For our model comparison, we used a locking depth of 20 km, an elastic
+layer thickness of 40 km, an earthquake recurrence time of 200 years, a
+shear modulus of 30 GPa, and a viscosity of $2.37\times 10^{19}$
+Pa-s. We used a relative plate velocity of 2 cm/year, implying a
+coseismic offset of 4 m every 200 years. The viscosity and shear modulus
+values yield a viscoelastic relaxation time of 50 years, and a $tau_0$
+value of 4.
+
+The problem is simulated in PyLith using a 3D mesh with (x, y, z)
+dimensions of (2000 km, 1000 km, 400 km). Velocity boundary conditions
+in the y-direction were applied to the -x and +x faces, and
+x-displacements were also constrained to zero on these faces. The
+z-displacements were constrained to zero on the -z face. Finally, we
+constrained the x-displacements to zero on the -y and +y faces to more
+accurately represent an infinite fault. For comparison with analytical
+results, we extracted the numerical results along an x-profile at y = 0,
+corresponding to the center of the mesh in the y-direction. We solved
+the problem using both trilinear hexahedral cells as well as linear
+tetrahedral cells, using two different resolutions for each cell
+type. In our coarsest hexahedral mesh we used a uniform resolution of
+20 km. In our higher resolution hexahedral mesh we refined an inner
+region (x-dimension = 480 km, y-dimension = 240 km, z-dimension = 100
+km) by a factor of 3, yielding a resolution near the center of the fault
+of 6.7 km. For the tetrahedral meshes, we maintained the same resolution
+near the center of the fault (20 km and 6.7 km); however, for these
+meshes we linearly increased the cell size to the outer edges of the
+mesh. This resulted in cells with maximum dimensions of approximately 60
+km for the coarser mesh and 40 km for the higher resolution mesh. Note
+that for both the hexahedral and tetrahedral coarse meshes, the cell
+size on the fault is the maximum allowable size that still allows us to
+represent the fault locking depth as a sharp boundary.
+
+Since this is a viscoelastic problem, it is necessary to 'spin-up' the
+solution for several earthquake cycles until a near steady-state is
+obtained. There is an additional issue for the numerical solution. In
+the analytical solution, steady plate motion is simply
+superimposed. This is not possible with the numerical solution; however,
+after several earthquake cycles we approach this state. We run both the
+analytical and numerical simulations for 10 earthquake cycles, for a
+total simulation time of 2000 years. For the PyLith solution we use a
+constant time step size of 5 years. Note that this time step size is
+1/10 of the viscoelastic relaxation time, and should therefore provide a
+rigorous test of the accuracy of the viscoelastic solution for
+moderately large time steps. The fault slip is specified using two
+different kinematic slip functions: the default StepSlipFn, which allows
+us to specify slip for each rupture event, and the ConstRateSlipFn,
+which allows us to specify steady slip on the creeping portion of the
+fault.
+
+A comparison of the analytical and numerical results is shown in Figure
+~\ref{fig:savage:prescott:profiles}. To examine the differences very
+close to the fault, we have used a logarithmic scale for the x
+axis. Further from the fault, all models show virtually identical
+results. The top portion of the figure shows the results during the
+third earthquake cycle, while the bottom portion shows the results
+during the tenth earthquake cycle. Close to the fault, all the
+simulations show similar behavior between the third and tenth
+cycle. Further from the fault, although all of the numerical solutions
+predict identical displacements, the viscoelastic solution has not yet
+achieved steady state, and thus under-predicts the displacement.
+
+By the tenth earthquake cycle, all of the solutions predict nearly
+identical far-field displacements. The differences all occur within one
+elastic layer thickness of the fault. The two higher resolution
+numerical models predict nearly identical results, and provide a very
+close fit to the analytical solution. Close to the fault, the effects of
+inadequate discretization become apparent for the coarse meshes, as the
+details of the complex fault slip cannot be accurately represented.
+
+Benchmarks such as this can be quite helpful when designing meshes for
+real-world problems. As seen in this benchmark, a coarse discretization
+can still yield accurate simulations as long as only far-field results
+are desired. If the region of interest is close to the fault, a finer
+mesh will be needed.
+
 \begin{itemize}
 \item Savage and Prescott
   \begin{itemize}
@@ -1527,7 +1627,7 @@ MGK acknowledges partial support from NS
   \noindent\includegraphics[width=84mm]{figs/savageprescott_soln}
   \caption{Deformation (exaggerated by a factor of 5000) 95\% of the
     way through earthquake cycle 10 of the Savage and Prescott
-    benchmark involving viscoelastic relaxation over multiple
+    benchmark, which involves viscoelastic relaxation over multiple
     earthquake cycles on a vertical, strike-slip fault. The
     coordinates are in units of locking depth and the displacement
     field is in units of coseismic slip. The locking depth is one-half
@@ -1546,7 +1646,7 @@ MGK acknowledges partial support from NS
     cycles 3 and 10. The displacements values shown are
     relative to the values at the beginning of the earthquake cycle to
   facilitate comparison between the analytical solution and the
-  numerical models which require spin-up to reach the steady-state
+  numerical models, which require spin-up to reach the steady-state
   solution. Both the hexahedral (Hex8) and tetrahedral (Tet4)
   discretizations resolve the viscoelastic deformation and display
   excellent agreement with the steady-state solution by the tenth