Implications of stress concentration on a strike-slip fault in an elastic plate subject to basal shear stress

Summary We consider a long strike-slip fault in a lithosphere modelled as an elastic slab. To the base of the slab a shear stress distribution is applied which simulates the viscous drag exerted by the asthenosphere. The resulant stress on the fault plane may directly fracture the lithosphere in its brittle upper portion; alternatively it may give rise at first to a stable aseismic sliding in the lower portion. In the latter case, stress concentration due to the deep aseismic slip is the relevant feature of the pre-seismic stress acting on the upper section of the lithosphere. The two cases are examined by use of dislocation theory and their observable effects compared. Different depths of the aseismic slip zone and the presence or absence of a uniform friction on the seismic fault are allowed for. If the model is applied to the San Andreas fault region, where a steady sliding condition actually seems to be present at shallow depth, it turns out that the slip amplitudes commonly associated with large earthquakes are consistent with average basal stress values which can be substantially lower than a few bars, a value often quoted as the steady state basal stress due to a velocity gradient in the upper asthenosphere.     

A Bayesian approach to estimating tectonic stress from seismological data

Earthquakes are conspicuous manifestations of tectonic stress, but the non-linear relationships between the stresses acting on a fault plane, its frictional slip, and the ensuing seismic radiation are such that a single earthquake by itself provides little information about the ambient state of stress. Moreover, observational uncertainties and inherent ambiguities in the nodal planes of earthquake focal mechanisms preclude straightforward inferences about stress being drawn on the basis of individual focal mechanism observations. However, by assuming that each earthquake in a small volume of the crust represents a single, uniform state of stress, the combined constraints imposed on that stress by a suite of focal mechanism observations can be estimated. Here, we outline a probabilistic (Bayesian) technique for estimating tectonic stress directions from primary seismological observations. The Bayesian formulation combines a geologically motivated prior model of the state of stress with an observation model that implements the physical relationship between the stresses acting on a fault and the resultant seismological observation. We show our Bayesian formulation to be equivalent to a well-known analytical solution for a single, errorless focal mechanism observation. The new approach has the distinct advantage, however, of including (1) multiple earthquakes, (2) fault plane ambiguities, (3) observational errors and (4) any prior knowledge of the stress field. Our approach, while computationally demanding in some cases, is intended to yield reliable tectonic stress estimates that can be confidently compared with other tectonic parameters, such as seismic anisotropy and geodetic strain rate observations, and used to investigate spatial and temporal variations in stress associated with major faults and coseismic stress perturbations.     

Source investigation of a small event using empirical Green's functions and simulated annealing

We propose a two-step inversion of three-component seismograms that (1) recovers the far-field source time function at each station and (2) estimates the distribution of co-seismic slip on the fault plane for small earthquakes (magnitude 3 to 4). The empirical Green's function (EGF) method consists of finding a small earthquake located near the one we wish to study and then performing a deconvolution to remove the path, site, and instrumental effects from the main-event signal.
The deconvolution between the two earthquakes is an unstable procedure: we have therefore developed a simulated annealing technique to recover a stable and positive source time function (STF) in the time domain at each station with an estimation of uncertainties. Given a good azimuthal coverage, we can obtain information on the directivity effect as well as on the rupture process. We propose an inversion method by simulated annealing using the STF to recover the distribution of slip on the fault plane with a constant rupture-velocity model. This method permits estimation of physical quantities on the fault plane, as well as possible identification of the real fault plane.
We apply this two-step procedure for an event of magnitude 3 recorded in the Gulf of Corinth in August 1991. A nearby event of magnitude 2 provides us with empirical Green's functions for each station. We estimate an active fault area of 0.02 to 0.15 km² and deduce a stress-drop value of 1 to 30 bar and an average slip of 0.1 to 1.6 cm. The selected fault of the main event is in good agreement with the existence of a detachment surface inferred from the tectonics of this half-graben.     

Dynamic stress field of a kinematic earthquake source model with k-squared slip distribution

