Bayesian Inference of Lower-Crustal Viscosity Near the Kunlun Fault Based on Geologic, Geomorphic, and Geodetic Data

George E. Hilley, Roland Bürgmann, Pei-Zhen Zhang (State Key Laboratory of Earthquake Dynamics, Beijing, China), and Peter Molnar (University of Colorado, Boulder, USA)

Introduction

Bayesian methods provide the means of integrating geologic, paleoseismic, seismic, and geodetic data to improve estimates of fault zone and lower-crustal properties (e.g., Segall, 2002) that are important for predicting the long-term deformation of plate boundary zones (e.g., Shen et al., 2001), understanding deformation and stress transfer following earthquakes, and estimating future seismic hazard (e.g., Hardebeck, 2004). Herein, we develop a Bayesian methodology that integrates geologic, geomorphic, and geodetic data to provide probabilistic estimates of fault-zone and lower-crustal properties.

We will apply this methodology to the Kunlun Fault in northern Tibet (Figure 24.1), where estimates of lower-crustal viscosity are generally lacking. Here, geologic and geomorphic information constrain the range of permissible long-term slip rates and coseismically generated offsets. We combined these a priori estimates with GPS velocities using a Bayesian implementation of an elastic-viscoelastic earthquake cycle model (Savage and Prescott, 1978) to estimate fault-zone and lower crustal properties in the area. The probabilistic nature of these models allows straightforward assessments of the uncertainties within, and covariance between model parameters. Our example shows that these methods may aid in elucidating the active tectonics of other areas, and improve seismic hazard assessments by formally assimilating disparate data types (e.g., geologic, geomorphic, seismic, or geodetic datasets) into crustal deformation studies.

Bayesian Modeling Framework

The Bayesian approach that we use to analyze data from northern Tibet uses geologic and geomorphic estimates of fault slip rates and observed offsets in conjunction with geodetic data to refine estimates of fault slip rate, schizosphere thickness, recurrence time of events, and viscous relaxation time of the plastosphere (``chizosphere'' refers to the portion of the upper crust that deforms elastically, while ``plastosphere'' refers to areas that deform elastically over short time-scales but undergo viscous stress relaxation over longer time-scales; Scholz, 1988). Our method is based on Bayes' Rule (Bayes, 1763), which allows quantitative refinement of initial estimates of model parameters (in this case, fault zone and lower crustal properties) given the information provided by the geodetic observations (Segall, 2002; Johnson and Segall, 2004):

$$P(m_{i}\vert x)=\frac{P(x\vert m_{i})\times P(m_{i})}{\sum_{i=1}^{n} P(x\vert m_{i})\times P(m_{i})}$$ (25.1)

In the context of this paper, $m_{i}$ denotes a vector of model parameters, $x$ is a vector of the observed geodetic velocities, $P(m_{i}\vert x)$ denotes the probability that the set of model parameters explains the geodetic velocities, $P(x\vert m_{i})$ is the probability of observing the geodetic velocities given a combination of model parameters $m_{i}$, $P(m_{i})$ is the probability that the chosen combination of model parameters actually occur, and the denominator normalizes the probability to all possible combinations of model parameters. In this context, $P(x\vert m_{i})$ encapsulates the goodness of fit between observed geodetic velocities and those that would be expected based on the physical model that relates the surface velocity distribution to fault-zone and lower crustal properties (described below). Geodetic inversions that find the set of model parameters that best explain the observed data maximize $P(x\vert m_{i})$. In contrast, the Bayesian approach not only considers the goodness of fit of the velocity to model predictions, but also $P(m_{i})$, which describes in a probabilistic sense, some prior knowledge about what sets of model parameters are likely to occur. This term may be used to incorporate information such as slip rates along faults and recurrence times of earthquakes whose range may be estimated based on geologic and geomorphic considerations (e.g., Buck et al., 1996; Biasi et al., 2002; Segall, 2002). Therefore, the Bayesian modeling strategy has a two-fold advantage to conventional geodetic inversions: 1) Uncertainties in model parameters and covariance between model parameters may be straightforwardly obtained for complicated models (e.g., Hargreves and Annan, 2002); and 2) Other types of data (e.g., geologic, geomorphic, and paleoseismic) may be used in conjunction with the geodetic data to improve estimates of model parameters. Explicit evaluation of Equation 24.1 may be untenable for potentially large, multidimensional parameter spaces that may arise even in simple physical models such as the model described below. Therefore, we solve Equation 1 using the Metropolis-Hastings variant of the Markov-Chain Monte Carlo (MCMC) simulation methods that allows approximation of the joint distribution $P(m_{i}\vert x)$ without an exhaustive search of the parameter space (Metropolis et al., 1953).

Figure 24.1: Location of the Kunlun Fault in central Asia. The Kunlun fault is one of a series of strike-slip faults that may accommodate a portion of the Indo-Eurasian collision. White line shows approximate location of GPS profile used in this study.

