DERIVATION OF LAYER PARAMETERS OF AN ELASTIC MEDIUM FROM REFLECTION COEFFICIENT MATRICES*

In this paper we develop a recursive algorithm to obtain the layer parameters of an elastic medium (density, P-wave velocity, S-wave velocity) from reflection coefficient matrices in terms of energy flux ratios for a non-normal incidence case. We define a layer impedance matrix, analogous to the impedance of an acoustic medium. Next we derive a matrix relationship between the layer impedance matrix of the n+ 1st layer and the reflection coefficient and parameter matrices of the nth layer. This relationship leads to recursively computing the parameters of the subsurface. We show that the elastic case—unlike the acoustic case—allows one to recover the layer parameters from the impedance matrix for non-normal incidence. The results of this work play a key role in the solution of the inverse problem with non-normal-incidence plane-wave seismic data when using a downward continuation technique.     

Dynamic stiffness matrices for two circular foundations

A systematic procedure is presented for generating dynamic stiffness matrices for two independent circular foundations on an elastic half-space medium. With the technique reported in References 1–3, the analytic solution of three-dimensional (3D) wave equations satisfying the prescribed traction due to the vibration of one circular foundation can be found. Since there are two analytic solutions for two prescribed tractions due to the vibrations of two circular foundations, the principle of superposition must be used to obtain the total solution. The interaction stresses (prescribed tractions) are assumed to be piecewise linear in the r-directions of both cylindrical co-ordinates for the two circular foundations. Then, the variational principle and the reciprocal theorem are employed to generate the dynamic stiffness matrices for the two foundations. In the process of employing the variational principle, a co-ordinate transformation matrix between two cylindrical co-ordinate systems is introduced. Some numerical results of dynamic stiffness matrices for the interaction of two identical rigid circular foundations are presented in order to show the effectiveness and efficiency of the present method, and some elaborations for its future extensions are also discussed.     

Lanczos method for dynamic analysis of damped structural systems

In this paper we extend the Lanczos algorithm for the dynamic analysis of structures⁷ to systems with general matrix coefficients. The equations of dynamic equilibrium are first transformed to a system of first order differential equations. Then the unsymmetric Lanczos method is used to generate two sets of vectors. These vectors are used in a method of weighted residuals to reduce the equations of motion to a small unsymmetric tridiagonal system. The algorithm is further simplified for systems of equations with symmetric matrices. By appropriate choice of the starting vectors we obtain an implementation of the Lanczos method that is remarkably close to that in Reference 7, but generalized to the case with indefinite matrix coefficients. This simplification eliminates one of the sets of vectors generated by the unsymmetric Lanczos method and results in a symmetric tridiagonal, but indefinite, system. We identify the difficulties that may arise when this implementation is applied to problems with symmetric indefinite matrices such as vibration of structures with velocity feedback control forces which lead to symmetric damping matrices. This approach is used to evaluate the vibration response of a damped beam problem and a space mast structure with symmetric damping matrix arising from velocity feedback control forces. In both problems, accurate solutions were obtained with as few as 20 Lanczos vectors.     

On recurrence relations for multipole coefficients

Summary General recurrence relations between the coefficients in thenth and (n+1)th order spherical harmonic multipole expansions are derived. The particular application presented here is the derivation of the equations concerned with representing the geomagnetic field by magnetic multipoles. The equations up to the 3rd order multipole are given as an example of the method. The main advantage in using these recurrence relations rather than other methods is that the mathematics is reduced to merely a matter of successive substitutions and this allows a fast step by step generation of the required equations, in a form for which there is a simple numerical program for solution.     

Fast algorithms for automatic mapping with space-limited covariance functions