Source models such as the k -squared stochastic source model with k -dependent rise time are able to reproduce source complexity commonly observed in earthquake slip inversions. An analysis of the dynamic stress field associated with the slip history prescribed in these kinematic models can indicate possible inconsistencies with physics of faulting. The static stress drop, the strength excess, the breakdown stress drop and critical slip weakening distance D _c distributions are determined in this study for the kinematic k -squared source model with k -dependent rise time. Several studied k -squared models are found to be consistent with the slip weakening friction law along a substantial part of the fault. A new quantity, the stress delay, is introduced to map areas where the yielding criterion of the slip weakening friction is violated. Hisada's slip velocity function is found to be more consistent with the source dynamics than Boxcar, Brune's and Dirac's slip velocity functions. Constant rupture velocities close to the Rayleigh velocity are inconsistent with the k -squared model, because they break the yielding criterion of the slip weakening friction law. The bimodal character of D _c / D _tot frequency–magnitude distribution was found. D _c approaches the final slip D _tot near the edge of both the fault and asperity. We emphasize that both filtering and smoothing routinely applied in slip inversions may have a strong effect on the space–time pattern of the inferred stress field, leading potentially to an oversimplified view of earthquake source dynamics.     

A kinematic viscoelastic model of the San Francisco earthquake of 1906

Summary. We construct a model of the San Andreas fault zone based on a rectangular fault in an elastic layer overlying a viscoelastic half-space. We alllow both steady and episodic aseismic slip at depth on the fault as well as a large-scale relative plate driving force. We use the model to explain the aseismic changes in geodetic triangulation angles observed during the 40 years following the 1906 San Francisco earthquake. The most important results are that viscoelastic relaxation can explain the data very well, and that the driving force of relative plate motion can be characterized by a horizontal distance scale perpendicular to the plate boundary of hundreds of kilometres.     

An improved method for the determination of the tectonic stress field from focal mechanism data

All conventional stress inversion methods, when applied to earthquake focal mechanism data, suffer from uncertainty as to which plane is the true fault plane. This paper deals with several problems in stress inversion brought about by this uncertainty. Our analysis shows that the direction of shear stress on the auxiliary plane does not coincide with the hypothetical slip direction unless the B -axis is parallel to one of the three principal stress directions. Based on this simple fact, we propose a new algorithm dealing with the ambiguity in fault/auxiliary plane identification. We also propose a method to handle the inhomogeneity problem of data quality, which is common and unique for focal mechanism data. Different inversion methods and algorithms are applied to two sets of 'focal mechanism' data simulated from field fault-slip measurement data. The inversion results show that, among the four stress parameters inverted, the stress ratio suffers the most from the ambiguity in fault/auxiliary plane identity, whereas the solutions for the principal stress directions are surprisingly good. The errors in inversion solutions resulting from the fault/auxiliary plane ambiguity can be significantly reduced by controlling subjectively the sample variance of the measurement errors. Our results also suggest that the fault plane cannot be distinguished correctly from the auxiliary plane with a high probability on the basis of the stress inversion alone.     

Earthquake mechanisms of the Adriatic Sea and Western Greece: implications for the oceanic subduction-continental collision transition

We present 21 focal solutions (magnitude > 5.5) reliably computed by body-wave modelling for the western Hellenic arc from Yugoslavia to the southern Peloponnese. Mechanisms located within the Aegean show normal faulting, the T-axis trending N-S in the centre and parallel to the active boundary in the external part. Mechanisms associated with the Keffalinia fault are consistent with dextral strike-slip motion. Reverse mechanisms located along the active boundary are remarkably consistent and do not depend on the nature of the active boundary (continental collision or oceanic subduction). The consistency in azimuth of the slip vectors and of the GPS velocity relative to Africa, all along the active boundary, suggests that the deformation is related to the same motion. The discrepancy between seismic-energy release and the amount of shortening confirms that the continental collision is achieved by seismic slip on faults but the oceanic subduction is partially aseismic. The northward decrease in velocity between continental collision and oceanic subduction suggests the continental collision to be a recent evolution of the active subduction.     

