Research overview 2002

The MAGMA center was created, supported by the European Commission It is expected to significantly affect the development of the whole department in the next three years. The year 2002 was unique also in a very large number of publications (see the list below) and awards: Professor V. Červený, whose 70th anniversary was celebrated in 2002, received the Gold Commemorative Medal of Charles University, the Professional Medal in Physics of the First Degree from the Association of Czech Mathematicians and Physicists, the Commemorative Medal of the Slovak Academy of Sciences, and the Honorary Membership of EAGE. Medals of the Slovak Academy of Sciences were awarded also to O. Novotný and J. Zahradník. K. Žáček received the Award of Merit from SEG for the best student paper at the 2001 SEG annual meeting. O. Šrámek obtained the French-Czech PhD scholarship.

Similarly to 2000 and 2001, research at the Department of Geophysics was carried out in three directions:

  1. Geodynamics
  2. Theory of seismic waves
  3. Earthquake studies


(reported by O. Čadek, Z. Martinec and C. Matyska)

Dynamic geoid

The investigation of the non-hydrostatic geoid concentrated on the role of lateral viscosity variations in the boundary layers, namely in the top 300 km of the mantle and in the core-mantle boundary region. As in the last years, the research was carried out in close cooperation with Luce Fleitout (ENS, Paris). The main results obtained for the upper mantle are summarized in the paper by Čadek and Fleitout (GJI, to appear in March 2003). For the lowermost part of the mantle, Čadek and Fleitout attempted to answer the question whether the free-slip boundary condition, commonly imposed in geodynamic modeling at the core-mantle boundary, really corresponds to the physical situation in the lowermost part of the mantle. Solving the inverse problem for laterally dependent viscosity in D", they found large-scale regions of a very high viscosity at the CMB that well correlate with the surface distribution of the hotspots. This result can be interpreted as indication of a thin layer of denser-then-average material of crustal origin located at the top of the core. This material, apparently associated with high-viscosity regions but in reality corresponding to a no-slip boundary condition, contains higher-than-average amount of radioactive elements which gives rise to mantle plumes and hotspot activity.

Postglacial rebound

Hanyk et al. (2002) published the method providing the whole Laplacian spectrum of postglacial viscoelastic relaxation of spherical Earth models. The first step toward incorporating the postseismic deformation modeling into the initial-value approach was performed (Matyska, 2002). The code for glacial isostatic modelling was supplemented by the possibility to compute the mass redistribution of ocean water and ice (the MSc Thesis of Ondřej Šrámek). The problem is governed by the sea-level equation which has been implemented in both spectral and spatial forms. A benchmark of the codes for implementing the sea-level equation was organized by O. Čadek and Z. Martinec; Jan Hagedorn, a PhD. student from GeoForschungsZentrum in Potsdam participated, too. The benchmark was successfully completed with an excellent agreement between two different approaches of implementing sea-level equation. It was also demonstrated that the present-day geoid can be significantly (up to 25%) influenced by the postglacial rebound, in particular in the case of highly viscous lower mantle. One of the most important activities of the group dealing with the postglacial rebound was also another benchmark, involving the codes for modeling the viscoelastic relaxation of the Earth with laterally dependent viscosity (Martinec, 2002a). The benchmark was organized in close cooperation with Luce Fleitout (ENS, Paris), Giorgio Spada (University Urbino), Carlo Giunchi and Spina Gianetti (ING, Rome). The goal of the benchmark was to test the spectral codes developed in Prague, and the finite-element codes used in other laboratories, and to built a standard that could be used by other postglacial rebound modelers to test their own 2D and 3D codes. The benchmark is still in progress and it will be finished during 2003.

Subduction of oceanic lithosphere

Subduction of oceanic lithosphere was studied both in elastic and viscous approximations. The elastic solution concentrated on the effects of pre-stress, phase transitions and thermal stress (Slancová and Čadek, submitted). The research of viscous slabs focused on the subduction initiation, as well as on the interaction between the subducting plate, the continental plate, and the 660 km discontinuity (Čížková et al., 2002). Kukačka and Matyska (2002a,b) showed that a strong lithospheric weakening is necessary to obtain realistic behavior of subducting lithosphere.

Thermal convection

To fill the gap between the tomographic models and computer simulations of thermal convection, Matyska and Yuen (2002) compared magnitudes of the 3-D local Bullen parameter obtained from thermal convection modeling with those provided by the seismic tomography. The agreement between both estimates points to a potential importance of the Bullen parameter in further studies of mantle dynamics. Mistr et al. (2002) dealt with applications of finite differences to thermal convection modeling in rotating systems with finite Prandtl number. Stability of convecting layers with different densities was studied in (Čížková and Matyska, 2002). They demonstrated that the density jump of only 3% is sufficient to keep the stability of the interface between the layers.

