The MAGMA center, i.e. the Prague Center of Mathematical Geophysics, Meteorology, and their Applications, supported by the European Commission, has started in 2003 its first year of existence with at a quite exciting "rate". The geophysical group profited from short visits of 10 senior researchers, as well as from 4 stays of PhD students and post-doctoral researchers, each stay lasting more than 1 month (7 man-months together). Moreover, two international meetings were co-organized in our country (the geophysical one being described below). For more details, see http://geo.mff.cuni.cz/magma. Simply speaking, MAGMA helps to shape our Department into a truly European dimension. Indeed, in a close relation with MAGMA, we successfully finished our work package of the 5^th framework EC project PRESAP, oriented to practical aspects of earthquake-hazard prediction. We also participated in several 6^th framework EC project proposals focused on the Earth research, its environmental aspects, and, in particular, on natural hazards (NERIES, 3HAZ-Corinth), as well as the proposal improving mobility and continuing education of young scientists (the Marie Curie research training networks SPICE and AEGEAN-QUAKE). Intensive cooperation with major oil companies has been strengthened, too.

Similarly to the previous years, research at the Department of Geophysics was carried out in three directions: Geodynamics, Theory of seismic waves, and Earthquake and structural studies.

The 8th European Workshop on Numerical Modeling of Mantle Convection and Lithospheric Dynamics was organized by our department (O. Čadek), in cooperation with the Geophysical Institute, Academy of Sciences, and took part in the Castle of Hrubá Skála, Czech Republic, from September 13 to September 18, 2003. There were 96 participants from 12 countries (including overseas), and 37 from them were supported by the MAGMA project. Four main topics of the workshop were:

Recent low-orbit satellite missions Oersted and CHAMP offer an opportunity to improve models of electrical conductivity of the Earth's mantle. Forward modeling of electromagnetic induction response of the Earth to Dst signal was carried out by Velímský et al. (2003). In cooperation with Mark Everett from Texas A&M University, the time-domain, spectral finite-element approach to transient 2-D geomagnetic induction in a spherical heterogeneous Earth was created and published (Martinec et al., 2003). Another paper, describing the time-domain spherical harmonic finite-element approach to transient 3-D geomagnetic induction was submitted by Velímský and Martinec. Special attention was devoted to the response of the realistic Earth's models to simulated geomagnetic storms (Everett and Martinec, 2003). The effect of heterogeneous surface conductance under ionospheric Sq excitation on satellite-observed responses was also studied. Preliminary results were submitted by Velímský and Everett. The cooperation with Mark Everett resulted in the stay of Jakub Velímský at Texas A&M University as a Post-Doc.

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, 2003). 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. The ellipsoidal Stokes boundary-value problem was studied to compute the geoidal heights. The disturbing potential was split into a low-degree reference potential known from global geopotential models, and a higher-degree potential obtained by solving the ellipsoidal Stokes problem in the form of surface integrals (Ardestani and Martinec, 2003a). The far-zone contribution was added in the other paper (Ardestani and Martinec, 2003b). The approach to the solution of the Dirichlet and the Stokes boundary-value problems on the Earth's ellipsoid was commented by Grafarend and Martinec (2003). The theory of the Bouger anomaly in spherical geometry for determination of precise geoid was developed in a submitted paper by Vanicek et al. (2003).