Source tomography by simulated annealing using broad-band surface waves and geodetic data: application to the Mw= 8.1 Chile 1995 event

We invert surface-wave and geodetic data for the spatio-temporal complexity of slip during the M _w =8.1 Chile 1995 event by simulated annealing. This quasi-global inversion method allows for a wide exploration of model space, and retains the non-linearity of the source tomography problem. Complex source spectra are obtained from 5 to 45 mHz from first- and second-orbit fundamental-mode Rayleigh waves using an empirical Green's function cross-correlation technique. Coseismic displacement vectors were measured at 10 GPS sites near Antofagasta. They are part of a French-Chilean experiment which monitors the Northern Chile seismic gap. The spectra, together with the geodetic data, are inverted for the moment distribution on a 2-D dipping fault, under the physical constraints of slip positivity and causality. Marginal a posteriori distributions of the model parameters are obtained from several independently inverted solutions. In general, features of the slip model are well resolved. Data are well fitted by a purely unilateral southward rupture with a nearly uniform velocity around 2.5–3.0 km s⁻¹, and a total duration of 65 s. Several regions of moment release were imaged, one near the hypocentre, a major one 80 km south of it and a minor one 160 km south of it. The major patch of moment release seemed to have propagated to relatively shallow depths near the trench, 100 km SSW of the epicentre. The region of major slip is located updip of the 1987, M _w =7.5 earthquake, suggesting a causal relationship. Most of the slip occurred updip of the hypocentre (36 km), but the entire coupled plate interface (20–40 km) ruptured during the Chile 1995 event.     

CMT data inversion using a Bayesian information criterion to estimate seismogenic stress fields

We developed an inversion method to estimate the stress fields related to earthquake generation (seismogenic stress fields) from the centroid moment tensors (CMT) of seismic events by using Akaike's Bayesian information criterion (ABIC). On the idea that the occurrence of an earthquake releases some part of the seismogenic stress field around its hypocentre, we define the CMT of a seismic event by a weighted volume integral of the true but unknown seismogenic stress field. Representing each component of the seismogenic stress field by the superposition of a finite number of 3-D basis functions (tri-cubic B-splines), we obtain a set of linear observation equations to be solved for the expansion coefficients (model parameters). We introduce prior constraint on the roughness of the seismogenic stress field and combine it with observed data to construct a Bayesian model with hierarchic, highly flexible structure controlled by hyper-parameters. The optimum values of the hyper-parameters are objectively determined form observed data by using ABIC. Given the optimum values of the hyper-parameters, we can obtain the best estimates of model parameters by using a maximum likelihood algorithm. We tested the validity of the inversion method through numerical experiments on two synthetic CMT data sets, assuming the distribution of fault orientations to be aligned with the maximum shear stress plane in one case and to be random in the other case. Then we applied the inversion method to actual CMT data in northeast Japan, and obtained the pattern of the seismogenic stress field consistent with geophysical and geological observations.     

A practical regularization for seismic tomography

Regularization is usually necessary in solving seismic tomographic inversion problems. In general the equation system of seismic tomography is very large, often making a suitable choice of the regularization parameter difficult. In this paper, we propose an algorithm for the practical choice of the regularization parameter in linear tomographic inversion. The algorithm is based on the types of statistical assumptions most commonly used in seismic tomography. We first transfer the system of equations into a Krylov subspace by using Lanczos bidiagonalization. In the transformed subspace, the system of equations is then changed into the form of a standard damped least squares normal equation. The solution to this normal equation can be written as an explicit function of the regularization parameter, which makes the choice of the regularization computationally convenient. Two criteria for the choice of the regularization parameter are investigated with the numerical simulations. If the dimensions of the transformed space are much less than that of the original model space, the algorithm can be very computationally efficient, which is practically useful in large seismic tomography problems.     

On the use of satellite altimetry to infer the earthquake rupture characteristics: application to the 2004 Sumatra event