Electromagnetic induction

Forward modeling of electromagnetic induction response of the Earth to transient ionospheric- and magnetospheric-current excitation was carried out. New satellite data require the electromagnetic induction modeling to be performed directly in time domain. We managed to construct a semi-implicit time integration scheme for magnetic-diffusion equation and showed that this scheme is numerically stable (Velímský et al., submitted). Our cooperation with Mark Everett from Texas A&M University was further extended on modeling of satellite magnetic induction data (Everett and Martinec, in press). The cooperation resulted in the invitation of Jakub Velímský to stay at Texas A&M University as a post-doc for the next two years.

Gradiometric data modeling

The research was completed on finding an adequate mathematical tool for inverting satellite gradiometric data into information on the external gravitational potential of the Earth (Martinec, in press). These data will be available since 2006 when the GOCE gradiometric mission will be launched. We managed to solve gradiometric boundary-value problem in terms of Green's functions that were expressed in spectral form as series of tensor spherical harmonics. This form of the solution can be applied to develop the gravitational field in terms of spherical harmonics from the GOCE data. Alternatively, by means of the addition theorems for spherical harmonics and the infinite-sum formulae for Legendre polynomials, the spectral forms were converted to closed spatial forms. These forms can be used to construct the upward- and downward-continuation operators that can be applied to transform the measured gradiometric data from a non-spherical satellite orbit to a mean-orbit sphere.

Theory of seismic waves

reported by L. Klimeš)

Derivatives and perturbations of travel time

The equations for calculation of the third-order and higher-order partial derivatives of travel time (Klimeš submitted-b). The equations for calculation of the second-order and higher-order perturbations of travel time in both isotropic and anisotropic media (Klimeš 2002f) have been derived. The first-order perturbations of both real-valued and complex-valued rays and travel times have also been studied (Červený 2002b; Červený & Pšenčík 2002b).

Coupling ray theory

The equations for the second-order perturbations of travel time have been applied to the estimation of the errors due to the common-ray approximations of the coupling ray theory (Klimeš & Bulant 2002, submitted; Bulant & Klimeš 2002d). The errors due to the common-ray approximations of the coupling ray theory and the errors due to other quasi-isotropic approximations of the coupling ray theory have been demonstrated on numerical examples (Bulant & Klimeš 2002a, 2002d; Klimeš & Bulant 2002, submitted; Bucha & Bulant 2002). Various kinds of the coupling ray theory for weakly anisotropic models have been studied and compared with the exact solution derived for the "twisted crystal" model (Bulant & Klimeš 2002c, 2002e, submitted).

Ray methods and paraxial ray approximation

Fermat's variational principle for anisotropic inhomogeneous media and its relation to the Finslerian metric and to the Hamiltonian have been studied (Červený 2002a). The expressions for the curvatures of the slowness and ray-velocity surfaces in terms of the corresponding wave-propagation Finslerian metric tensor have been found (Klimeš 2002g, 2002k). The dependence of the paraxial-ray propagator matrix on the Hamiltonian have been demonstrated (Klimeš 2002j). The equations for the linear paraxial approximation of the polarization vectors and for the variation of the polarization vectors with a velocity perturbation have been published (Klimeš 2002a). The formulae for the ray-theory amplitudes in laterally varying layered isotropic 2-D models have been specified in detail (Červený & Pšenčík 2002a).

Reflection/transmission coefficients

Explicit equations for approximate linearized reflection/transmission coefficients at a generally oriented weak-contrast interface separating generally anisotropic media have been derived (Klimeš 2002h, 2002i, submitted-a). The study of the reflection/transmission coefficients in viscoelastic media continued (Brokešová & Červený 2002a, 2002b).

Lyapunov exponents, Sobolev scalar products, and model smoothing

The equations derived for the estimation of the average Lyapunov exponents (Klimeš 2002b), describing the ray chaos due to heterogeneities in the velocity model, have been applied to the construction of velocity models suitable for ray tracing and other high-frequency asymptotic methods (Bulant 2002a; Bulant & Klimeš 2002b). The designed algorithm of constructing velocity models have been tested on various 2-D and 3-D synthetic structures (Bulant 2002b, 2002c, submitted; Žáček 2002a; Bucha & Bulant 2002).

Inversion of seismic data

Particular attention have been devoted to the resolution of seismic inversion techniques, model fitting, conversion and smoothing with the application of medium correlation functions and Sobolev scalar products (Klimeš 2002c, 2002d, 2002e; Bulant & Klimeš 2002b).