In this paper we discuss a fast Bayesian extension to kriging algorithms which has been used successfully for fast, automatic mapping in emergency conditions in the Spatial Interpolation Comparison 2004 (SIC2004) exercise. The application of kriging to automatic mapping raises several issues such as robustness, scalability, speed and parameter estimation. Various ad-hoc solutions have been proposed and used extensively but they lack a sound theoretical basis. In this paper we show how observations can be projected onto a representative subset of the data, without losing significant information. This allows the complexity of the algorithm to grow as O(n m ²), where n is the total number of observations and m is the size of the subset of the observations retained for prediction. The main contribution of this paper is to further extend this projective method through the application of space-limited covariance functions, which can be used as an alternative to the commonly used covariance models. In many real world applications the correlation between observations essentially vanishes beyond a certain separation distance. Thus it makes sense to use a covariance model that encompasses this belief since this leads to sparse covariance matrices for which optimised sparse matrix techniques can be used. In the presence of extreme values we show that space-limited covariance functions offer an additional benefit, they maintain the smoothness locally but at the same time lead to a more robust, and compact, global model. We show the performance of this technique coupled with the sparse extension to the kriging algorithm on synthetic data and outline a number of computational benefits such an approach brings. To test the relevance to automatic mapping we apply the method to the data used in a recent comparison of interpolation techniques (SIC2004) to map the levels of background ambient gamma radiation.

---

Partial eigensolution of damped structural systems by Arnoldi's method

An efficient numerical algorithm is developed to solve the quadratic eigenvalue problems arising in the dynamic analysis of damped structural systems. The algorithm can even be applied to structural systems with non-symmetric matrices. The algorithm is based on the use of Arnoldi's method to generate a Krylov subspace of trial vectors, which is then used to reduce a large eigenvalue problem to a much smaller one. The reduced eigenvalue problem is solved and the solutions are used to construct approximate solutions to the original large system. In the process, the algorithm takes full advantage of the sparseness and symmetry of the system matrices and requires no complex arithmetic, therefore, making it very economical for use in solving large problems. The numerical results from test examples are presented to demonstrate that a large fraction of the approximate solutions calculated are very accurate, indicating that the algorithm is highly effective for extracting a number of vibration modes for a large dynamic system, whether it is lightly or heavily damped.     

COMPUTATION OF INTERVAL VELOCITIES FROM COMMON REFLECTION POINT MOVEOUT TIMES FOR n LAYERS WITH ARBITRARY DIPS AND CURVATURES IN THREE DIMENSIONS WHEN ASSUMING SMALL SHOT-GEOPHONE DISTANCES*

It is well known that interval velocities can be determined from common-reflection-point moveout times. However, the mathematics becomes complicated in the general case of n homogeneous layers with curved interfaces dipping in three dimensions. In this paper the problem is solved by mathematical induction using the second power terms only of the Taylor series which represents the moveout time as a function of the coordinate differences between shot and geophone points. Moreover, the zero-offset reflection times of the nth interface in a certain area surrounding the point of interest have to be known. The n—I upper interfaces and interval velocities are known too on account of the mathematical induction method applied. Thus, the zero-offset reflection raypath of the nth interface can be supposed to be known down to the intersection with the (n—1)th interface. The method applied consists mainly in transforming the second power terms of the moveout time from one interface to the next one. This is accomplished by matrix algebra. Some special cases are discussed as e.g. uniform strike and small curvatures.     

On the informative value of the largest sample element of log-Gumbel distribution

Extremes of stream flow and precipitation are commonly modeled by heavytailed distributions. While scrutinizing annual flow maxima or the peaks over threshold, the largest sample elements are quite often suspected to be low quality data, outliers or values corresponding to much longer return periods than the observation period. Since the interest is primarily in the estimation of the right tail (in the case of floods or heavy rainfalls), sensitivity of upper quantiles to largest elements of a series constitutes a problem of special concern. This study investigated the sensitivity problem using the log-Gumbel distribution by generating samples of different sizes (n) and different values of the coefficient of variation by Monte Carlo experiments. Parameters of the log-Gumbel distribution were estimated by the probability weighted moments (PWMs) method, method of moments (MOMs) and maximum likelihood method (MLM), both for complete samples and the samples deprived of their largest elements. In the latter case, the distribution censored by the non-exceedance probability threshold, F _T , was considered. Using F _T instead of the censored threshold T creates possibility of controlling estimator property. The effect of the F _T value on the performance of the quantile estimates was then examined. It is shown that right censoring of data need not reduce an accuracy of large quantile estimates if the method of PWMs or MOMs is employed. Moreover allowing bias of estimates one can get the gain in variance and in mean square error of large quantiles even if ML method is used.     