Standard data and methods, such as the inversion of seismic and GPS data, have been used extensively to infer the details of the 2004 December 26 earthquake. The unprecedented large size of this event gave the opportunity to modern altimeters to provide the first clear records of a tsunami in deep ocean, therefore allowing us to study the rupture history from an independent perspective. We invert the Jason-1 and Topex–Poseidon altimetry records, considering the new constraints available on the geometry of the fault plane, and taking them into account in a 3-D rupture model. The data are corrected for the non-negligible effect of satellite motion during measurements. Our results show that the rupture propagated over the 1500 km of subduction zone initially identified by the aftershock distribution, with a magnitude of   M _w= 9.1  . Our solution compares well with the latitudinal distribution of slip inferred from other data sets, with a maximum of energy release north of Sumatra, and two other slip patches near the Nicobar and Andaman islands. Based on waveform comparison, we assert that the shallow portion of the megathrust offshore Banda Aceh had slip amplitudes of more than 20 m. Also, we find that significant amounts of slip (about 10 m) concentrated below the Andaman islands and did not propagate on the shallow portion of the interface. Although synthetic tests tend to show less resolution in the northern part of the rupture, this solution is compatible with the near-field data (GPS, coral heads and imagery), and would allow one to explain the apparent paradox between the large local displacements and the moderate tsunami observed locally. Finally, we demonstrate the rapidly dominating effect of propagation and slip distribution over the rupture velocity, and how it precludes the direct estimate of this latter parameter.     

Crosshole seismic waveform tomography – II. Resolution analysis

In an accompanying paper, we used waveform tomography to obtain a velocity model between two boreholes from a real crosshole seismic experiment. As for all inversions of geophysical data, it is important to make an assessment of the final model, to determine which parts of the model are well-resolved and can confidently be used for geological interpretation. In this paper we use checkerboard tests to provide a quantitative estimate of the performance of the inversion and the reliability of the final velocity model. We use the output from the checkerboard tests to determine resolvability across the velocity model. Such tests can act as good guides for designing appropriate inversion strategies. Here we discovered that, by including both reference-model and smoothing constraints in initial inversions, and then relaxing the smoothing constraint for later inversions, an optimum velocity image was obtained. Additionally, we noticed that the performance of the inversion was dependent on a relationship between velocity perturbation and checkerboard grid-size: larger velocity perturbations were better-resolved when the grid-size was also increased. Our results suggest that model assessment is an essential step prior to interpreting features in waveform tomographic images.     

Sequential integrated inversion of refraction and wide-angle reflection traveltimes and gravity data for two-dimensional velocity structures

A new algorithm is presented for the integrated 2-D inversion of seismic traveltime and gravity data. The algorithm adopts the 'maximum likelihood' regularization scheme. We construct a 'probability density function' which includes three kinds of information: information derived from gravity measurements; information derived from the seismic traveltime inversion procedure applied to the model; and information on the physical correlation among the density and the velocity parameters. We assume a linear relation between density and velocity, which can be node-dependent; that is, we can choose different relationships for different parts of the velocity–density grid. In addition, our procedure allows us to consider a covariance matrix related to the error propagation in linking density to velocity. We use seismic data to estimate starting velocity values and the position of boundary nodes. Subsequently, the sequential integrated inversion (SII) optimizes the layer velocities and densities for our models. The procedure is applicable, as an additional step, to any type of seismic tomographic inversion.
We illustrate the method by comparing the velocity models recovered from a standard seismic traveltime inversion with those retrieved using our algorithm. The inversion of synthetic data calculated for a 2-D isotropic, laterally inhomogeneous model shows the stability and accuracy of this procedure, demonstrates the improvements to the recovery of true velocity anomalies, and proves that this technique can efficiently overcome some of the limitations of both gravity and seismic traveltime inversions, when they are used independently.
An interpretation of field data from the 1994 Vesuvius test experiment is also presented. At depths down to 4.5 km, the model retrieved after a SII shows a more detailed structure than the model obtained from an interpretation of seismic traveltime only, and yields additional information for a further study of the area.     

