Modelling of ground deformation and gravity fields using finite element method: an application to Etna volcano

Elastic finite element models are applied to investigate the effects of topography and medium heterogeneities on the surface deformation and the gravity field produced by volcanic pressure sources. Changes in the gravity field cannot be interpreted only in terms of gain of mass disregarding the ground deformation of the rocks surrounding the source. Contributions to gravity changes depend also on surface and subsurface mass redistribution driven by dilation of the volcanic source. Both ground deformation and gravity changes were firstly evaluated by solving a coupled axisymmetric problem to estimate the effects of topography and medium heterogeneities. Numerical results show significant discrepancies in the ground deformation and gravity field compared to those predicted by analytical solutions, which disregard topography, elastic heterogeneities and density subsurface structures. With this in mind, we reviewed the expected gravity changes accompanying the 1993–1997 inflation phase on Mt Etna by setting up a fully 3-D finite element model in which we used the real topography, to include the geometry, and seismic tomography, to infer the crustal heterogeneities. The inflation phase was clearly detected by different geodetic techniques (EDM, GPS, SAR and levelling data) that showed a uniform expansion of the overall volcano edifice. When the gravity data are integrated with ground deformation data and a coupled FEM modelling was solved, a mass intrusion could have occurred at depth to justify both ground deformation and gravity observations.     

Effects of mechanical layering on volcano deformation

The migration and accumulation of magma beneath volcanoes often causes surface displacements that can be measured by geodetic techniques. Usually, deformation signals are explained using models with uniform mechanical properties. In this paper, we study surface displacements due to magma chamber inflation, using heterogeneous finite element models. We first present a systematic analysis of the influence of mechanical layering, showing that the stiffness contrast significantly affects the entity and the pattern of vertical and radial displacements. Second, as an example we apply the models to interpret ground displacements at Darwin volcano (Galápagos Islands) as revealed by InSAR data in the period 1992–1998. The considered models suggest that geodetic data interpreted using homogeneous models leads to underestimation of the source depth and volume change. Thus, we propose correction factors for the source parameters estimated by homogeneous models, in order to consider a range of variation due to mechanical layering as analysed in this study. The effect of the mechanical heterogeneities affects the correct understanding of geodetic data and also influences the evaluation of a volcanic hazard potential.     

Ground deformation in a viscoelastic medium composed of a layer overlying a half-space: a comparison between point and extended sources

We obtain and compare analytical and numerical solutions for ground displacement caused by an overpressurized magma chamber placed in a linear viscoelastic medium composed of a layer over a half-space. Different parameters such as size, depth and shape of the chamber, crustal rheology and topography are considered and discussed. Numerical solutions for an axisymmetric extended source are computed using a finite element method (FEM). Analytical solutions for a point source are obtained using the dislocation theory and the propagator matrix technique. In both cases, the elastic solutions are used together with the correspondence principle of linear viscoelasticity to obtain the solution in the Laplace transform domain. Viscoelastic solutions in the time domain are derived inverting the Laplace transform using the Prony series method. The differences between the results allow us to constrain the applicability of the point source and the flat surface hypothesis, which are usually implicitly assumed when analytical solutions are derived. The effect of the topography is also considered. The results obtained show that neglecting the topographic effects may, in some cases, introduce an error greater than that implicit in the point-source hypothesis. Therefore, for an adequate modelling and interpretation of the time-dependent displacements, topography must be considered.     

Inverse scattering of surface waves: imaging of near-surface heterogeneities

An efficient inverse scattering method is developed for imaging near-surface heterogeneities using scattered surface waves. Three dimensional elastodynamic wave propagation and scattering in a laterally invariant embedding medium is considered. The Born Approximation is used and the scattered wavefield is expressed as a domain type integral representation. The computation time of Green's tensor elements is reduced by considering the radial symmetry of the medium. The method is validated by numerical tests. Ultrasonic laboratory data obtained from a scale model experiment are used for imaging the near-surface inhomogeneities caused by an epoxy-filled hole in the surface of an aluminum block. Both synthetic and the scale model tests show that the location, the actual density contrast and the depth of the inhomogeneities are reasonably well estimated.     