Modeling of the non-hydrostatic geoid focused 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 (Čadek and Fleitout, 2003). The authors solved the inverse problem for laterally dependent viscosity in D", and found large-scale regions of a very high viscosity at the CMB that correlate well with the surface distribution of hotspots. The stability of layered thermal convection models with hypothetical chemical interface at a depth of 1000 km was studied by Čížková and Matyska (special issue of Phys. Earth Planet. Inter., in press). They found that a density increase of about 3% is sufficient to maintain the layered convection for a time of hundreds of million of years. Kukačka and Matyska (special issue of Phys. Earth Planet. Inter., in press) incorporated a rheological weakening into subduction modeling in a new self-consistent manner. The output of their non-linear models points to a crucial role of the rheological weakening to the dip angle and heat dissipation during subduction. Matyska (2003) created a new way of describing the seismic wave field generation and coseismic static deformation due to finite-extent source in depth-dependent elastic models. His approach is based on a direct decomposition of the displacement field, without employing the usual Green function formalism. (The latter is an important step towards linking together the geodynamic and seismic research at our department, one of major and most difficult goals of the MAGMA project.) Hanyk et al. (in preparation) studied energy jumping into the Earth's mantle generated by glacial forcing, and found that for abrupt changes of the forcing function the dissipation of heat below the lithosphere can reach a magnitude substantially higher than the radiogenic heating, namely under the presence of a low viscosity layer (asthenosphere). Košál (MSc thesis, 2003) studied the mathematical properties of the free oscillations of the Earth's models with maxwellian rheology, and demonstrated why these models are numerically ill-posed. Inovecký (MSc thesis, 2003) created a new weak formulation of equations describing the postglacial rebound, and showed the numerical applicability of finite elements to this formulation, for models with strong lateral variations of viscosity.

Several other papers were submitted in 2003: Employment of the knowledge of the full magnetic induction vector from satellite missions to get a model of mantle electrical conductivity with higher resolution than models obtained from surface data is described in Martinec and McCreadie (2003). Viscoelastic relaxation of the Earth for a glacial loading of Fennoscandian extent was studied in the framework of 2-D models of mantle viscosity, enabling to capture the influence of continental roots (Martinec and Wolf, 2003). Several spatial and temporal contributions resulting in changes of the geoid above Greenland, namely those due to ice loading in North America and those resulting from the Greenland Ice Sheet, were evaluated by Fleming et al. (2003).

Most of work reported in this section belongs to the consortium "Seismic Waves in Complex 3-D Structures" (coordinated by Prof. V. Červený). This is a framework of our long-lasting cooperation with the research departments of major oil companies, worldwide.

Seismic waves in viscoelastic anisotropic media

Considerable attention has been devoted to harmonic plane waves in viscoelastic anisotropic media (Červený 2003a, 2003b, 2003c, 2004; Červený & Pšenčík 2003a, 2003b).

Explicit equations for approximate linearized reflection/transmission coefficients at a generally oriented weak-contrast interface separating generally anisotropic media have been derived and published (Klimeš 2003a).

The equations for calculating the first-order and higher-order spatial and perturbation derivatives of amplitude in both isotropic and anisotropic media have been derived (Klimeš 2003d).

The caustic identification algorithm for isotropic media has been generalized to anisotropic media, and the rules for the phase shift of the anisotropic-ray-theory wavefield due to caustics have been derived (Klimeš 2004b).

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 (Bulant & Klimeš 2003c; Klimeš & Bulant 2003, 2004). 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š 2003a, 2003b, 2004; Klimeš & Bulant 2004). Isotropic ray theory, anisotropic ray theory and various kinds of the coupling ray theory for weakly anisotropic models have been studied and compared with the exact solution derived for the "simplified twisted crystal" and "oblique twisted crystal" models (Bulant & Klimeš 2003a, 2003b, 2004; Bulant, Klimeš, Pšenčík & Vavryčuk 2003; Klimeš 2004a). Equations for the numerical common S-wave ray tracing and for the corresponding dynamic ray tracing in a smooth elastic anisotropic medium have been proposed (Klimeš 2003e).

The equations derived for the estimation of the average Lyapunov exponents, 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 2004). The designed algorithm of constructing velocity models has been compared with the smoothing methods of other authors and tested on various 2-D and 3-D synthetic structures (Bucha, Bulant & Klimeš 2003a, 2003b, Bulant 2004). Capabilities of the ray tracing software developed at our department have been further extended (Bucha & Bulant 2003a, 2003b; Klimeš 2003e). Numerical ray tracing has been tested on various 3-D synthetic structures, including the smoothed SEG/EAGE Salt Model (Bucha 2003; Bucha & Bulant 2003a, 2003b).

Particular attention has been devoted to the sensitivity of seismic waves to the structure (Klimeš 2003b, 2003c; Žáček & Klimeš) and to the physical meaning of velocity macro models (Bucha, Bulant & Klimeš 2003a, 2003b).