High-frequency radiation from crack (stress drop) models of earthquake faulting

We present a theory for the radiation of high-frequency waves by earthquake faults. We model the fault as a planar region in which the stress drops to the kinematic friction during slip. This model is entirely equivalent to a shear crack. For two-dimensional fault models we show that the high frequencies originate from the stress and slip velocity concentrations in the vicinity of the fault's edges. These stress concentrations radiate when the crack expands with accelerated motion. The most efficient generation of high-frequency waves occurs when the rupture velocity changes abruptly. In this case, the displacement spectrum has an ω^-2 behaviour at high frequencies. The excitation is proportional to the intensity of the stress concentration near the crack tips and to the change in the focusing factor due to rupture velocity. We extend these two-dimensional results to more general three-dimensional fault models in the case when the rupture velocity changes simultaneously on the rupture front. Results are similar to those described for two-dimensional faults. We apply the theory to the case of a circular fault that grows at constant velocity and stops suddenly. The present theory is in excellent agreement with a numerical solution of the same problem.
Our results provide upper bounds to the high-frequency radiation from more realistic models in which rupture velocity does not change suddenly. The ω^-2 is the minimum possible decay at high frequencies for any crack model of the source.     

Dynamic rupture and stress change in a normal faulting earthquake in the subducting Cocos plate

A large nearly vertical, normal faulting earthquake ( M _w = 7.1) took place in 1997 in the Cocos plate, just beneath the ruptured fault zone of the great 1985 Michoacan thrust event ( M _w = 8.1). Dynamic rupture and resultant stress change during the 1997 earthquake have been investigated on the basis of near-source strong-motion records together with a 3-D dynamic model.
Dynamically consistent waveform inversion reveals a highly heterogeneous distribution of stress drop, including patch-like asperities and negative stress-drop zones. Zones of high stress drop are mainly confined to the deeper, southeastern section of the vertical fault, where the maximum dynamic stress drop reaches 280 bars (28 MPa). The dynamically generated source time function varies with location on the fault, and yields a short slip duration, which is caused by a short scalelength of stress-drop heterogeneities. The synthetic seismograms calculated from the dynamic model are generally consistent with the strong-motion velocity records in the frequency range lower than 0.5 Hz.
The pattern of stress-drop distribution appears, in some sense, to be consistent with that of coseismic changes in shear stress resulting from the 1985 thrust event. This consistency suggests that the stress transfer from the 1985 event to the subducting plate could be one of the possible mechanisms that increased the chance of the occurrence of the 1997 earthquake.     

Elastic full-waveform inversion for earthquake source parameters

2-D full-waveform inversion of double-couple earthquake sources is implemented. Temporally and spatially extended sources are represented by superposition of double-couples. The source parameters solved for are the spatial location, origin time, amplitude and orientation of each double-couple. The velocity and density distribution and source time function are assumed to be known a priori but may be arbitrarily complicated. The non-linear inverse problem is solved by iterative linear approximation. The Jacobian matrix elements for source depth and rupture angle are computed by wavefield extrapolation forward in time, while those for origin time and amplitude are computed analytically. A smoothing technique that results in faster convergence and avoids local minima associated with cycle skipping is applied at each iteration. A spatial sampling interval, between discrete sources, of one-quarter wavelength of the dominant shear wave is optimal for inversion if high uniqueness of the result is desired. The presence of a fault is inferred from the spatial continuity of the rupture solution, rather than being imposed a priori. The method is illustrated by successful application to three synthetic source models: a single double-couple, a single extended rupture and a double extended rupture. The resolutions of the source depth and origin time are higher, and their posterior covariances are lower than those of the amplitude and rupture angle at each source point. Source depth, origin time and amplitude are primarily determined by the data; the rupture angle is more strongly influenced by the a priori information.     

Active tectonics of an apparently aseismic region: distributed active strike-slip faulting in the Hangay Mountains of central Mongolia

