Department of Geophysics, belonging to the Faculty of Mathematics and Physics, Charles University in Prague, has its roots seated as deep as in the 20's of the last century. Its structure and priorities have evolved from nearly pure seismology to the present-day broad scope covering nearly all main branches of the physics of the Earth and, most recently, comprising also planetary aspects. Research is tightly coupled with education at the bachelor, master and doctoral levels. In 2007, there were 19 staff members at the department (counting permanent and temporary positions), and 15 PhD students (13 of them supervised by the staff members).

In 2007, the Department of Geophysics participated in four EC projects: SPICE, C2C, 3HAZ and IMAGES. SPICE (2004-2007) and C2C (2007-2010) are the pan-European Marie Curie Research Training Networks involving more than 20 universities and specialized in theory of seismic wave propagation and the fate of subducting material. 3HAZ-CORINTH (2004-2007) is a research project specifically targeted on the 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.

Research project "Seismic Waves in Complex 3-D Structures" (SW3D), coordinated at the Department of Geophysics since 1993, has continued successfully in 2007. The project has been supported by 6 companies (BP America Inc., U.S.A.; Chevron U.S.A. Inc., U.S.A.; Observatorio Nacional, Brasil; Petrobras - CENPES, Brazil; Rock Solid Images, U.S.A.; Schlumberger Cambridge Research Limited, U.K.) in the framework of the international SW3D consortium.