Surface Deformation due to Loading of a Layered Elastic Half-space: Constructing the Solution for a General Polygonal Load

We describe an algorithm for rapidly computing the surface displacements induced by a general polygonal load on a layered, isotropic, elastic half-space. The arbitrary surface pressure field is discretized using a large number, n, of equally-sized circular loading elements. The problem is to compute the displacement at a large number, m, of points (or stations) distributed over the surface. The essence of our technique is to reorganize all but a computationally insignificant part of this calculation into an equivalent problem: compute the displacements due to a single circular loading element at a total of m n stations (where m n is the product m × n). We solve this “parallel” problem at high computational speed by utilizing the sparse evaluation and massive interpolation (SEMI) method. Because the product m n that arises in our parallel problem is normally very large, we take maximum possible advantage of the acceleration achieved by the SEMI algorithm.     

Validation of methods to estimate design discharge flow rates for dam spillways with large regulating capacity

Abstract

The estimation and review of discharge flow rates in hydraulic works is a fundamental problem in water management. In the case of dams with large regulating capacity, in order to estimate return periods of discharge flow rates from the spillways, it becomes necessary to consider both peak flow and volume of the incoming floods. In this paper, the results of the validation for several methods of assessing design floods for spillways of dams with a large flood control capacity are presented; the validation is performed by comparing the maximum outflows (or the maximum levels reached in the reservoir) obtained from the routing of the design floods with those obtained from the routing of the historical annual maximum floods. The basin of Malpaso Dam, Mexico, is used as the case study.

Editor D. Koutsoyiannis

Citation Domínguez, M.R. and Arganis, J.M.L., 2012. Validation of methods to estimate design discharge flow rates for dam spillways with large regulating capacity. Hydrological Sciences Journal, 57 (3), 460–478.     

A practical numerical method for large strain liquefaction analysis of saturated soils

Accurate prediction of the liquefaction of saturated soils is based on strong coupling between the pore fluid phase and soil skeleton. A practical numerical method for large strain dynamic analysis of saturated soils is presented. The u–p formulation is used for the governing equations that describe the coupled problem in terms of soil skeleton displacement and excess pore pressure. A mixed finite element and finite difference scheme related to large strain analysis of saturated soils based on the updated Lagrangian method is given. The equilibrium equation of fluid-saturated soils is spatially discretized by the finite element method, whereas terms associated with excess pore pressure in the continuity equation are spatially discretized by the finite difference method. An effective cyclic elasto-plastic constitutive model is adopted to simulate the non-linear behavior of saturated soils under dynamic loading. Several numerical examples that include a saturated soil column and caisson-type quay wall are presented to verify the accuracy of the method and its usefulness and applicability to solutions of large strain liquefaction analysis of saturated soils in practical problems.     

Reflection and transmission coefficients for transition layers

Summary Formulae are derived for the reflection and transmission coeficients of plane elastic waves for a transition layer. Haskell's technique and the so-called delta matrices[5, 7] are used for this purpose. No problems are encountered in deriving the reflections and transmission coefficients from Haskell's matrices[3]. However, in some cases Haskell's matrices do not guarantee the accuracy required. For this reason attention is mainly devoted to deriving the reflection and transmission coefficients from the delta matrices. In deriving the transmission coefficients use is made of the fact that some3×3 subdeterminants of the delta matrices are squares of the3×3 subdeterminants of Haskell's matrices.     

Asymptotic treatment of the stability of a rotating layer of fluid with rigid boundaries

Abstract