The mechanical model we employ (evaluated through $P(x\vert m_{i})$) idealizes earthquake-cycle deformation as resulting from elastic strain release during rupturing events followed by viscoelastic relaxation of the lower crust and upper mantle (e.g., Savage and Prescott, 1978). The crust is idealized as two-dimensional in cross-section, and so movement along the major strike-slip faults occurs in the out-of-plane dimension. The rheology of the crust is treated as a two-layered medium in which an elastically deforming schizosphere overlies a linear, visco-elastically deforming plastosphere (e.g., Savage and Prescott, 1978; Savage, 2000; Segall, 2002; Dixon et al., 2003). The deformation during the earthquake cycle is defined by the long-term, average strike-slip velocity of the modeled fault ($s$), the Maxwell relaxation time of the plastosphere ( $\tau = 2\nu/\mu$, where $\nu$ is the viscosity and $\mu$ is the shear modulus), the thickness of the schizosphere ($H$), the time since the last earthquake ($t$), and the average recurrence interval ($T$). Following Savage and Prescott (1978) and Segall (2002), we non-dimensionalize the model by introducing four groups: $t^{*} = t/T$, $\tau^{*} = \tau/T$ , $s^{*} = s/s_{ref}$, and $H^{*} = H/x_{ref}$, where $s_{ref}$ and $x_{ref}$ are arbitrary velocity and length scales set to 1 mm/yr and 1 km, respectively. Although this is a simple representation of the earthquake cycle, this model is appropriate for inferring slip rates along the Kunlun Fault (and active strike-slip faults in Tibet in general) because: 1) The model captures the basic observation that the crustal deformation during the seismic cycle consists of coseismic, elastic, deformation followed by postseismic viscoelastic adjustment of the lower crust (e.g., Savage and Prescott, 1978; Segall, 2002); and 2) Several previous studies have successfully employed these types of simple models to provide first-order estimates of fault zone and lower-crustal properties. Given the density of geodetic measurements from this portion of Tibet, the earthquake cycle modeling we employ allows us to consider velocity variations during the earthquake cycle that cannot be assessed using the types of dislocation models currently employed to interpret interseismic geodetic velocities from this area (i.e., Wallace et al., 2004; Zhang et al., 2004).

Application to the Kunlun Fault

Along the Kunlun Fault, numerous studies show consistency in estimated slip rates over a variety of different time-scales (10-100 kyr time-scales), and so the primary goal of the specific example is to estimate lower-crustal properties beneath the Kunlun fault based on geologic, geomorphic, and geodetic datasets. Along the Kunlun fault, 1) Geologic and geomorphic studies indicate that long-term slip rates are between 9-16 mm/yr (Kidd and Molnar, 1988; van der Woerd et al., 2001); 2) Geomorphic studies and historic ruptures suggest that the slip during each event is approximately 10 m (van der Woerd et al., 2001); and 3) The lack of historic seismicity places a minimum bound on the recurrence time of earthquakes along the fault. Using this information, we will construct prior joint probability densities and apply our Bayesian methodology to refine estimates of each of these parameters and the Maxwell relaxation time of the lower crust.

References

Bayes, T., An essay towards solving a problem in the doctrine of chances, Phil. Trans. R. Soc. Lond., 53, 370-418, 1763.

Biasi, G. P., R. J. Weldon II, T. E. Fumal, and G. G. Seitz, Paleoseismic event dating and the conditional probability of large earthquakes on the southern San Andreas Fault, California, Bull. Seism. Soc. Am., 92, 2761-2781, 2002.

Buck, C. E., W. G. Cavanagh, and C. D. Litton, The Bayesian Approach to Interpreting Archaeological Data, Wiley, Chichester, 382 pp, 1996.

Dixon, T. H., E. Norabuena, and L. Hotaling, Paleoseismology and Global Positioning System: Earthquake-cycle effects and geodetic versus geologic fault slip rates in the Eastern California shear zone, Geology, 31, 55-58, 2003.

Hardebeck, J. L., Stress triggering and earthquake probability estimates, J. Geophys. Res., B109, B04310, doi:10.1029/2003JB002437, 2004.

Hargreaves, J. C., and J. D. Annan, Assimilation of paleo-data in a simple Earth system model, Clim. Dynam., 19, 371-381, 2002.

Johnson, K. M., and P. Segall, Viscoelastic earthquake cycle models of deep stress-driven creep along the San Andreas Fault system, J. Geophys. Res., B, in press, 2004.

Kidd, W. S. F. and P. Molnar, Quaternary and active faulting observed on the 1985 Academia Sinica-Royal Society Geotraverse of Tibet, Phil. Trans. R. Soc. Lond., 327, 337-363, 1988.

Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equations of state calculations by fast computing machines, J. Chem. Phys., 21, 1087-1091, 1953.

Savage, J. C., Viscoelastic-coupling model for the earthquake cycle driven from below, J. Geophys. Res., B105, 25,525-25,532, 2000.

Savage, J. C. and W. H. Prescott, Asthenosphere readjustment and the earthquake cycle, J. Geophys. Res., 83, 3369-3376, 1978.

Scholz, C.H., The mechanics of earthquakes and faulting, Cambridge University Press, Cambridge, 439 pp, 1988.

Segall, P., Integrating geologic and geodetic estimates of slip rate on the San Andreas fault system, Int. Geol. Rev., 44, 62-82, 2002.

Shen, F., L. H. Royden, and B. C. Burchfiel, Large-scale crustal deformation of the Tibetan Plateau, J. Geophys. Res., B106, 6793-6816, 2001.

van der Woerd, J., P. Tapponnier, F. J. Ryerson, A.-S. Meriaux, B. Meyer, Y. Gaudemer, R. C. Finkel, M. W. Caffee, G.h. Zhao, and Z.q. Xu, Uniform postglacial slip-rate along the central 600 km of the Kunlun Fault (Tibet), from 26Al, 10Be, and 14C dating of riser offsets, and climatic origin of the regional morphology, Geophys. J. Int., 148, 356-388, 2001.

Wallace, K., G.h.Yin, and R. Bilham, Inescapable slow slip on the Altyn Tagh fault, Geophys. Res. Lett., 31, L09613, doi:10.1029/2004GH019724, 2004.

Zhang, P.-Z., Z. Shen, M. Wang, W. Gan, R. Bürgmann, P. Molnar, Z. Niu, J. Sun, Q. Wang, J. Wu, H. Sun, and X. You, Continuous deformation of the Tibetan Plateau from GPS, Geology, in press, 2004.