Effects of topography and crustal heterogeneities on the source estimation of LP event at Kilauea volcano

The main goal of this study is to improve the modelling of the source mechanism associated with the generation of long period (LP) signals in volcanic areas. Our intent is to evaluate the effects that detailed structural features of the volcanic models play in the generation of LP signal and the consequent retrieval of LP source characteristics. In particular, effects associated with the presence of topography and crustal heterogeneities are here studied in detail. We focus our study on a LP event observed at Kilauea volcano, Hawaii, in 2001 May. A detailed analysis of this event and its source modelling is accompanied by a set of synthetic tests, which aim to evaluate the effects of topography and the presence of low velocity shallow layers in the source region. The forward problem of Green's function generation is solved numerically following a pseudo-spectral approach, assuming different 3-D models. The inversion is done in the frequency domain and the resulting source mechanism is represented by the sum of two time-dependent terms: a full moment tensor and a single force. Synthetic tests show how characteristic velocity structures, associated with shallow sources, may be partially responsible for the generation of the observed long-lasting ringing waveforms. When applying the inversion technique to Kilauea LP data set, inversions carried out for different crustal models led to very similar source geometries, indicating a subhorizontal cracks. On the other hand, the source time function and its duration are significantly different for different models. These results support the indication of a strong influence of crustal layering on the generation of the LP signal, while the assumption of homogeneous velocity model may bring to misleading results.     

Small-scale mantle flows and induced lithospheric stress near island arcs

Summary. This paper explores the middle ground between complex thermally-coupled viscous flow models and simple corner flow models of island arc environments. The calculation retains the density-driven nature of convection and relaxes the geometrical constraints of corner flow, yet still provides semianalytical solutions for velocity and stress. A novel aspect of the procedure is its allowance for a coupled elastic lithosphere on top of a Newtonian viscous mantle. Initially, simple box-like density drivers illustrate how vertical and horizontal forces are transmitted through the mantle and how the lithosphere responds by trench formation. The flexural strength of the lithosphere spatially broadens the surface topography and gravity anomalies relative to the functional form of the vertical flow stresses applied to the plate base. I find that drivers in the form of inclined subducting slabs cannot induce self-driven parallel flow; however, the necessary flow can be provided by supplying a basal drag of 1–5 MPa to the mantle from the oceanic lithosphere. These basal drag forces create regional lithospheric stress and they should be quantifiable through seismic observations of the neutral surface. The existence of a shallow elevated phase transition is suggested in two slab models of 300 km length where a maximum excess density of 0.2 g cm⁻³ was needed to generate an acceptable mantle flow. A North New Hebrides subduction model which satisfies flow requirements and reproduces general features of topography and gravity contains a high shear stress zone (75 MPa) around the upper slab surface to a depth of 150 km and a deviatoric tensional stress in the back arc to a depth of 70 km. The lithospheric stress state of this model suggests that slab detachment is possible through whole plate fracture.     

Born scattering of longperiod body waves

The Born approximation is applied to the modelling of the propagation of deeply turning longperiod body waves through heterogeneities in the lowermost mantle. We use an exact Green's function for a spherically symmetric earth model that also satisfies the appropriate boundary conditions at internal boundaries and the surface of the earth. The scattered displacement field is obtained by a numerical quadrature of the product of the Green's function, the exciting wavefield and structural perturbations. We study three examples: scattering of longperiod P waves from a plume rising from the coremantle boundary (CMB), generation of longperiod precursors to PKIKP by strong, localized scatterers at the CMB, and propagation of corediffracted P waves through largescale heterogeneities in D". The main results are as follows: (1) the signals scattered from a realistic plume are small with relative amplitudes of less than 2 per cent at a period of 20 s, rendering plume detection a fairly difficult task; (2) strong heterogeneities at the CMB of appropriate size may produce observable longperiod precursors to PKIKP in spite of the presence of a diffraction from the PKP B caustic; (3) corediffracted P  waves ( P _diff) are sensitive to structure in D" far off the geometrical ray path and also far beyond the entry and exit points of the ray into and out of D"; sensitivity kernels exhibit ringshaped patterns of alternating sign reminiscent of Fresnel zones; (4) P _diff also shows a nonnegligible sensitivity to shear wave velocity in D"; (5) down to periods of 40 s, the Born approximation is sufficiently accurate to allow waveform modelling of P _diff through largescale heterogeneities in D" of up to 5 per cent.     