Chandrasekhar (1961) has summarized the stability results of Bénard convection in a rotating fluid for the cases where the boundary surfaces are both rigid and free, and for both exchange of stabilities and overstability. His analysis provides very accurate results for a limited range of Taylor number J. Bisshopp and Niiler (1965) presented an asymptotic analysis of the rigid boundary problem for exchange of stabilities which is valid for very large Taylor number. The present paper makes use of modern rotating fluid theory to develop an approximate scheme for evaluating the Rayleigh number and other parameters and variables. Known asymptotic results for the free boundary problem at large J are used and an expansion in powers of E^1/6 (the Ekman number, E = 2J ^?½) yields a sequence of equations and appropriate boundary conditions for the rigid boundary problem. After the algorithm for the calculation is developed, results are given for the problem to second order in the expansion parameter for the case of exchange of stabilities and to first order in the expansion parameters for the overstable case. Ekman boundary layers are important in the development as one might expect. However, an additional, diffusive boundary layer of thickness E^? is necessary to provide the details of the temperature field. This boundary layer is the thermal response in the vertical direction to the horizontal spacing of the cells which is also order E^?. The horizontal spacing of the cells is essentially a series of contiguous, Stewartson (1957) layers of thickness E^?.     

An explicit method for numerical simulation of wave equations: 3D wave motion

In this paper, an explicit method is generalized from 1D and 2D models to a 3D model for numerical simulation of wave motion, and the corresponding recursion formulas are developed for 3D irregular grids. For uniform cubic grids, the approach used to establish stable formulas with 2M-order accuracy is discussed in detail, with M being a positive integer, and is illustrated by establishing second order (M=1) recursion formulas. The theoretical results presented in this paper are demonstrated through numerical testing.     

Double parameterized regularization inversion method for migration velocity analysis in transversely isotropic media with a vertical symmetry axis

Simultaneous estimation of velocity gradients and anisotropic parameters from seismic reflection data is one of the main challenges in transversely isotropic media with a vertical symmetry axis migration velocity analysis. In migration velocity analysis, we usually construct the objective function using the l₂ norm along with a linear conjugate gradient scheme to solve the inversion problem. Nevertheless, for seismic data this inversion scheme is not stable and may not converge in finite time. In order to ensure the uniform convergence of parameter inversion and improve the efficiency of migration velocity analysis, this paper develops a double parameterized regularization model and gives the corresponding algorithms. The model is based on the combination of the l₂ norm and the non‐smooth l₁ norm. For solving such an inversion problem, the quasi‐Newton method is utilized to make the iterative process stable, which can ensure the positive definiteness of the Hessian matrix. Numerical simulation indicates that this method allows fast convergence to the true model and simultaneously generates inversion results with a higher accuracy. Therefore, our proposed method is very promising for practical migration velocity analysis in anisotropic media.     

A deep-learning-based approach for seismic surface-wave dispersion inversion (SfNet) with application to the Chinese mainland

Surface-wave tomography is an important and widely used method for imaging the crust and upper mantle velocity structure of the Earth. In this study, we proposed a deep learning (DL) method based on convolutional neural network (CNN), named SfNet, to derive the v_S model from the Rayleigh wave phase and group velocity dispersion curves. Training a network model usually requires large amount of training datasets, which is labor-intensive and expensive to acquire. Here we relied on synthetics generated automatically from various spline-based v_S models instead of directly using the existing v_S models of an area to build the training dataset, which enhances the generalization of the DL method. In addition, we used a random sampling strategy of the dispersion periods in the training dataset, which alleviates the problem that the real data used must be sampled strictly according to the periods of training dataset. Tests using synthetic data demonstrate that the proposed method is much faster, and the results for the v_S model are more accurate and robust than those of conventional methods. We applied our method to a dataset for the Chinese mainland and obtained a new reference velocity model of the Chinese continent (ChinaVs-DL1.0), which has smaller dispersion misfits than those from the traditional method. The high accuracy and efficiency of our DL approach makes it an important method for v_S model inversions from large amounts of surface-wave dispersion data.     

Determination of focal coordinates of near earthquakes in an area of unknown crustal parameters

Summary A new computation method is described which determines the focal coordinates, theP _n-andP _G-velocity, the focal time and the Moho-depth in the case of near earthquakes. The method is provided especially for areas where the crustal parameters are insufficiently known. It presumes that the standard error ofP _n-andP _G-arrival times is a minimum. The problem turns out to represent the solution of a set of at least seven equations with seven unknown parameters. Four of these unknowns can be separated in advance. This leads to a considerable simplification. The three remaining unknowns are the focal coordinates. Their best fitting values are found by the gradient method. An example is presented which shows how errors can arise if the initial trial coordinate is chosen at an unfavourable position.     

