Professor Vlastislav Cerveny was honored by Society of Exploration Geophysicits with the Maurice Ewing Medal in recognition of his lifetime internationally recognized work in the ray theory of wave propagation.

The MAGMA center, i.e. the Prague Center of Mathematical Geophysics, Meteorology, and their Applications, supported by the European Commission, has completed in 2004 its second successful year of existence. The geophysical group profited from short visits as well as from more than 1 month lasting stays of 13 senior researchers, 7 post-doctoral researchers and 13 PhD students (23 person-months altogether). For more details, see http://geo.mff.cuni.cz/magma. In a close relation with MAGMA, we entered into three new 6^th-framework European projects: SPICE, 3HAZ-CORINTH and IMAGES. SPICE (2004-2007) is the pan-European Marie Curie Reseach Training Network involving 14 universities and specialized in theory of seismic wave propagation. 3HAZ-CORINTH (2004-2006) is a targeted research project oriented to three main natural hazards in the Gulf of Corinth, Greece, viz earthquakes, landslides and tsunamis. IMAGES (2005-2008) aims at transfer of knowledge between seismologists and applied geophysicists (Schlumberger Cambridge Research), studying microearthquakes induced by oil drilling. Intensive cooperation with several major oil companies, lasting since 1993, 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.

For a spherically symmetric viscoelastic Earth model, the movement of the rotation vector due to surface and internal mass redistribution during the Pleistocene glaciation cycle has conventionally been computed in the Laplace-transform domain (Martinec and Hagedoorn, submitted). The new method offers the possibility to model the rotational response of the Earth induced by surface glacial loading by numerical time integration of the linearized Liouville equation. The theory extends the conventional approach based on the second-degree load Love numbers to general 3-D viscoelastic earth models. The time-domain solution of both the glacial isostatic adjustment and the induced rotational response of the Earth is easy to combine with a time-domain solution of sea-level equation with a time-varying shoreline geometry. An axisymmetric viscoelastic model was employed to explain Fennoscandian relaxation data (Martinec and Wolf, in press).

Moreover, we have studied the possibility of short time-scale energy transfer from the ice sheet loading and unloading processes to the Earth's interior via viscous dissipation associated with the transient viscoelastic flow in the mantle (Hanyk et al., submitted). We have focussed on the magnitude of glacially-induced deformations and the corresponding shear heating for an ice sheet of the spatial extent of Laurentide region in Maxwellian viscoelastic compressible models with a Newtonian viscosity. We have found that shear heating from the transient viscoelastic flow can represent a non-negligible mantle energy source with cryogenic origins.

We also investigated global geoid change resulting from fluctuations in the spatial and temporal configuration of several regions of major present-day glaciation (Fleming et al., 2004). The areas of primary interest are Greenland, Antarctica, Patagonia, Alaska and Iceland. The contribution to geoid displacement due to the ongoing isostatic adjustment of the Earth to past ice-load changes was examined, in particular those following the Last Glacial Maximum, but also more recent events, e.g. the late Holocene and Little Ice Age. For present-day predictions, we employ recently published descriptions of changes in these glaciated areas that are a product of a variety of data types. Based on an optimistic accuracy estimate of the new satellite mission GRACE, the geoid displacement in these areas is detectable up to degree and order 48, however, this is still dependent upon the final accuracy achieved.

Earth's mantle convection, lithospheric subduction Studies on generation of slab-like downwellings in the global layered mantle convection models and influence of the zone of weakness on the dip angle of subducting lithosphere, submitted in 2003, were finally published (Cizkova and Matyska, 2004; Kukacka and Matyska, 2004). Combined effects of kinematic and free-slip boundary conditions, age of the plate and trench migration rate on the stresses within the slabs in 2-D Cartesian numerical models of lithospheric subduction have been modeled in a broad parameter space and publication of results is under preparation. Convection models with a hypothetical interface at a 1000 km depth have been also used as synthetic models in 2-D and 3-D synthetic tomographic inversions (Behounkova et al., submitted). The resolving power of tomography in different parts of the mantle has been described in detail. We have found that the role of projection (discretization) error, usually omitted in standard resolution studies, may be substantial.

Dynamic consequences of the recently discovered post-perovskite phase in the deep mantle were studied in (Matyska and Yuen, in press; Matyska and Yuen, submitted). They have investigated the impact arising from the interaction of temperature-dependent and depth-dependent viscosity with radiative thermal conductivity on the dynamics of both the ascending and descending flows in the presence of both the phase change at 670 km depth and post-perovskite transition at 2650 km depth. The importance of radiative thermal conductivity in maintaining the coherent structures of large upwellings in the lower mantle was demonstrated. These results also revealed that a greater degree of asymmetry is produced in the vertical flow structures of the mantle by the phase transitions. Therefore, mass and heat-transfer between the upper- and lower mantle will deviate substantially from the traditional whole-mantle convection model.

