We extend the full wavefield modeling with forward scattering theory and Volterra Renormalization to a vertically varying two-parameter (velocity and density) acoustic medium. The forward scattering series, derived by applying Born-Neumann iterative procedure to the Lippmann-Schwinger equation (LSE), is a well known tool for modeling and imaging. However, it has limited convergence properties depending on the strength of contrast between the actual and reference medium or the angle of incidence of a plane wave component. Here, we introduce the Volterra renormalization technique to the LSE. The renormalized LSE and related Neumann series are absolutely convergent for any strength of perturbation and any incidence angle. The renormalized LSE can further be separated into two sub-Volterra type integral equations, which are then solved noniteratively. We apply the approach to velocity-only, density-only, and both velocity and density perturbations. We demonstrate that this Volterra Renormalization modeling is a promising and efficient method. In addition, it can also provide insight for developing a scattering theory-based direct inversion method.

We develop a super-grid modeling technique for solving the elastic wave equation in semi-bounded two- and three-dimensional spatial domains. In this method, waves are slowed down and dissipated in sponge layers near the far-field boundaries. Mathematically, this is equivalent to a coordinate mapping that transforms a very large physical domain to a significantly smaller computational domain, where the elastic wave equation is solved numerically on a regular grid. To damp out waves that become poorly resolved because of the coordinate mapping, a high order artificial dissipation operator is added in layers near the boundaries of the computational domain. We prove by energy estimates that the super-grid modeling leads to a stable numerical method with decreasing energy, which is valid for heterogeneous material properties and a free surface boundary condition on one side of the domain. Our spatial discretization is based on a fourth order accurate finite difference method, which satisfies the principle of summation by parts. We show that the discrete energy estimate holds also when a centered finite difference stencil is combined with homogeneous Dirichlet conditions at several ghost points outside of the far-field boundaries. Therefore, the coefficients in the finite difference stencils need only be boundary modified near the free surface. This allows for improved computational efficiency and significant simplifications of the implementation of the proposed method in multi-dimensional domains. Numerical experiments in three space dimensions show that the modeling error from truncating the domain can be made very small by choosing a sufficiently wide super-grid damping layer. The numerical accuracy is first evaluated against analytical solutions of Lamb’s problem, where fourth order accuracy is observed with a sixth order artificial dissipation. We then use successive grid refinements to study the numerical accuracy in the more complicated motion due to a point moment tensor source in a regularized layered material.

In this paper, we propose a strong stability-preserving predictor-corrector (SSPC) method based on an implicit Runge-Kutta method to solve the acoustic- and elastic-wave equations. We first transform the wave equations into a system of ordinary differential equations (ODEs) and apply the local extrapolation method to discretize the spatial high-order derivatives, resulting in a system of semi-discrete ODEs. Then we use the SSPC method based on an implicit Runge-Kutta method to solve the semi-discrete ODEs and introduce a weighting parameter into the SSPC method. On top of such a structure, we develop a robust numerical algorithm to effectively suppress the numerical dispersion, which is usually caused by the discretization of wave equations when coarse grids are used or geological models have large velocity contrasts between adjacent layers. Meanwhile, we investigate the performance of the SSPC method including numerical errors and convergence rate, numerical dispersion, and stability criteria with different choices of the weighting parameter to solve 1-D and 2-D acoustic- and elastic-wave equations. When the SSPC is applied to seismic simulations, the computational efficiency is also investigated by comparing the SSPC, the fourth-order Lax-Wendroff correction (LWC) method, and the staggered-grid (SG) finite difference method. Comparisons of synthetic waveforms computed by the SSPC and analytic solutions for acoustic and elastic models are given to illustrate the accuracy and the validity of the SSPC method. Furthermore, several numerical experiments are conducted for the geological models including a 2-D homogeneous transversely isotropic (TI) medium, a two-layer elastic model, and the 2-D SEG/EAGE salt model. The results show that the SSPC can be used as a practical tool for large-scale seismic simulation because of its effectiveness in suppressing numerical dispersion even in the situations such as coarse grids, strong interfaces, or high frequencies.