Thin elastic shells with variable thickness for lithospheric flexure of one-plate planets

Planetary topography can either be modelled as a load supported by the lithosphere, or as a dynamic effect due to lithospheric flexure caused by mantle convection. In both cases the response of the lithosphere to external forces can be calculated with the theory of thin elastic plates or shells. On one-plate planets the spherical geometry of the lithospheric shell plays an important role in the flexure mechanism. So far the equations governing the deformations and stresses of a spherical shell have only been derived under the assumption of a shell of constant thickness. However, local studies of gravity and topography data suggest large variations in the thickness of the lithosphere. In this paper, we obtain the scalar flexure equations governing the deformations of a thin spherical shell with variable thickness or variable Young's modulus. The resulting equations can be solved in succession, except for a system of two simultaneous equations, the solutions of which are the transverse deflection and an associated stress function. In order to include bottom loading generated by mantle convection, we extend the method of stress functions to include loads with a toroidal tangential component. We further show that toroidal tangential displacement always occurs if the shell thickness varies, even in the absence of toroidal loads. We finally prove that the degree-one harmonic components of the transverse deflection and of the toroidal tangential displacement are independent of the elastic properties of the shell and are associated with translational and rotational freedom. While being constrained by the static assumption, degree-one loads can deform the shell and generate stresses. The flexure equations for a shell of variable thickness are useful not only for the prediction of the gravity signal in local admittance studies, but also for the construction of stress maps in tectonic analysis.     

Local and regional components of western Aegean deformation extracted from 100 years of geodetic displacement measurements

Our objectives are as follows. First, we wish to develop a methodology to recover the long-term component of deformation from any set of distributed, time-averaged geodetic strain measurements that were subject to seismic disturbance, given a catalogue of local seismicity that occurred during the measurement period. Second, using seismic and geodetic data sets that span approximately 100 years, we apply this technique in the western Aegean to assess the role of local seismicity in regional deformation. The methodology is developed using a model for crustal deformation constructed from a long-term, smooth regional strain field combined with instantaneous, local perturbations from upper-crustal earthquakes approximated by static elastic dislocations. By inverting geodetic displacements for the smooth field while simultaneously floating influential but uncertain earthquake source parameters, an estimate of the regional component of deformation that is approximately independent of the seismicity can be made. In the western Aegean we find that the horizontal component of regional deformation can be described with minor inaccuracy by a quadratic relative displacement field. The principal horizontal extensional axes calculated from the regionally smooth displacement field agree in orientation with the T-axes of earthquakes in the region. These observations indicate that the instantaneous elastic strain of the 10 km thick seismogenic layer is driven by a stress field that is smooth on the scale of the geodetic network as a whole, 200-300 km.     

The onset of extension during lithospheric shortening:a two-dimensional thermomechanical model for lithospheric unrooting

We model the evolution of the lithosphere during its shortening and consequent gravitational collapse with special emphasis on the induced variations in the surface stress regime and dynamic topography. In particular, we analyse the conditions leading, immediately after lithospheric failure, to local extension, eventually coeval with compression. Different crustal rheologies and kinematic conditions as well as thermally imposed mechanical rupture are considered. Numerical calculations have been performed by using a 2-D finite element code that couples the thermal and mechanical equations for a Newtonian rheology with a temperature-dependent viscosity. The results show that, after the failure of a gravitationally unstable lithospheric root, the replacement of lithospheric mantle by warmer asthenospheric material induces a considerable variation in the dynamic topography and in the surface stress regime. The occurrence of local extension, its intensity and its spatial distribution depend mainly on whether convergence continues throughout the process or ceases after or before the lithospheric failure. Similarly, uplift/subsidence and topographic inversion are controlled by kinematic conditions and crustal rheology. Mechanical rupture produces drastic changes in the surface stress regime and dynamic topography but only for a short time period, after which the system tends to evolve like a continuous model.     