Gaussian-packet prestack depth migration

A new algorithm of the optimization of the shape of Gaussian beams and packets for prestack migrations in complex models have been derived, including smoothing of the beam parameters over the Hamiltonian hypersurface in the phase space (Žáček 2002b). The theoretical and numerical study of the decomposition of the time sections into optimized Gaussian packets, which is of key importance in the Gaussian packet migration, has been commenced (Žáček 2002c).

Earthquake studies

reported by J. Zahradník)

Seismic stations

Seismic stations of the Charles University, jointly operated with the University of Patras since 1997 in Greece were upgraded. The present network comprises three stations in the western part of the Corinth Gulf, each one equipped with a weak-motion broad-band velocigraph CMG3-T and a strong-motion accelerograph CMG5-T. The stations are Sergoula, Mamousia, and the University of Patras. Under testing is a fourth station (CMG5-T) in the eastern part of the Corinth Gulf at Loutraki. The continuous and triggered three component records are sampled at 20 Hz, and 100 Hz, respectively. The selected data are available from, updated every 4 months.

3D earthquake modeling

3D earthquake modeling based on a hybrid combination of the source, path and site effects was published as a methodical paper (Opršal a Zahradník, 2002). It allows fast calculations up to frequencies of engineering interest (10-20 Hz) on standard personal computers. Not only hybrid coupling of the 3D finite differences for local site effects with discrete wavenumbers for 1D crustal propagation was tested (Opršal and Zahradník, 2002), but also the coupling of the 3D finite difference method with the ray method. The latter will allow inclusion of the lateral heterogeneities of the crust (Opršal et al., 2002).

Strong motion simulation methods

Three alternative methods have been developed for numerical modeling of the strong ground motions due to finite-extent sources, using combined deterministic and stochastic approaches: (a) The composite-source method with equal sized subevents, in which higher frequencies are obtained by stochastic perturbation and extrapolation of the deterministic low-frequency solution, the co-called PEXT method (Zahradník a Tselentis, 2002). (b) The composite model based on fractal distribution of subevent sizes, with subevents treated as point sources or finite sources (MSc Thesis of Jan Burjánek), (c) The kinematic method of the k**(-2) stochastic slip distribution, with k-dependent rise time, with or without asperities (MSc Thesis of František Gallovič); see also two submitted papers (Gallovič and Brokešová). The three methods were compared with each other, as well as with the 1999 Athens data (accelerograms and empirical attenuation relations). The main difference between the methods is in the amount of their source directivity in the high frequency range, a lower directivity in the composite models, and a higher directivity in the kinematic models. Since the data do not simply prefer one of them, we suggest the ground motion predictions to be performed with all the three methods, simultaneously.

Empirical Green's function method

Modification of the empirical Green's function method, allowing unequal focal mechanisms of the mainshock and aftershock was finished and published with an application to the 1999 Athens earthquake (Plicka a Zahradník, 2002).

The Mw=6.5 Skyros earthquake

The Mw=6.5 Skyros earthquake, Aegean Sea, was modeled using regional broad-band stations in Greece. Inversion of the amplitude spectra of complete waveforms (0.05-0.08 Hz) was used to retrieve the focal mechanism and source depth. Using the aftershock distribution, the earthquake was interpreted as a left-lateral strike slip, not very common in the studied area. The fault length was estimated by combining observed mainshock spectra and synthetic spectra of a weak event, representing impulse response of the medium. Similar results were obtained by means of a true aftershock. The waveform modeling (0.05-0.20 Hz) was performed for the 20 x 10 km finite-extent fault. It showed that the rupture propagation was predominantly unilateral, from NW to SE (Zahradník, 2002).

The Egion Mw=4.3 earthquake

The Egion Mw=4.3 earthquake, Corinth Gulf An integrated study of the earthquake sequence of February-July 2001 in the western part of the Corinth Gulf was performed. Several location methods were applied using the regional network PATNET (Janský et al., submitted). The preferred grid search location suggested that the main cluster of the sequence and the mainshock had a depth of 16 km. The amplitude spectra of complete waveforms of the mainshock at three local stations were inverted between 0.1 and 0.2 Hz for double-couple focal mechanism and depth, validated by forward waveform modeling. The optimum solution provided the depth of 8 km. To resolve the remarkable difference between the two depth estimates, the local short-period network of the EC Corinth Rift Laboratory was additionally employed, which clearly confirmed the 8 km depth. The T axis of this earthquake, is consistent with the regional direction of extension ~N10°. However, none of the nodal planes can be associated to an active structure seen at the surface. The relationship of this earthquake with deeper faults remains open (Zahradnik et al., submitted).