Apart from the seismic station PRA (created in Prague in 1924), four seismic stations in Greece have been operated by the department in cooperation with the Patras University since 1997 (

Three PhD theses were defended at the department in 2007 (J. Burjanek, M. Behounkova and O. Sramek).

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

Theory of seismic waves (reported by V. Cerveny and L. Klimes)

Sensitivity of seismic waves to the structure

It has been determined how the perturbations of a generally heterogeneous isotropic or anisotropic structure manifest themselves in the wave field, and which perturbations can be detected within a limited aperture and a limited frequency band (Klimes, 2007a). Perturbations of elastic moduli and density may be decomposed into Gabor functions. The wave field scattered by the perturbations is then composed of waves scattered by individual Gabor functions. If a short-duration broad-band incident wave field with a smooth frequency spectrum is considered, the wave scattered by one Gabor function is composed of only several Gaussian packets with uniquely defined frequencies and directions of propagation. We can thus evaluate to which properties of the structure the recorded wave field is sensitive. The derived approximate solution of the forward scattering problem may be utilized in designing the migration algorithm based on true linearized inversion of the complete set of seismograms recorded for all shots.

Seismic waves in anisotropic viscoelastic media

Considerable attention has been devoted to homogeneous and inhomogeneous harmonic plane waves propagating in anisotropic viscoelastic media, and to the corresponding reflection/transmission laws (Cerveny, 2007a). In the papers by Cerveny & Psencik (2007a, in press 2007), weakly inhomogeneous plane waves propagating in weakly dissipative anisotropic media are discussed. The first-order perturbations are used to derive simple expressions for the propagation and attenuation vectors and for the polarization vectors. It is shown that the scalar product of the attenuation vector with the energy-velocity vector represents an intrinsic attenuation factor, which does not depend on the inhomogeneity of the plane wave under consideration. Actually, the intrinsic attenuation factor represents 1/Q, where Q is the generalization of the quality factor for anisotropic media. The accuracy of the perturbation method is studied by comparison of approximate and exact results.

The papers by Cerveny & Psencik (2007b, submitted 2007) are devoted to the quality factor Q for an anisotropic dissipative medium. Exact and approximate quality factors of homogeneous and inhomogeneous plane waves, propagating in generally anisotropic dissipative media, are derived from a classical definition of Buchen (1971). The approximate expression for 1/Q corresponds to the intrinsic attenuation factor, derived in Cerveny & Psencik (2007a). The same value of 1/Q is obtained no matter whether the considered wave is homogeneous or weakly inhomogeneous. Attenuation coefficients of the wave along different straight line profiles are also studied. To obtain the intrinsic attenuation factor 1/Q, the profile parallel to the energy-velocity vector must be used. The results are generalized to waves generated by point sources and propagating in heterogeneous media. The computation of Q then requires a quadrature along the reference ray. In weakly dissipative anisotropic media, Q can be expressed in terms of weak-attenuation parameters.

Time-averaged and time-dependent energy-related quantities of harmonic waves in inhomogeneous viscoelastic anisotropic media have been studied by Cerveny & Psencik (2007c). Several forms of the energy-balance equation have been derived.

In the papers by Cerveny, Klimes & Psencik (2007, submitted 2007), the attenuation vector of seismic body waves propagating in heterogeneous weakly dissipative anisotropic media is studied using the travel-time perturbation theory by Klimes (2002). The final expressions for the attenuation vector require simple quadratures along the reference ray, constructed in a relevant reference perfectly elastic anisotropic medium. The results are applicable to any type of seismic source, including a point source. It is shown that the waves in heterogeneous weakly dissipative anisotropic media are, in general, inhomogeneous. The exception are waves generated by a point source in a homogeneous isotropic medium, which are always homogeneous.

Anisotropic ray theory

Great attention has been devoted to dynamic ray tracing and ray-propagator matrices in layered inhomogeneous anisotropic media.

Cerveny & Moser (2007) and Moser & Cerveny (2007) studied the relations between the 6x6 ray-propagator matrices in Cartesian and ray-centred coordinates, and found useful transformations between them. They also derived the 6x6 and 4x4 interface propagator matrices in ray-centred coordinates.

Cerveny (2007b) compared the dynamic ray tracing systems, consisting of four linear ordinary differential equations of the first order, in ray-centred coordinates and in local wavefront orthonormal coordinates. The systems, known form literature and derived independently, consist of equations expressed in different forms. Cerveny showed that both systems are fully equivalent and may be used alternatively.

Cerveny & Moser (submitted 2007) derived initial conditions for dynamic ray tracing system for heterogeneous anisotropic media. The initial conditions correspond to an arbitrarily curved initial surface, to a point source and to an arbitrary curved initial line.

Coupling ray theory

Coupling ray theory is required for modelling propagation of S waves in heterogeneous weakly anisotropic media (Bulant & Klimes, 2007d).

Klimes (2007b) proposed the coupling ray series, which yields the coupling ray theory in a similar way as the standard ray series yields the standard ray-theory solution of the elastodynamic equation.

Bulant & Klimes (2007b, 2007e) numerically tested the equations for the anisotropic-common-ray approximation of the coupling ray theory in three smooth 1-D velocity models of differing degree of anisotropy. They calculated the synthetic seismograms obtained by both the isotropic-common-ray and anisotropic-common-ray approximations of the coupling ray theory, and compared these seismograms with the more accurate coupling-ray-theory synthetic seismograms simulated by the second-order perturbation expansion of travel time from the anisotropic common rays.

Velocity macro models and numerical ray tracing

Bulant & Klimes (2007a, 2007c, submitted 2007) showed how to calculate the sonic-log travel times in the geological structure from the sonic-log velocities while taking into account the effects of the non-vertical propagation of seismic waves due to the source offset and due to the heterogeneous velocity in the structure, together with the effects of deviation of the well trajectory from strictly vertical. These travel times serve for comparison of sonic logs with vertical seismic profiling.

Bulant, Eisner, Psencik & Le Calvez (2007) paid attention to the well deviation surveys. Not performing accurate well deviation surveys for hydraulic fracture monitoring and neglecting the effects of the well trajectory results in significant errors in the calculated fracture azimuth and other parameters. For common geometries, a 5-degree deviation uncertainty of monitoring or treatment wells can cause more than 40-degree uncertainty in inverted fracture azimuths. The errors caused by the deviation uncertainty should not be compensated by artificially adjusting the velocity model.

CD-ROM with SW3D software, data and papers

Compact disk SW3D-CD-11 (Bucha & Bulant, 2007) contains the revised and updated versions of the software developed within the consortium research project "Seismic Waves in Complex 3-D Structures" (SW3D), together with input data used in various calculations. Compact disk SW3D-CD-11 also contains over 270 complete papers from journals and from SW3D consortium research reports.

Earthquake and structural studies (reported by J. Zahradnik)

Studies of earthquakes in Greece

The international project 3HAZ-CORINTH of the 6th framework program of EU has been completed in 2007 (Bernard et al., 2006). It played a fundamental role in strengthening our collaboration with Greek and French scientists focused on earthquakes, landslides and tsunamis in the Gulf of Corinth. Three seismic stations were operated jointly by the Charles University and the University of Patras in the framework of the project (Sergoula, Mamousia and Loutraki). The latter, LTK station, equipped with satellite telemetry, became in 2007 an on-line data contributor to the international data center ORFEUS. Some remaining technical problems related to the data flow through NAQS server will be solved in 2008.

The program package ISOLA designed for inversion of regional waveforms into multiple point source model was accepted for publication (Sokos and Zahradnik). In cooperation with the Seismological Laboratory, University of Patras, the package has been routinely used to determine moment tensors of earthquakes in western Greece, and to post the results automatically via e-mail onto web pages of the European Seismological Center (EMSC). The ISOLA software has been continuously upgraded during the year, with the most recent version freely available at the Patras-Prague joint web page ( Particular attention was paid to the non-double-couple components arising from the earthquake complexity, studied both in synthetic tests and real data of three M5 Zakynthos 2006 earthquakes (Zahradnik et al., 2007).

Modeling the dynamic stress field

An important first step to enrich the scope of the department by studies of the source dynamics has been made by J. Burjanek, who defended his PhD thesis in 2007 (Burjanek, 2007). His analysis of dynamic stress field associated with the slip history prescribed in the advanced kinematic models was published in GJI (Burjanek and Zahradnik, 2007). It is a detail numerical study of the full dynamic stress history on a fault by the boundary element method, formulated in spectral domain. Basic characteristics of the resulting stress field were discussed in relation to recent dynamic source models of real earthquakes. The study included the static stress drop, the strength excess and the dynamic stress decrease distributions corresponding to the kinematic k-square stochastic source model with k-dependent rise time. A new quantity, the stress delay, was introduced.

Modeling strong ground motions

Gallovic and Burjanek (2007) compared two common approaches to strong ground motion simulations due to finite extent seismic sources: the integral and composite approach. It has been shown that the integral model gives larger scatter of peak ground acceleration than empirical models, which is explained by the overestimated directivity effect. In Gallovic and Brokesova (2007a) a new hybrid approach was proposed and applied to modeling of the 1997 M6.1 Kagoshima earthquake, Japan, and the 1999 M5.9 Athens earthquake, Greece. The hybrid modeling has been also applied to the probabilistic aftershock hazard assessment (PAHA), see Gallovic and Brokesova (2007b, c), accepted for publication in the J. Seismology. Further, an influence of non-planar faults on strong ground motions, so far, one of the first studies of this kind, was studied in the paper by Kaeser and Gallovic (2007), accepted for publication in GJI.

Modeling 3D local site effects (finite-difference method)

Local 3D site effects due to surface topography and complex geometry of the bedrock interface were investigated by hybrid method developed in recent years by I. Oprsal. His software, able to manage realistic source-path-site configurations and wavefields of the engineering interest (up to 10 Hz), has attracted a number of scientists and extended the collaboration of the department quite significantly. The numerical investigation of the Roman city of Augusta Raurica, Switzerland that probably suffered from a strong earthquake in the year 250 AD, was published (Oprsal and Faeh, 2007). Parallel version of the code, developed in cooperation of I. Oprsal and A. Caserta (the external PhD student of the Charles Univ.) has been developed, further tested and submitted in Computers and Geosciences.

I. Oprsal, on leave from Charles University since 2007/01 at the Disaster Prevention Research Institute DPRI, (Section of Strong Motion Seismology, Division of Earthquake Disaster Prevention), Kyoto University, Japan, participated in the investigation of the long-period ground motions in the Osaka basin that are important for understanding and predicting effects of large subduction earthquakes of the Nankai trough. His 3D modeling of an M7 intermediate-depth earthquake of Oct. 12, 1993 proved the existence of the basin generated surface waves; complex pattern of the source and basin effects will need more investigation in 2008.

Structural studies

We participated in structural studies in several regions of the world. Spicak et al. 2007) investigated Wadati-Benioff zone of Central America. Jansky et al. (2007a,b) constructed new 1D models of the upper crust in the Corinth Gulf, Greece. The abundant data of the CRL network allow also nonstandard assessment of the velocity structure based on nonlinear inversion of the P and S wave travel times by the Neighborhood Algorithm. Kolinsky and Brokesova (2007) used surface waves to study structure in Western Bohemia, where they applied sophisticated new software for the frequency-time analysis. Kampfova et al. (2007) detected strong Moho reflections in the Ore Mountains, Czech Republic, and determined the ratio of the P- to S-wave velocities in the crust of the region. Hanzlova et al. (2007) located weak seismic events at the Orlik water reservoir region. Holub and Novotny (2007) suggested to approximate the travel-time curves of refracted waves by rational functions (quotients of two polynomials), and then to use the classical Wiechert Herglotz method to interpret the smoothed travel-time curves.

Geodynamics (reported by C. Matyska)

Slab dynamics modelling

Cizkova et al. (2007) studied the deformation of the old slabs in the transition zone and especially the stress distribution within the slabs. They concluded that the stress field is mainly affected by bending forces and to a lesser extent by the effects of upper mantle phase transitions. The stress distribution is characterized by a layered structure similar to the structure of double-seismic zones reported in shallow portions of some slabs. Behounkova and Cizkova (2008) attempted to determine the basic mechanical properties of subducting slabs in the lower mantle and to answer the question under what conditions the thick, blob-like structures as those observed by tomography could be produced in the numerical models of subduction. The results indicate that the slab thickening is observed only if the relatively weak slabs (stress limit 0.5 GPa) are considered and a viscosity increases by a factor of 10 or 30 in the lower mantle. Depending on the strength of the crust, two types of thickened structures are observed in the lower mantle. In the case of a relatively weak crust ensuring a low friction on the contact of the subducting and overriding plates, buckled instabilities form in the lower mantle. In case of a strong crust, slabs tend to thicken by compression and conductive cooling.

Geoid over the subducting plates

Tosi and Martinec (2008) tested various viscosity models of the mantle and different shapes of slabs in the lower mantle and compared the predicted surface observables (dynamic geoid and topography) with the topography and gravity signals around subduction zones. They concluded that the lower mantle viscosity controls the low-degree geoid to the first order, ensuring a positive geoid signal only when it is at least 50-times greater than the upper mantle viscosity. Very localized LVV associated with the cold slab play a minor role on geoid amplitudes, suggesting that slabs do not need to remain stiff and do not efficiently transmit stresses from deep mantle to the Earth surface. Finally, the presence of anomalous slab density in the lower mantle is always needed to obtain geoid amplitudes comparable with the observations, favoring the picture of slabs that penetrate the transition zone and sink into the deep mantle.

Thermal convection in the mantle with post-perovskite phase change and radiative heat transfer, lower mantle viscosity structure

In cooperation with David A. Yuen from Minnesota Supercomputing Institute we presented the results of numerical mantle convection models demonstrating that dynamical effects induced by variable mantle viscosity, depth-dependent thermal expansivity, radiative thermal conductivity at the base of the mantle, the spinel to perovskite phase change and the perovskite to post-perovskite phase transition in the deep mantle can result in multiscale mantle plumes: stable lower-mantle superplumes are followed by groups of small upper-mantle plumes. Both radiative thermal conductivity at the base of the lower mantle and a strongly decreasing thermal expansivity of perovskite in the lower mantle can help induce partially layered convection with intense shear heating under the transition zone, which creates a low-viscosity zone and allows for the production of secondary mantle plumes emanating from this zone. Large-scale upwellings in the lower mantle, which are induced mainly by both the style of lower-mantle viscosity stratification and decrease of thermal expansivity, control position of central uppermantle plumes of each group as well as the upper-mantle plume-plume interactions.

From the depth-dependent thermal expansivity, we deduced the 3-D density anomalies from the seismic velocity anomalies inferred from seismic tomographic inversion. Using these density distributions, we calculated the viscous responses of the Earth due to these density anomalies for a given viscosity structure. We then focused on the lateral viscosity variations of the deep mantle on the solution of the inverse problem involving the inferences of the viscosity from the long-wavelength geoid. Our solution for the large-scale lateral viscosity structure in the lowermost mantle shows that the region underneath hot spots have significantly higher viscosity in the deep mantle than the region below subduction regions.

Slab interaction with the post-perovskite phase transition

Using a simplified numerical model of mantle convection and assuming a non-Newtonian nature of the post-perovskite, Cizkova et al. (2007) demonstrated that a viscosity paradox mentioned in the previous section is dynamically feasible. Assuming reasonable rheological parameters of the lower mantle perovskite and post-perovskite they get very weak post-perovskite lens in the cold slab material, while the viscosity in the neighbouring perovskite plumes is higher despite their considerably higher temperature. Further, the dynamics of the core-mantle boundary region is considerably influenced by the assumed rheological properties of the post-perovskite.

Anisotropy of the asthenosphere

Marquart, Schmeling and Cadek (2007) used different types of numerical models of mantle convection to predict the distribution of seismic anisotropy in the asthenosphere in the North Atlantic region. The predictions were compared with observed seismic anisotropy patterns and discussed in connection with the evolution of the Mid-Atlantic Ridge and the Iceland Plume.

Electrical conductivity of the Earth’ mantle

Velimsky et al. (submitted to GJI) have developed a new approach to the 3-D inverse problem in time-domain. The method employs the conjugate gradient algorithm with efficient evaluation of the gradient of misfit in the model space by adjoint approach. Mathematical formulation has been derived for three different boundary conditions and the code has been validated using synthetic checkerboard models.

Ice sheet modelling

Soucek & Martinec (2008) developed a novel technique, called the SIA-I algorithm, for computing the Stokes solution of glacier flow, which is based on the traditional scaling assumptions of the Shallow Ice Approximation. The algorithm significantly reduces the CPU-time for computing the Stokes glacier flow if it is compared to other full-Stokes or higher-order numerical solvers. Moreover, CPU-time grows only linearly with the number of computed field variables. In near future, we will aim at incorporating this fast solver into large-scale glaciological models.

The SIA-I algorithm was successfully tested for both sliding and non-sliding boundary conditions at the ice bedrock in the frame of the international benchmark project ISMIP-HOM for glaciological models of higher-order accuracy (Pattyn et al., 2008). With the cooperation with O. Rybak from Alfred Wegener Institute fur Polar und MeereForschung, the SIA-I was used to compute the Stokes glacier flow in the region of Dronning Maud Land in Antarctica. For a highly irregular ice-bedrock topography with abrupt changes in sliding velocities and with the rheology strongly dependent on temperature, the SIA-I techbique provides sufficiently accurate results if compared to second-order models, but with substantially lower (1-2 orders of magnitude) computational costs.

Visualization of postglacial rebound

Hanyk et al. (2007) applied the 3D visualization package Amira to depict both the external and internal deformation histories of the transient viscoelastic flow inside the mantle induced by postglacial uplift. Of particularly great interest are the transient displacement fields and shear heating inside the mantle. They also integrated the visualization results into the Google Earth virtual globe by combining this scheme with the Amira package to provide a better geographical and dynamical context.

Earth gravitational field from satellite data

Sasgen, Martinec and Fleming (2007) inferred regional mass changes in Antarctica using 4 years of Gravity Recovery and Climate Experiment (GRACE) data. They decomposed the time series of the Stokes coefficients into their linear as well as annual and semi-annual components by a least-squares adjustment and applied a statistical reliability test to the Stokes potential-coefficients' linear temporal trends. Mass changes in three regions of Antarctica that display prominent geoid-height change were determined by adjusting predictions of glacier melting at the tip of the Antarctic Peninsula and in the Amundsen Sea Sector, and of the glacial-isostatic adjustment (GIA) over the Ronne Ice Shelf. The mass-change estimates derived from all GRACE releases and time intervals lie within 20% (Amundsen Sea Sector), 30% (Antarctic Peninsula) and 50% (Ronne Ice Shelf region) of the bootstrap-estimated mean, demonstrating the reliability of results obtained using GRACE observations.

Planetary physics

Planetological research at the department is represented by two papers prepared in a close cooperation with researchers of the Planetology Institute in Nantes, France. Choblet, Cadek, Couturier and Dumoulin (2007) tested a new finite volume code to simulate thermal convection in planetary mantles and its topographic and gravitationa response. Tobie, Cadek, and Sotin, (Icarus, in press) deals with the anomalous heat flux observed on the South Pole of Saturn's moon Enceladus. Assuming a viscoelastic model, the authors demonstrated that the high heat flux can be explained by a tidal deformation of Enceladus' icy mantle provided that a liquid water is present at the base of the mantle.