FUNDAMENTAL FUNCTIONS FOR HORIZONTALLY STRATIFIED EARTH*

The potential distribution and the wave propagation in a horizontally stratified earth is considered and the analogy of the mathematical expression for seismic transfer function, electromagnetic and electric kernel functions, and magnetotelluric input impedance is discussed. Although these specific functions are conveniently treated by a separate expression in each method, it is indicated that the function for seismic and electromagnetic methods is mathematically the same with a change in the physical meaning of the variables from one method to the other. Similarly, the identity of the mathematical expressions of the resistivity kernel function and magnetotelluric input impedance is noticed. In each method a specific geophysical function depends on the thickness and the physical properties of the various layers. Every specific function involves two interdependent fundamental functions, that is P_n and Q_n, or P_n and P*_n, having different physical meaning for different methods. Specific functions are expressible as a ratio P_n/Q_n or P*_n/P_n. Fundamental functions may be reduced to polynomials. The fundamental polynomials Q*_n and P*_n describing the horizontally stratified media are a system of polynomials orthogonal on the unit circle, of first and second order, respectively. The interpretation of geophysical problems corresponds to the identification of the parameters of a system of fundamental orthogonal polynomials. The theorems of orthogonal polynomials are applied to the solution of identification problems. A formula for calculating theoretical curves and direct resistivity interpretation is proposed for the case of arbitrary resistivity of the substratum. The basic equation for synthetic seismograms is reformulated in appendix A. In appendix B a method is indicated for the conversion of the seismic transfer function from arbitrary to perfectly reflective substratum.     

Fragility curves for reinforced concrete structures in Skopje (Macedonia) region

A method for the development of earthquake intensity–damage relations, given as fragility curves and damage probability matrices is proposed in this paper. The proposed method is applied on reinforced-concrete frame-wall structures. Two sets of fragility curves and damage probability matrices are developed. The first one is for reinforced-concrete frame structures lower than 10 stories. For this purpose, a six-story frame structure is used. The other set is defined for reinforced-concrete frame-wall structures higher than 10 stories. A 16-story frame-wall structure was chosen as a sample. The sample structures were designed according to Macedonian design code. The conditions of the local seismic hazard were the subject of special concern for the development of earthquake intensity–damage relations. Because of the limited number of real time histories from the Skopje region, a set of 240 synthetic time histories were generated. Geological dates from the Skopje region were used. Response of the sample structures under earthquake excitation was defined performing nonlinear dynamic analysis. Modeling of the nonlinear behavior of the structural elements was completed according to state-of-the-art methods in this field. A modified Park and Ang damage model was chosen as a measure of the structure's response to earthquake excitation. Five damage states were defined to express the condition of damage. As a result of the analytical research, the values of the global damage index corresponding to each damage state were determined. Using the dates from the nonlinear dynamic analysis of the sample structures under all 240 synthetic time histories, the two sets of fragility curves and damage probability matrices were defined.     

NON-NORMAL INCIDENCE INVERSION: EXISTENCE OF SOLUTION*

An inverse problem is one in which the parameters of a model are determined from measured seismic data. Important to the solution of inverse problems is the issue of whether or not a solution exists. In this paper we show, in a constructive manner, that a solution does exist to the specific inverse problem of determining the parameters of a horizontally stratified, lossless, isotropic and homogeneous layered system that is excited by a non-normal incidence (NNI) plane wave. Mode conversion between P- and S-waves is included. We develop a seven-step layer-recursive procedure for determining all of the parameters for layer j. These parameters are P-wave and S-wave velocities and angles of incidence, density, thickness, traveltimes, and reflection- and transmission-coefficient matrices. Downward continuation of data from the top of one layer to the top of the next lower layer is an important step in our procedure, just as it is in normal incidence (NI) inversion. We show that, in order to compute all parameters of layer j, we need to (and can) compute some parameters for layer j+ 1. This is a non-causal phenomenon that seems to be necessary in NNI inversion but is not present in NI inversion.