A hybrid solution for wave propagation problems in regular media with bounded irregular inclusions

Summary. This paper describes a method that enables one to calculate the effects of localized heterogeneities on the wavefield in an otherwise regular medium. It does so by connecting a finite element solution for a heterogeneous inclusion to any type of solution for the regular medium, e.g. a reflectivity solution for a layered medium or an analytical solution for a simple half-space. Once the Green's functions for the regular medium are determined, the method reduces to a coupled set of algebraic equations for the wavefield, with the incident field and/or body forces as known variables. An efficient numerical scheme is derived for the solution of these equations.     

A critical assessment of viscous models of trench topography and corner flow

Summary. We obtain stresses for Newtonian viscous flow in simple geometries (e.g. corner flow, bending flow) in order to study the effect of imposed velocity boundary conditions. Stress for a delta function velocity boundary condition decays as 1/ r ²; for a step function velocity, stress goes as 1/ r ; for a discontinuity in curvature, the stress singularity is logarithmic. For corner flow, which has a discontinuity of velocity at a certain point, the corresponding stress has a 1/ r singularity. However, for a more realistic circular-slab model, the stress singularity becomes logarithmic. Thus the stress distribution is very sensitive to the boundary conditions, and in evaluating the applicability of viscous models of trench topography it is essential to use realistic geometries.
Topography and seismicity data from northern Honshu, Japan, were used to construct a finite element model, with flow assumed constant speed and tangent to the top of the grid, for both Newtonian and non-Newtonian flow (power law 3 rheology). Normal stresses at the top of the grid are compared to the observed trench topography. There is poor agreement. Purely viscous models of subducting slabs with simple, geometrically consistent velocity boundary conditions do not predict normal stress patterns compatible with observed topography. Elasticity and plasticity appear to be important in determining trench topography.     

Extension of the Thomson-Haskell method to non-homogeneous spherical layers

Summary. Amplitude spectra of Rayleigh and Love waves in a layered non-gravitating spherical earth have been obtained using as a source, displacement and stress discontinuities. In each layer elastic parameters and density follow specified functions of radial distance and the solutions of the equations of motion are obtained in terms of exponential functions. The Thomson—Haskell method is extended to this case. The problem reduces to simple calculations as in a plane-layered medium. Numerical results of phase and group velocities up to periods of 300 s in various earth models when compared with earlier results (obtained by numerical integration) show that the present method can be used with sufficient accuracy. The differences in phase velocity, group velocity and amplitude (also surface ellipticity in the case of Rayleigh waves) between spherical- and flat-earth models have been investigated in the range 20–300–s period and expressed in polynomials in the period.     

Seismic radiation from dynamic coalescence, and the reconstruction of dynamic source parameters on a planar fault

The dynamic coalescence of two mode II cracks on a planar fault is simulated here using the elastodynamic boundary integral equation method. We focus on the complexity of the resultant slip rate and seismic radiation in the crack coalescence model (CCM) and on the reconstruction of a single crack model (SCM) that can reproduce the CCM waveforms from heterogeneous source parameters rather than coalescence. Simulation results reveal that localized higher slip rates are generated by coalescence as a result of stress interaction between the approaching crack tips. The synthesized seismic radiation exhibits a distinct coalescence phase that has striking similarities to stopping phases in the radiation and propagation properties. The corresponding SCM yields a singular increase in the stress drop distribution, which is accompanied by a sudden decrease in it across the point of coalescence in the CCM. This implies that the generation of high-frequency radiation is more efficient from coalescence than from stopping, although both phenomena exhibit the same strong  ω⁻²  -type displacement spectra.     

Inference of core-mantle boundary topography from ISC PcP and PKP traveltimes