The decomposition of the time sections into optimized Gaussian packets is of key importance in the Gaussian packet migration. The equations for the decomposition have been derived and the decomposition was numerically tested (Žáček 2003a, 2003b).

The PhD dissertation by Plicka (2003), devoted to modeling of finite-extent seismic sources by empirical Green's functions, was successfully defended. In MSc thesis by Kolínský (2003), devoted to the dispersion of seismic surface waves, various methods of the frequency-time analysis were discussed, and a new computer code for computing group velocities was presented. The code includes also retrieval of the filtered time series, corresponding to a selected branch of the dispersion curve. The method was applied to surface waves along some Eurasian paths recorded at the seismic station Praha.

The stations, jointly operated with the University of Patras since 1997, were further upgraded. The present network comprises four sites, each one equipped with a weak-motion broad-band velocigraph CMG3-T and a strong-motion accelerograph CMG5-T. The stations Sergoula and Mamousia (existing in the previous period, too) are situated on the northern and southern coast of the Corinth Gulf, in its western part. They will be used, besides other, in a new EC project 3HAZ-Corinth, coordinated by P. Bernard. The third (new) site is Loutraki, at the eastern edge of the Corinth Gulf. The fourth site (also new) is Pylos, close to Kalamata city, on the south-west of the Pelloponesos. Thanks to massive investment of the Seismological Laboratory, Patras University, The Loutraki and Pylos stations have got the satellite data transmission to Patras (Libra system, Nanometrics company). The selected data are available from http://seis30.karlov.mff.cuni.cz, updated every 4 months.

The 3D modeling based on a hybrid combination of the source, path and site effects, methodically developed in 2002, has been applied in practice. The focus of our attention was one of a few site in Athens (the Ano Liosia suburb), where the 1999 earthquake had the largest intensity, IX. The source and path effects were modeled by the composite-source model and the discrete wavenumber method (method PEXT of J. Zahradnik), while the local sedimentary basin was modeled by the 3D finite-difference method (I. Opršal). The hybrid modeling proved to be an efficient tool up to frequency of engineering interest. As such, it enabled explanation of the damaging ground motions at the Ano Liosoa site as a combined effect of the source directivity and lateral heterogeneity of the site. The true ground motion during mainshock was not recorded at the site, but the calculations indicate that the horizontal acceleration was as large as > 0.5g. It is essential that the model is at lest partly justified, since the same source model provided synthetic accelerograms which are in a very good agreement with the measured strong motions in the other sites of Athens. The work was presented in Japan (Opršal et al., 2003), where the first author became a member of the IASPEI/IAEE working group on Effects of Surface Geology, as successor of J. Zahradník who worked in that W.G. since its commencement. After improving the geological model, the paper will be prepared for publication jointly with the co-authors from Greece. Present version is available at http://geo.mff.cuni.cz/students/oprsal/geophys/geophys.htm

All three methods developed in 2002 for numerical modeling of the strong ground motions were during 2003 applied within the EC project PRESAP, in which our team was responsible for the corresponding Work Package (W.P. 4). The project culminated by a blind experiment aimed at testing predictability of strong aftershocks in connection to Coulomb stress transfer, and simulation of the strong motions due to the aftershock. The results were collectively reported at AGU in San Francisco (Scotti et al., 2003). The kinematic method based on representation theorem, the k**(-2) stochastic slip distribution with k-dependent rise time, optionally including also asperities was re-submitted after revision, thus two papers can be expected to appear soon (Gallovič and Brokešová, submitted a,b). The composite source model of Jan Burjánek, based on fractal distribution of subevent sizes, came through further development. Significantly different directivity was found when treating subevents as point sources or finite sources. The Green function interpolation scheme over the fault plane was also developed, and increased numerical efficiency of the calculations. It was also shown how to benchmark synthetics and their peak values against the empirical attenuation curves with objective to determine a free parameter of the model, viz. the maximum slip velocity, or the effective stress drop. Significant progress in both methods was made while J. Burjánek and F. Gallovič were at the University of Kyoto, upon the invitation of Prof. K. Irikura. The work was presented both at IASPEI in Japan, and later (including results obtained during the stay in Japan, with Prof. Irikura as a co-author) at AGU in USA (Burjánek et al. 2003; Gallovič et al., 2003). A specific question, how to compare an observed strong motion record with the whole suite of the synthetics, provided by any method including a stochastic part (hence including a variety of "realizations") was addressed in a journal paper by Gallovič and Zahradník (submitted).