We develop a stable finite difference approximation of the three-dimensional viscoelastic wave equation. The material model is a super-imposition of N standard linear solid mechanisms, which commonly is used in seismology to model a material with constant quality factor Q. The proposed scheme discretizes the governing equations in second order displacement formulation using 3N memory variables, making it significantly more memory efficient than the commonly used first order velocity-stress formulation. The new scheme is a generalization of our energy conserving finite difference scheme for the elastic wave equation in second order formulation [SIAM J. Numer. Anal., 45 (2007), pp. 1902-1936]. Our main result is a proof that the proposed discretization is energy stable, even in the case of variable material properties. The proof relies on the summation-by-parts property of the discretization. The new scheme is implemented with grid refinement with hanging nodes on the interface. Numerical experiments verify the accuracy and stability of the new scheme. Semi-analytical solutions for a half-space problem and the LOH.3 layer over half-space problem are used to demonstrate how the number of viscoelastic mechanisms and the grid resolution influence the accuracy. We find that three standard linear solid mechanisms usually are sufficient to make the modeling error smaller than the discretization error.

To finite-difference model elastic wave propagation in a combined structure with solid, fluid and porous subregions, a set of modified Biot’s equations are used, which can be reduced to the governing equations in solids, fluids as well as fluid-saturated porous media. Based on the modified Biot’s equations, the field quantities are finite-difference discretized into unified forms in the whole structure, including those on any interface between the solid, fluid and porous subregions. For the discrete equations on interfaces, however, the harmonic mean of shear modulus and the arithmetic mean of the other parameters on both sides of the interfaces are used. These parameter averaging equations are validated by deriving from the continuity conditions on the interfaces. As an example of using the parameter averaging technique, a 2-D finite-difference scheme with a velocity-stress staggered grid in cylindrical coordinates is implemented to simulate the acoustic logs in porous formations. The finite-difference simulations of the acoustic logging in a homogeneous formation agree well with those obtained by the analytical method. The acoustic logs with mud cakes clinging to the borehole well are simulated for investigating the effect of mud cake on the acoustic logs. The acoustic logs with a varying radius borehole embedded in a horizontally stratified formation are also simulated by using the proposed finite-difference scheme.

Spectral element methods are well established in the field of wave propagation, in particular because they inherit the flexibility of finite element methods and have low numerical dispersion error. The latter is experimentally acknowledged, but has been theoretically shown only in limited cases, such as Cartesian meshes. It is well known that a finite element mesh can contain distorted elements that generate numerical errors for very large distortions. In the present work, we study the effect of element distortion on the numerical dispersion error and determine the distortion range in which an accurate solution is obtained for a given error tolerance. We also discuss a double-grid calculation of the spectral element matrices that preserves accuracy in deformed geometries.

This is a review of progress made since [R. McKibbin, “An attempt at modelling hydrothermal eruptions”, Proc. 11th New Zealand Geothermal Workshop 1989 (University of Auckland, 1989), 267–273] began development of a mathematical model for progressive hydrothermal eruptions (as distinct from “blasts”). Early work concentrated on modelling the underground process, while lately some attempts have been made to model the eruption jet and the flight and deposit of ejected material. Conceptually, the model is that of a boiling and expanding two-phase fluid rising through porous rock near the ground surface, with a vertical high-speed jet, dominated volumetrically by the gas phase, ejecting rock particles that are then deposited on the ground near the eruption site. Field observations of eruptions in progress and experimental results from a laboratory-sized model have confirmed the conceptual model. The quantitative models for all parts of the process are based on the fundamental conservation equations of motion and thermodynamics, using a continuum approximation for each of the components.

In this paper a limit theorem is proved that establishes conditions under which the distribution of the difference between the size of the initial event in a random sequence, modeled as a finite point process, and the largest subsequent event approaches a limiting form independent of the size of the initial event. The underlying assumptions are that the sizes of the individual events follow an exponential distribution, that the expected total number of events increases exponentially with the size of the initial event, and that the structure of the sequence approximates that of a Poisson process. Particular cases to which the results apply include sequences of independent and identically distributed exponential variables, and the epidemic-type aftershock (ETAS) branching process model in the subcritical case. In all these cases the form of the limit distribution is shown to be that of a double exponential (type-I extreme value distribution). In sampling from a family of aftershock sequences, with possibly different underlying parameters, the limit distribution is a mixture of such double exponential distributions. The conditions for the simple limit to exist relate to the approximation of the distribution of the number of events by a Poisson distribution. One such condition requires the coefficient of variation (ratio of standard deviation to mean) of the number of events to converge to 0 as the mean increases. The results provide a statistical background to Båth's law in seismology, which asserts that in an aftershock sequence the magnitude of the main shock is commonly around 1.2 units higher than the magnitude of the largest aftershock.