We invert ISC PcP and PKP absolute and differential traveltimes in an attempt to infer the long-wavelength topography of the core-mantle boundary (CMB). The data selection and processing methods are described and evaluated. These travel-time data are very noisy and the geographic distribution of the data is highly non-uniform, inhibiting reliable inference of CMB topography. Spatial averaging enhances the coherent component of the residual variance (related to heterogeneity), however, the random component of the variance is much larger than the coherent component. We show that for PcP data the coherent signal due to mantle heterogeneity overshadows that arising from the CMB, and that the effects of mantle heterogeneity are mapped into our inferred CMB solutions. The PcP data are not correlated across the spatial averaging bins and seem to have a strong bias due to small-scale structure and/or noise. The non-uniform geographic sampling of the data plays a role in the mapping of mantle heterogeneity onto the CMB. Spatial patterns of CMB models inferred from different phases do not agree. Amplitudes of seismically inferred CMB undulations vary greatly. The sensitivity of inferred CMB models to the processing, spatial averaging procedure, and inversion techniques are investigated. Topographic amplitudes increase strongly with increasing input residual variance. The power spectrum of inferred topography indicates that there are unmodelled heterogeneities that must be described with spherical harmonics of degree 6 and higher. Based on this work, we conclude that reliable inference of long-wavelength CMB topography is not likely with the current ISC data set or with a spherical harmonic expansion truncated to degree and order 6.     

大变形有限元法分析边坡稳定性

边坡稳定性分析的重要任务是确定边坡潜在滑面的位置及相应的安全系数,鉴于坡体位移可在现场监测中直接得到,提出边坡潜在滑面确定的位移判据,采用大变形有限元首先计算出边坡的位移场,然后在开挖面的不同高度作塑性区内的水平位移参考线,依次连接参考线上的位移突变点即得到边坡失稳的潜在滑面位置,分析滑面上的应力以及岩体抗剪强度即得到边坡的整体安全系数。     

The relationship between seismic deformation and deep-seated gravitational movements during the 1997 Umbria–Marche (Central Italy) earthquakes

This paper re-evaluates the origin of some peculiar patterns of ground deformation in the Central Apennines, observed by space geodetic techniques during the two earthquakes of the Colfiorito seismic sequence on September 26th, 1997. The surface displacement field due to the fault dislocation, as modelled with the classic Okada elastic formulations, shows some areas with high residuals which cannot be attributed to non-simulated model complexities. The residuals were investigated using geomorphological analysis, recognising the geologic evidence of deep-seated gravitational slope deformations (DSGSD) of the block-slide type. The shape and direction of the co-seismic ground displacement observed in these areas are correlated with the expected pattern of movement produced by the reactivation of the identified DSGSD. At least a few centimetres of negative “Line of Sight” ground displacement was determined for the Costa Picchio, Mt. Pennino, and Mt. Prefoglio areas. A considerable horizontal component of movement in the Costa Picchio DSGSD is evident from a qualitative analysis of ascending and descending interferograms. The timing of the geodetic data indicates that the ground movement occurred during the seismic shaking, and that it did not progress appreciably during the following months. This work has verified the seismic triggering of DSGSD previously hypothesized by many researchers. A further implication is that in the assessment of DSGSD hazard seismic input needs to be considered as an important cause of accelerated deformation.     

吉林省德惠市土壤特性空间变异特征与格局

Based on regionalized variable theory, semivariograms of geo-statistics were used to research the spatial variation of soil properties quantitatively. The results showed that the semivariogram of soil organic matter is best described by spherical model, the best model for semivariograms of soil total N and available K is exponential models and that of available P belongs to linear with sill model. Those soil properties have different spatial correlations respectively, the lag of organic matter is the highest and that of available P is the lowest, the spatial correlation of N and available K belongs to moderate degree. Spatial heterogeneities are different too, the degree of organic matter and total N are higher, the degree of available K is in the next place and that of available P is the lowest. Influenced by the shape, topography and soil of the study area, all isotropies of available P are obvious in all directions while anisotropies of others are manifested. According to the analytical results, supported by GIS, Kriging and IDW methods are applied to describe and analyze the spatial distribution of soil properties. The results indicate that soil organic matter, total N and available K are distributed regularly from northeast to southwest, while available P is distributed randomly.