In order to overcome a lack of data that could constrain mantle flow models, a new insight into dynamics of the asthenosphere is now expected from comparing the plate-driven flow models with the data on seismic and electrical anisotropy in the sublithospheric mantle. The research, combining different types of mantle flow modeling with both electromagnetic and seismic data, has been recently started by O. Cadek. It is focussed on the Australian plate, which is relatively fast and for which a unique set of seismic and electromagnetic observations is available.

Effects of lateral viscosity variations in the lower mantle on the long-wavelength geoid and surface topography have been studied (Cadek and Fleitout, submitted). Information on gravitational and topographic signal was also used to investigate the internal structure of Venus (Pauer, 2004). The main results is that the viscosity in the Venusian mantle must increase with depth by several orders of magnitude. This conclusion differs from previous studies by other authors who mostly argued for only a week increase of viscosity with depth.

Forward modeling of electromagnetic induction was implemented in the time domain for satellite data (Martinec and McCreadie, 2004; McCreadie and Martinec, 2004). The input of the model consists of the X-component of the magnetic induction vector recorded by the CHAMP and/or Oersted vector magnetometer along individual night-time, mid-latitude satellite tracks during intense geomagnetic storms. The output is the Z-component of the magnetic induction vector at the satellite altitudes, which is used as misfit function in the inverse formulation. Conductivity in the mantle under the Pacific region is sought by the local 2-D inverse modeling of electromagnetic induction. The results of the inverse modeling suggest an increased mantle conductivity in the southern Pacific. Similar approach can be used below other areas to construct a preliminary 3-D conductivity model of the upper mantle. Fully 3-D responses of conductivity model are also verified (Velimsky and Martinec, in press).

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

Recent developments in seismic ray method

Important developments in seismic ray method achieved during the last 20 years have been described in an extensive review paper by Cerveny, Klimes & Psencik (submitted). The topics include ray histories, two-point ray tracing, controlled initial-value ray tracing, wavefront tracing, interpolation within ray cells, paraxial ray methods, third-order and higher-order spatial derivatives of travel time, second-order and higher-order perturbation derivatives of travel time, optimization of model updates during linearized inversion of travel times, coupling ray theory for S waves, quasi-isotropic approximations of the coupling ray theory, Gaussian beams, Gaussian packets, optimization of the shape of Gaussian beams or packets, asymptotic summation of Gaussian beams and packets, linear canonical transforms, coherent state transforms, Maslov methods, decomposition of a general wave field into Gaussian packets or beams, sensitivity of waves to heterogeneities, Gaussian packet migrations, higher-order ray-theory approximations, direct computation of first-arrival travel times, ray method with complex eikonal, hybrid methods, ray chaos, Lyapunov exponents and rotation numbers, models suitable for ray tracing, application of Sobolev scalar products to smoothing models.