The paper considers the superposition of modified Omori functions as a conditional intensity function for a point process model used in the exploratory analysis of earthquake clusters. For the examples discussed, the maximum likelihood estimates converge well starting from appropriate initial values even though the number of parameters estimated can be large (though never larger than the number of observations). Three datasets are subjected to different analyses, showing the use of the model to discover and study individual clustering features.

This paper uses the epidemic-type aftershock sequence (ETAS) point process model to study certain seismicity features of the Jiashi swarm of certain earthquakes, investigating in particular whether there is relative quiescence prior to the quite big events within the Jiashi sequence. The seven earthquake sequences studied occurred in the region of Jiashi, south of Tianshan Mountain, Xinjiang, China. The particular ETAS model that is developed is consistent with the reality of seismic activity. The various features of Jiashi swarm activity can be described as focusing in different stages. There is obvious precursory quiescence prior to most big events with Ms ≥ 6.0 within the Jiashi swarm. Thus, checking for relative quiescence can be use for earthquake prediction.

For many years the modified Mercalli (MM) scale has been used to describe earthquake damage and effects observed at scattered locations. In the next stage of an analysis involving MM data, isoseismal lines based on the observations have been added to maps by hand, i.e. subjectively. However a few objective methods have been proposed (by e.g. De Rubeis et al., Brillinger, Wald et al. and Pettenati et al.). The work presented here develops objective methods further. In particular the ordinal character of the MM scale is specifically taken into account. Numerical smoothing is basic to the approach and methods involving splines, local polynomial regression and wavelets are illustrated. The approach also allows the inclusion of explanatory variables, for example site effects. The procedure is implemented for data from the 17 October 1989 Loma Prieta earthquake.

A point process procedure can be used to study reservoir-induced seismicity (RIS), in which the intensity function representing earthquake hazard is a combination of three terms: a constant background term, an ETAS (epidemic-type aftershock sequence) term for aftershocks, and a time function derived from observation of water levels of a reservoir. This paper presents the results of such a study of the seismicity in the vicinity of the Tarbela reservoir in Pakistan. Making allowance for changes in detection capability and the background seismicity related to tectonic activity, earthquakes of magnitude ≥ 2.0, occurring between May 1978 and January 1982 and whose epicentres were within 100 km of the reservoir, were used in this analysis. Several different intensities were compared via their Akaike information criterion (AIC) values relative to those of a Poisson process. The results demonstrate that the seismicity within 20 km of the reservoir correlates with water levels of the reservoir, namely, active periods occur about 250 days after the appearance of low water levels. This suggests that unloading the reservoir activates the seismicity beneath it. Seasonal variations of the seismicity in an area up to 100 km from the reservoir were also found, but these could not be adequately interpreted by an appropriate RIS mechanism.

The paper reviews the formulation of the linked stress release model for large scale seismicity together with aspects of its application. Using data from Taiwan for illustrative purposes, models can be selected and verified using tools that include Akaike's information criterion (AIC), numerical analysis, residual point processes and Monte Carlo simulation.

The annual earthquake predictions of the China Seismological Bureau (CSB) are evaluated by means of an R score (an R score is approximately 0 for completely random guesses, and approximately 1 for completely successful predictions). The average R score of the annual predictions in China in the period 1990–1998 is about 0.184, significantly larger than 0.0. However, background seismicity is higher in seismically active regions. If a ‘random guess' prediction is chosen to be proportional to the background seismicity, the expected R score is 0.123, and the nine-year mean R score of 0.184 as observed is only marginally higher than this background value. Monte Carlo tests indicate that the probability of attaining an R score of actual prediction by background seismicity based on random guess is about . It is concluded that earthquake prediction in China is still in a very preliminary stage, barely above a pure chance level.