Two papers devoted to this event, and submitted in 2002, were accepted after re-write. The depth instability of location methods without S waves was clearly demonstrated. Trade-off between origin time and depth was eliminated by an innovating technique, the so-called station-difference method (Janský et al., in press). Problems with first-motion polarity method to retrieve the focal mechanism resulted in a recommendation to use gradient models (instead the usual models of constant velocity layers), and to employ the waveform inversion as much as possible. The resulting focal mechanism has its T axis consistent with the regional direction of extension ~N10°. However, none of the two nodal planes can be associated to the known tectonic structures. The paper, in which we considerably profited from co-operation with Greek and French colleagues, will appear soon in 2004 (Zahradník et al., in press).

The moment tensor inversion for multiple point sources, based on Kikuchi and Kanamori (1991), was extended to full waveform data at regional (or local) distances. The new code proved to be efficient for retrieving major source contributions of the 2003 Lefkada, Greece earthquake. The source model was derived from five 3-component regional stations (epicentral distances < 140 km), at periods 10-20 seconds. The model consists of two fault segments, well explaining two aftershock clusters: one at the Lefkada Island, and the other one at the Cefalonia Island, nearly 40 km apart. The Cefalonia segment started its rupture 14 seconds later. Each segment is described by two closely spaced asperities (at about 5 km from each other). The earthquake proved to be a complex rupture process, not only as regards its space-time development, but also as regards the focal mechanism. Large deviations from pure double couple (DC) were found, but interpreted as artifacts. The DC-constrained mechanism is predominantly right-lateral strike slip along of an SSW-NNE orientation, but the individual subevents had slightly different focal mechanisms. The uncertainty estimates were derived by repeatedly excluding one of the five stations. The new software name is ISOLA, to stress the interest in the Isolated Asperities. It includes not only the Fortran code (author J. Zahradník), but also the Matlab graphical user friendly interface (author E. Sokos). Dr. Sokos worked on his code with us in Prague for 6 weeks as a guest within the EC project MAGMA. The results were reported on the web pages of the European Mediterranean Seismological Center (EMSC), http://www.emsc-csem.org/ (Zahradník et al., 2003), and a journal paper is being prepared for submission.

A new section advertising earthquake software developed by our group has been established at http://seis30.karlov.mff.cuni.cz/. Although the software is only partially documented, it has yet attracted a significant attention.

The expansions proposed by Thomas (1965) for computing geodetic arc lengths and azimuths on a reference ellipsoid were analyzed, and a more accurate approach was proposed for computing azimuths (Novotný and Málek, 2003).

A review of previous deep seismic soundings in the Bohemian Massif was performed, and factors influencing accuracy of seismic interpretations were discussed (Novotný and Málek, in press). Special attention was paid to the use of quarry blasts as seismic sources, and to the application of the Wiechert-Herglotz method. It was recommended to interpret the observed data separately for the individual geological units, to apply topographic reductions, and to use suitable methods of smoothing the travel-time curves. The travel times of P-waves, measured in the West-Bohemian region during the CELEBRATION 2000 refraction experiment, were used to derive one-dimensional models for the following geological units: plutons, crystallinicum, and the Mariánské Lázně Complex. Characteristic features of the derived models are as follows: relatively low P-wave velocities at the surface, and prominent velocity increase at a depth of about one kilometre (Málek et al., submitted). It was demonstrated by Novotný et al. (submitted) that the Wiechert-Herglotz method can be used even in geologically complicated regions if the observed data are grouped according to geological units, and smoothed considerably. For the smoothing, polynomial and rational approximations were analyzed in detail. The procedure was applied to P-wave travel times from a refraction profile in western Bohemia.

Publications