[cig-commits] commit: Added conclusions. Edits to Savage and Prescott section.
Mercurial
hg at geodynamics.org
Fri May 4 11:23:45 PDT 2012
changeset: 119:c3fd6d615a85
user: Brad Aagaard <baagaard at usgs.gov>
date: Fri May 04 11:15:30 2012 -0700
files: faultRup.tex figs/savageprescott_profiles.pdf
description:
Added conclusions. Edits to Savage and Prescott section.
diff -r 77311386bf85 -r c3fd6d615a85 faultRup.tex
--- a/faultRup.tex Fri May 04 08:23:02 2012 -0700
+++ b/faultRup.tex Fri May 04 11:15:30 2012 -0700
@@ -1276,7 +1276,7 @@ is modified or optimized. We also rely o
is modified or optimized. We also rely on full-scale benchmarks to
verify that the code properly solves the numerical problem. These
benchmark problems include quasi-static strike-slip and reverse
-viscoelastic simulations and many of problems in the suite of
+viscoelastic simulations and various problems in the suite of
spontaneous dynamic rupture benchmarks developed by the Southern
California Earthquake Center (SCEC) and the United States Geological
Survey \citep{Harris:etal:SRL:2009}. In this section we focus on two
@@ -1292,129 +1292,98 @@ 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.
+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 infinitely
+long strike-slip fault in an elastic layer overlying a Maxwell
+viscoelastic half-space. Figure~\ref{fig:savage:prescott::solution}
+illustrates the geometry of the problem with an exaggerated view of
+the deformation during the tenth earthquake cycle. Between earthquakes
+the upper portion of the fault is locked, while the lower portion
+slips at the plate velocity. At regular intervals (the earthquake
+recurrence time) the upper portion of the fault ruptures such that the
+slip on the locked portion exactly complements the slip on the
+creeping portion so the cumulative slip over an earthquake cycle is
+uniform.
-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 in the y-direction (along-strike)
-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.
+This problem tests the ability the kinematic fault implementation to
+include steady aseismic creep and multiple earthquake ruptures complex
+along with viscoelastic relaxation. The analytical solution for this
+problem provides the along-strike component of surface displacement as
+a function of distance perpendicular to the fault. The solution is
+controlled by the ratio of the fault locking depth to the thickness of
+the elastic layer and the ratio of the earthquake recurrence time to
+the viscoelastic relaxation time, $\tau_0 = \mu T/2\eta$, where $T$ is
+the recurrence time, $\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.
+For this benchmark we use 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,
+and 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 $\tau_0 =4$.
+We employ a 3-D model (2000 km by 1000 km by 400 km) with Dirichlet
+boundary conditions enforcing symmetry to approximate an infinitely
+long strike-slip fault. We apply velocity boundary conditions in the
+y-direction to the -x and +x faces with zero x-displacement. We
+constrain the vertical displacements on the bottom of the
+domain. Finally, we fix the x-displacements on the -y and +y faces to
+enforce symmetry consistent with an infinitely long strike-slip fault.
-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
-represent more accurately 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
+We examine four difference numerical solutions considering the effects
+of cell type (hexahedral versus tetrahedral) and discretization
+size. In our coarse hexahedral mesh we use a uniform resolution of 20
+km. In our higher resolution hexahedral mesh we refine an inner region
+(480 km by 240 km by 100 km) by a factor of three, yielding a
+resolution near the center of the fault of 6.7 km. For the
+tetrahedral meshes, we match the discretization size of the hexehedral
+mesh near the center of the fault (20km or 6.7 km) while increasing
+the discretization size in a geometric progression at a rate of
+1.02. This results in a maximum discretization size 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 discretization 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
-achieved. 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 coseismic slip for each rupture event, and the
-ConstRateSlipFn, which allows us to specify steady slip on the creeping
-portion of the fault.
+In this viscoelastic problem the numerical model does not achieve
+steady-state behavior until after several earthquake cycles. The
+duration of this model spin-up depends on how close the initial
+conditions of the match the steady-state behavior and the relaxation
+time of the viscoelastic material. Shorter relaxation will converge to
+the steady-state behavior in fewer earthquake cycles. We simulate ten
+earthquake cycles for each of the numerical models for a total
+duration of 2000 years with a time step of five years. This time step
+corresponds to one tenth of the viscoelastic relaxation time; hence it
+tests the accuracy of the viscoelastic solution for moderately large
+time steps relative to the relaxation time.
-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
-second 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 second 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.
+Figure~\ref{fig:savage:prescott:profiles} compares the numerical
+results extracted on the ground surface along the center of the model
+perpendicular to the fault with the analytic solution. Using a
+logarithmic scale with distance from the fault facilitates examing the
+solution both close to the fault and in the far-field. In the
+far-field the simulations show virtually identical results; however,
+the viscoelastic solution has not yet achieved steady state and the
+simulations underpredict the displacement. By the tenth earthquake
+cycle, the simulations reach a steady state and match the analytical
+solution.
-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. Thus, with propert discretization size, PyLith is able to
-represent complex spatial and temporal kinematic slip distributions for
-quasi-static problems extremely accurately.
-
-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}
- \item Parameters (mesh, geometry, discretization size, physical properties)
- \item Figures
- \begin{itemize}
- \item Geometry
- \item profiles for cycles 3 and 10
- \end{itemize}
- \item Table of parameters
- \item Spin-up, compare hex8 and tet4 against analytic solution
- \end{itemize}
-\end{itemize}
-\brad{QUESTIONS FOR CHARLES: The hex8 versus tet4 comparison isn't very
- useful as both give essentially identical results. Would looking at
- another quantity (displacement along a vertical profile, or stress)
- provide a more stringent test? Is there an easy way to adjust the
- parameters to give something more numerically challenging
- (higher/lower viscosity, coarser mesh)?}
+Within about one elastic thickness of the fault the effect of the
+resolution of the numerical models becomes apparent. We find large
+errors for the coarse models, which have discretization sizes matching
+the elastic thickness. The finer resolution models (6.7 km
+discretization size) provide a close fit to the analytical
+solution. The 6.7 km hexahedral solution is indistinguishable from the
+analytical solution in Figure~\ref{fig:savage:prescott:profiles}(b);
+the 6.7 km tetrahedral solution slightly underpredicts the analytical
+solution. The greater accuracy of the hexahedral cells relative to
+the tetrahedral cells with the same nominal discretization size for
+quasi-static solutions is consistent with our findings for other
+benchmarks. The greater number of polynomial terms in the basis
+functions of the hexahedral allows the model to capture a more complex
+deformation field.
\subsection{Dynamic}
\label{sec:verification:dynamic}
@@ -1530,6 +1499,31 @@ supershear rupture and properly implemen
supershear rupture and properly implements a Drucker-Prager
elastoplastic bulk rheology and slip-weakening friction.
+% ------------------------------------------------------------------
+\section{Conclusions}
+\label{sec:conclusions}
+
+The domain decomposition approach to implementing fault slip in PyLith
+provides a flexible numerical implementation as evident by its
+application in quasi-static and dynamic simulations and kinematic
+(prescribed) slip and spontaneous rupture (frictional sliding) models
+of earthquake rupture. We find that algebraic multigrid
+preconditioners for elasticity combined with a custom preconditioner
+for the fault block associated with the Lagrange multipliers
+accelerates the convergence of the Krylov solver with the fewest
+number of iterations and the least sensitivity to problem
+size. \brad{Add something about poor scaling of AMG or good scaling of
+ ASM?} Benchmark tests demonstrate the accuracy of our implementation
+in PyLith with excellent agreement with the analytical solution for
+viscoelastic relaxation and strike-slip faulting over multiply
+earthquake cycle and excellent agreement with other codes for
+supershear spontaneous dynamic rupture on a dipping normal fault
+embedded in an elastoplastic domain. Consequently, we believe this
+methodology will provide a promising avenue for modeling the
+earthquake cycle through coupling of quasi-static simulations of the
+interseismic and postseismic deformation and dynamic rupture
+simulations of earthquake rupture propagation.
+
% ----------------------------------------------------------------------
% Notation -- End each entry with a period.
\begin{notation}
@@ -1567,14 +1561,11 @@ elastoplastic bulk rheology and slip-wea
% ------------------------------------------------------------------
\begin{acknowledgments}
-
-\begin{itemize}
- \item Funding from SCEC, CIG, GNS Science, USGS.
- \item SCEC contribution number.
- \end{itemize}
-
-MGK acknowledges partial support from NSF grant EAR-0949446.
-
+ Development of PyLith has been supported by the Earthquake Hazards
+ Program of the U.S. Geological Survey, the Computational
+ Infrastructure for Geodynamics (NSF grant EAR-0949446), GNS Science,
+ and the Southern California Earthquake Center. This is SCEC
+ contribution number ??.
\end{acknowledgments}
diff -r 77311386bf85 -r c3fd6d615a85 figs/savageprescott_profiles.pdf
Binary file figs/savageprescott_profiles.pdf has changed
More information about the CIG-COMMITS
mailing list