Considerable attention has been devoted to harmonic plane waves in viscoelastic anisotropic media (Cerveny 2004a; Cerveny & Psencik, 2004a, 2004b, submitted-a, submitted-b), and to the corresponding reflection/transmission coefficients (Cerveny, 2004b).

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 wave field due to caustics have been derived (Klimes, submitted).

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 (Klimes & Bulant, 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 & Klimes, 2004b; Klimes & 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 & Klimes, 2004b; Bulant, Klimes, Psencik & Vavrycuk, 2004; Klimes, 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, coded, numerically tested, and applied to the calculation of coupling-ray-theory seismograms (Bulant & Klimes, 2004a).

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. 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 (Bulant, 2004). Capabilities of the ray tracing software have further been extended (Bucha & Bulant, 2004; Bulant & Klimes, 2004a). Numerical ray tracing has been tested on various 3-D synthetic structures, including the smoothed SEG/EAGE Salt Model (Bucha, 2004a, 2004b, 2004c; Bucha & Bulant, 2004).

Equations for the propagation of Gaussian packets in smooth isotropic media have been derived in the form suitable for numerical applications (Klimes, 2004b). Klimes (2004c) and Klimes & Zacek (2004) demonstrated that the summation of Gaussian beams and packets is considerably comprehensive and flexible, and may be formulated in many ways. The form of the summation depends primarily on the specification of the wave field. Asymptotic summation of Gaussian beams or packets also includes the Maslov method and its various generalizations, based on canonical transforms or coherent state transforms, as special cases. Decomposition of a general wave field into Gaussian packets also includes the system of Gaussian packets scattered from Gabor functions forming medium perturbations.

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 (Zacek, 2004a). Gaussian-packet prestack depth migration of the decomposed time sections has been proposed and first numerical tests have been started (Zacek, 2004b, 2004c).

The present network, developed in cooperation with the Patras University, comprises four sites, each one equipped with a weak-motion broad-band velocigraph CMG 3-T and a strong-motion accelerograph CMG 5-T. The selected data are available from http://seis30.karlov.mff.cuni.cz, updated every 4 months. The stations Sergoula and Mamousia are situated on the northern and southern coast of the Corinth Gulf, respectively, both in its western part. They are operated as stand-alone stations, and (since 2004) belong to the new EC project 3HAZ-Corinth, coordinated by P. Bernard. The other two sites have satellite data transmission to Patras. It is Loutraki station, at the eastern edge of the Corinth Gulf, and Pylos station, close to Kalamata city, on the south-west of the Pelloponesos.

Main use of the data in 2004 was to clarify how the weak and strong motion instruments efficiently complement each other (Zahradnik, 2004). Nature of strange long period pulses on CMG 3-T broadband records due to nearby earthquakes. was also investigated. The pulses have been identified as normal response to a bit abnormal ground motion input, viz a sudden (step-like) horizontal acceleration, most likely connected with a local tilt provoked by the vibratory seismic motion in the immediate vicinity of the seismic instrument (Zahradnik and Plesinger, submitted). The paper also shows that similar effects (although less obvious) may significantly obscure the Le-3D/20s broadband records of the Greek national network. A simple way how to detect the long-period disturbances, and how to "clean" the records prior their use in seismic source studies has been suggested.

The 3-D modeling based on a hybrid combination of the source, path and site effects, methodically developed since 2002, has been applied in practice. The source and path effects are modeled by the composite-source model and the discrete wavenumber method (method PEXT of J. Zahradnik), while the local sedimentary basin are modeled by the 3-D finite-difference method (I. Oprsal). 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 Liosia site during the 1999 Athens earthquake as a combined effect of the source directivity and lateral heterogeneity of the site. Results were presented in Japan, Germany and Canada, and published (Oprsal et al., 2004).

The same method was successfully applied to predict earthquake effects of a hypothetical future strong earthquake in Basel, Switzerland (Oprsal et al., submitted).

Several versions of the strong-motion simulation technique have been tested and further developed.

The composite method, based on fractal distribution of subevent sizes, with Green function interpolation, has been upgraded by J. Burjanek.

The kinematic method based on representation theorem, the k**(-2) stochastic slip distribution with k-dependent rise time, optionally including also asperities, was published (Gallovic and Brokesova, 2004a). A parametric study and application to the 1999 Athens earthquake was also completed (Gallovic and Brokesova, 2004b). J. Gallovic concentrated not only on his attempts how to make the source directivity modeling realistic enough, e.g. by combining the composite and kinematic models, but also how to employ the finite-extent source numerical ground motion simulation into probabilistic hazard assessment schemes.

A quite new research line was followed by J. Burjanek, who developed his codes for calculating dynamic stress field on the fault whose slip is described by a kinematic source model. It is important for understanding rheology of the faulting process, e.g. for testing applicability of the broadly used slip weakening rheological model.

Grid search and a recently developed innovation of the seismic location method, the so-called station-difference method, were applied to the Egion M4.3 earthquake of 2001, Corinth Gulf (Jansky et al., 2004). Focal mechanism of the same earthquake was retrieved by inverting amplitude spectra at three local stations (Zahradnik et al., 2004).

A 1-D gradient model for the Egion region was retrieved by J. Jansky and V. Plicka from the Corinth Rift Laboratory (French-Greek) data using the Neighbourhood Algorithm.

A new inversion code to simultaneously locate earthquakes and to retrieve a 1-D crustal model was developed by O. Novotny, based on gradient methods. First attempts were made to apply the method in the Corinth Gulf.

Master-event relocation of weak earthquakes in Western Bohemia was published by (Jansky and Malek, 2004).

The moment tensor inversion for multiple point sources, based on full waveform data at regional distances (code ISOLA) was further developed by J. Zahradnik. The code was successfully applied to the M6.3 earthquake at Lefkada, Greece, 2003. 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 and 14 seconds later. The earthquake proved to be a complex rupture process, not only as regards its space-time development, but also as regards the focal mechanism. The paper came through revision and was accepted for publication (Zahradnik et al., in press). Ch. Benetatos from the University in Thessaloniki worked with us in Prague for one month as a guest within the EC project MAGMA, and ISOLA code was also applied to the M6.5 earthquake at Skyros Island, Greece, 2001.

We participated in a new phase of the structural studies of the Earth's crust in our country by the method of Deep Seismic Sounding method (DSS). We concentrated on the interpretation of the DSS data from seismo-active region of West Bohemia. The available experience and possibilities have been reviewed by Malek et al. (2004b). New methods of smoothing travel time curves, based on polynomial and rational approximations were developed (Novotny et al., 2004a). Data interpretation for three individual geological blocks in Western Bohemia was separately performed, and 1-D velocity models of the upper crust were derived (Malek et al., 2004a). Fast velocity increase in the uppermost 1 km was found, not detected in the previous DSS studies. Analogic structural studies based on body waves have been developed for the Moravia-Silesia region, too.

A theoretical approach based on dispersive surface waves was applied to measurements from the eastern and northeastern Moravia (Novotny et al., 2004b).

Publications