We identify and describe a series of east–west left-lateral strike-slip faults (named the Songino-Margats, the Hag Nuur, the Uliastay and the South Hangay fault systems) in the Hangay mountains of central Mongolia: an area that has little in the way of recorded seismicity and which is often considered as a rigid block within the India–Eurasia collision zone. The strike-slip faults of central Mongolia constitute a previously unrecognized hazard in this part of Mongolia. Each of the strike-slip faults show indications of late Quaternary activity in the form of aligned sequences of sag-ponds and pressure-ridges developed in alluvial deposits. Total bed-rock displacements of ∼3 km are measured on both the Songino-Margats and South Hangay fault systems. Bed-rock displacements of 11 km are observed across the Hag Nuur fault. Cumulative offset across the Uliastay fault systems are unknown but are unlikely to be large. We have no quantitative constraint on the age of faulting in the Hangay. The ≤20 km of cumulative slip on the Hangay faults might, at least in part, be inherited from earlier tectonic movements. Our observations show that, despite the almost complete absence of instrumentally recorded seismicity in the Hangay, this part of Mongolia is cut through by numerous distributed strike-slip faults that accommodate regional left-lateral shear between Siberia and China. Central Mongolia is thus an important component of the India–Eurasia collision that would be overlooked in models of the active tectonics based on the distribution of seismicity. We suggest that active faults such as those identified in the Hangay of Mongolia might exist in other, apparently aseismic, regions within continental collision zones.     

Simultaneous velocity and interface tomography of normal-incidence and wide-aperture seismic traveltime data

A tomographic inversion technique that inverts traveltimes to obtain a model of the subsurface in terms of velocities and interfaces is presented. It uses a combination of refraction, wide-angle reflection and normal-incidence data, it simultaneously inverts for velocities and interface depths, and it is able to quantify the errors and trade-offs in the final model. The technique uses an iterative linearized approach to the non-linear traveltime inversion problem. The subsurface is represented as a set of layers separated by interfaces, across which the velocity may be discontinuous. Within each layer the velocity varies in two dimensions and has a continuous first derivative. Rays are traced in this medium using a technique based on ray perturbation theory, and two-point ray tracing is avoided by interpolating the traveltimes to the receivers from a roughly equidistant fan of rays. The calculated traveltimes are inverted by simultaneously minimizing the misfit between the data and calculated traveltimes, and the roughness of the model. This 'smoothing regularization' stabilizes the solution of the inverse problem. In practice, the first iterations are performed with a high level of smoothing. As the inversion proceeds, the level of smoothing is gradually reduced until the traveltime residual is at the estimated level of noise in the data. At this point, a minimum-feature solution is obtained, which should contain only those features discernible over the noise.
The technique is tested on a synthetic data set, demonstrating its accuracy and stability and also illustrating the desirability of including a large number of different ray types in an inversion.     

Global crustal thickness from neural network inversion of surface wave data

We present a neural network approach to invert surface wave data for a global model of crustal thickness with corresponding uncertainties. We model the a posteriori probability distribution of Moho depth as a mixture of Gaussians and let the various parameters of the mixture model be given by the outputs of a conventional neural network. We show how such a network can be trained on a set of random samples to give a continuous approximation to the inverse relation in a compact and computationally efficient form. The trained networks are applied to real data consisting of fundamental mode Love and Rayleigh phase and group velocity maps. For each inversion, performed on a 2°× 2° grid globally, we obtain the a posteriori probability distribution of Moho depth. From this distribution any desired statistic such as mean and variance can be computed. The obtained results are compared with current knowledge of crustal structure. Generally our results are in good agreement with other crustal models. However in certain regions such as central Africa and the backarc of the Rocky Mountains we observe a thinner crust than the other models propose. We also see evidence for thickening of oceanic crust with increasing age. In applications, characterized by repeated inversion of similar data, the neural network approach proves to be very efficient. In particular, the speed of the individual inversions and the possibility of modelling the whole a posteriori probability distribution of the model parameters make neural networks a promising tool in seismic tomography.