Research overview 2011


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 and geomagnetism 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 2011, there were 19 staff members at the department (counting permanent, temporary and part-time positions), and 18 PhD students (11 of them supervised by the staff members).


Research project SW3D, Seismic Waves in Complex 3-D Structures, coordinated at the Department of Geophysics since 1993, has continued successfully in 2011.  The project has been supported by 5 companies or research institutes (BP America Inc., U.S.A.; Chevron U.S.A. Inc., U.S.A.; NORSAR, Norway; Petrobras, Brazil; Schlumberger Cambridge Research Limited, U.K.) in the framework of the international SW3D consortium. The project commisioned by ESA which aims to develop and test the Swarm Level 2 processing facility started in October 2010. The Charles University, as a member of the SMART consortium, is providing support in development of the time-domain chain for inversion of Swarm data in terms of 3-D mantle conductivity.


Seismic station PRA (created in Prague in 1924), equipped with CMG-3T, is linked with the international data center ORFEUS to provide on-line data transfer to Utrecht. Joint project of the Charles University in Prague and the University in Patras, initiated in 1997, has successfully continued. Current configuration of the jointly operated seismic stations is as follows: SERG (Sergoula), LTK (Loutraki), PYL (Pylos), PDO (Prodromos), ZKS (Zakynthos), ANX (Ano Chora). For details, see Each site is equipped with a pair of the broad-band and strong-motion instruments. All instruments are of Guralp manufacturer, except ZKS (Nanometrics). The stations became a part of the Hellenic Unified Seismic Network (HUSN), with data sharing by three universities (Athens, Thessaloniki and Patras) and the National Observatory of Athens. The participating institutions in Greece make everyday use of the stations in their earthquake locations.  Moreover, the broadband waveforms are used to routinely invert moment tensors at the University of Patras and report them to EMSC. Software ISOLA, continually upgraded (, is used for this purpose. Consultations to many external users of the software are provided, too. The department also participated in operation of the borehole (60 m deep) broadband CMG-3TD station GOPC located at the Geodetic observatory Pecny (GOPE).  


Operation of all these stations is partly supported by the project CzechGeo/EPOS, see . Within this project the department operates also the CzechGeo/EPOS Seismological Software Centre, see


The Editors or Associated Editors of international journals have been Vlastislav Červený (Russian Geology and Geophysics), Luděk Klimeš (Journal of Seismic Exploration), Ctirad Matyska (Journal of Geophysical Research – Solid Earth; Studia Geophysica et Geodaetica) and Jiří Zahradník (Journal of Seismology).  Ctirad Matyska has been the member of the editorial board of the journal Pokroky matematiky, fyziky a astronomie, which is devoted to popularization of mathematics, physics, astronomy and education in these fields. Jakub Velímský has been nominated an "Outstanding Reviewer 2011" by the Geophysical Journal International.


The department was happy to host visiting professor Craig Bina from Northwestern University Evanston during his three months research leave.


The  theses defended at the Department in 2010 included the following:

B.C.: Jakub Michálek, Ondřej Peisar

M.S.: Miroslav Kuchta, Helena Munzarová, Pavla Procházková, Libor Šachl

Ph.D.: Arrigo Caserta


As in previous years, research at the Department was carried out in three directions: Theory of Seismic Waves,  Earthquake and Structural Studies, and Geodynamics.


Earthquake sources (reported by Jiří Zahradník)


The following events were studied (Mw denotes the moment magnitude): Mw 6.3 Movri Mountain 2009 Greece, two earthquakes Mw>5 Efpalio 2010 Greece, Mw 6.3 L’Aquila 2009 Italy, Mw 6 Silakhor 2006 Iran, Mw 6.0 Horseshoe 2007 west of Gibraltar, Mw 9.0 Tohoku 2011 Japan, Mw 7.2 Van 2011 Turkey. Each earthquake was studied in one or more papers; see below. In case of Japan (Zahradník et al., 2011) and Turkey (Zahradník and Sokos, 2011) we presented research reports to the European Mediterranean Seismological Center (EMSC), posted on Internet in the framework of quick post-seismic investigations. We focused on the moment tensors, position and time of the centroid, as well as on the identification of the fault plane. We also tried to construct a kinematic model of the spatial-temporal evolution of the fault process. In case of the Efpalio 2010, Corinth Gulf, earthquake we also participated in the relocation of the whole sequence (Janský et al., 2012, in press) and in the seismotectonic interpretation (Sokos et al., 2012, in press); for the first time we also cooperated with geodesists in using the GPS measurements of co-seismic surface motion in Efpalio, Greece (ibid.).


We also focused on theoretical and practical problems of the inverse problem, when seismograms are used (inverted) to compute the kinematic model of the fault process. We proved (Gallovič and Zahradník, 2011) that in certain cases the result may significantly differ from the true process. It is because the inverse problem is ill-posed, and the necessary regularization can alter not only details of the spatial-temporal slip pattern, but also its gross features; the co-called artifacts arise. Such a solution (although wrong) may be very stable, which many authors erroneously interpret as indication of the correctness. We learned how to avoid such misinterpretations by a careful combination of real data and synthetic tests. The knowledge was applied to study the damaging earthquake in L’Aquila, Italy 2009 (Gallovič and Zahradník, 2012, in press; Ameri et al., 2012, in press).


The cooperation with INGV Milan further continued. Ameri et al. (2011) compared three ground-motion simulation techniques, including the broadband hybrid integral-composite technique of Gallovič and Brokešová (2007), and investigated their performances in near-fault strong-motion modeling. The test case focused on the 1980 M6.9 Irpinia earthquake, the strongest event recorded in Italy in the last 30 years. Important engineering parameters were analyzed, such as PGA, PGV and spectral accelerations.


A new method was developed to study the moment tensor resolvability under the assumption of a fixed position and time of the source. Matrix of this linear inverse problem can be analyzed in terms of its eigenvectors and eigenvalues, which makes it possible to determine the likely bounds of the source angular parameters (strike, dip and rake) for a given model of the Earth’s crust and a given source-station configuration. Data (seismograms) are not needed for this purpose, thus the method is suitable for a network design. This issue was analyzed in Zahradník and Custodio (2012, in press), where we presented theoretical resolvability maps for various configurations of seismic stations in Iberia (Portugal, Spain, Morocco), including the planned OBS instruments. The study included analysis of a non-standard case of a single-station inversion. It was shown why the single-station inversion can work (e.g. for the Horseshoe Mw 6.0 earthquake, 2007), and why in some other cases it may fail. The success or failure is determined by a (complex) joint effect of the source depth and the frequency range. This result shows that the standard (almost ‘mythical’) requirement of a good azimuthal coverage is not always a must, at least not for the double-couple part of the solution.


The ISOLA software for multiple point-source moment tensor inversion of local and regional seismograms,, continually upgraded in cooperation with E. Sokos (University of Patras, Greece), found many new users in 2011, e.g. in connection with a one-week course organized by the Universidad de Heredia, Costa Rica. The course had two teachers (J. Zahradník and E. Sokos) and was attended by 15 students from 8 countries of the Latin America. Numerous e-mail consultations were provided to other ISOLA users during 2011, too.


Independently, also the moment-tensor inversion in the spectral domain was developed by M. Pakzad, a former PhD student in Prague, who finally received his degree in Tehran. In 2011 we re-established the cooperation through investigation of the Mw 6 Silakhor 2006 Iran earthquake and its two aftershocks.


New instruments

A newly invented mechanical sensor system, called Rotaphone, for recording rotation rate components through measurements of spatial ground motion gradients was further developed. Two new prototypes were subjected to elaborate testing at the Albuquerque Seismological Laboratory (ASL), U. S. Geological Survey, having excellent facilities to test rotational seismic sensors over a wide frequency band. The results of these tests together with first application of the new instruments in monitoring local shallow earthquakes in Western Bohemia seismoactive region were reported by Brokešová et al. (2012). A new six-degree-of-freedom prototype for recording three translational and three rotational components was deployed in the vicinity of Provadia salt-works (Bulgaria) for several months. Collocated rotational and translational recordings written both by local induced seismic events and stronger earthquakes from Romania, Bulgaria and Western Turkey are studied with a focus on various seismic rotation effects. This study is still in progress and will continue during 2012.


Structural studies

Shear wave velocity-depth sections were investigated in four geological units in western part of the Bohemian Massif (Kolínský et al., 2011). The seismic models were retrieved by detailed analyses of phase velocities of fundamental Love wave modes, using three earthquakes from the Aegean Sea at periods above 12 seconds. Undulations of the phase velocity dispersion curves were discussed as caused by frequency-varying backazimuths and the surface-wave mode coupling. Corrections for the true backazimuths were performed to partly reduce the undulations. Four sub-profiles crossing three main geological units were studied and compared with previous investigations (Moldanubian, Teplá-Barrandian, and Saxothuringian units). Special attention was devoted to the seismically active area of the Eger Rift (the Teplá-Barrandian and Saxothuringian unit contact), which significantly differs from the other units. Low upper-crust velocities suggest sedimentary and volcanic filling of the rift, while high lower-crust velocities and weak or even missing Moho implies possible upper mantle updoming.


The analysis of unusual long-period (5 s) body waves, observed between the P- and S-wave onsets at near-regional stations in Greece for an event whose source duration was short (1 s), was performed in relation with the diploma (MSc) thesis of J. Vackář. Preliminary results reveal that this (PL) wave exists predominantly on the radial component of the velocity records. The dispersion curves were derived from real records by means of the frequency-time analysis (code SVAL of P. Kolínský). First attempts to interpret the PL waves and dispersion curves by means of synthetic seismograms were also made. Various effects, such as the source depth and/or anomalous Vp/Vs velocity ratios, were considered. This study will be finished in 2012.


Free oscilations of the Earth

Zábranová et al. (2012) analysed synthetic seismograms for several source fast solutions of the 2011 Tohoku earthquake obtained from surface waves and tested them against the observed gravity data from the superconducting gravimeter installed at the GOPE station. Good fit of the synthetic amplitude spectrum with the data up to 1.7 mHz was obtained without an additional increase of the moment magnitude. In this aspect, the 2011 Tohoku earthquake was similar to the 2010 Maule earthquake and not to the 2004 Sumatra-Andaman earthquake, where the free-oscillations studies resulted in an increase of the early Mw values.



Theory of seismic waves (reported by L. Klimeš)


Paraxial ray methods and Gaussian beams in anisotropic media

Cerveny and Psencik (2011c) summarized the main principles of seismic ray theory for heterogeneous isotropic and anisotropic media with curved interfaces, together with paraxial ray methods and some extensions and modifications of seismic ray theory in the invited review paper for the Encyclopedia of Solid Earth Geophysics.


Seismic waves in anisotropic attenuating media

A very suitable approximate method for calculating the complex-valued travel time in attenuating media is the perturbation from the reference travel time calculated along real-valued reference rays to the complex-valued travel time defined by the complex-valued Hamilton-Jacobi equation.  The accuracy of this approach strongly depends on the choice of the reference and perturbation Hamiltonian functions.  Klimes and Klimes (2011) proposed a novel general method for constructing the reference and perturbation Hamiltonian functions for any given complex-valued Hamiltonian function.


Paper by Cerveny and Psencik (2011a) is devoted to the computation of attenuation angles of inhomogeneous plane waves propagating in isotropic or anisotropic viscoelastic media and to the computation of their boundary values.


Cerveny and Psencik (2011b) investigated time-harmonic, homogeneous and inhomogeneous plane waves by the componental specification of their slowness vectors.  In the componental specification, the plane wave is specified, at an arbitrarily chosen plane, by a known projection of its slowness vector on the plane. Analogously to the mixed specification of slowness vectors of plane waves, the componental specification leads to the solution of an algebraic equation of the sixth degree with complex-valued coefficients.  Some special cases are discussed in a greater detail, in which the algebraic equation factorizes and allows analytical solution.  The componental specification plays an important role in the problem of reflection and transmission of plane waves on a plane interface between the viscoelastic halfspaces.


Anisotropic ray theory

Klimes (2011b) simplified the derivation of the zero-order ray-theory Green tensor in a heterogeneous anisotropic medium.


Coupling ray theory

Bulant, Psencik, Farra and Tessmer (2011a, 2011b) compared three different approximations of the coupling ray theory for S waves with the very accurate synthetic seismograms calculated by the Fourier pseudo-spectral method. The comparison is performed in six velocity models considerably differing in their anisotropy.


Ray-based Born approximation

Sachl (2011a, 2011b) tested the behaviour and accuracy of the ray-based first-order Born approximation for simple perturbations of a homogeneous background model, and of a heterogeneous background model. He also derived and applied the correction enabling to calculate 3-D seismograms by the 2-D Born approximation in a 2-D configuration.


Structural seismology

Bulant and Martakis (2011) smoothed given P-wave and S-wave propagation velocities in order to obtain a smooth 2-D velocity model suitable for ray tracing. They also smoothed given positions of structural interfaces and constructed a 2-D velocity model with interfaces by smoothing the same given P-wave and S-wave velocities separately in individual blocks. In the constructed model with interfaces, they calculated synthetic seismograms which simulate a reflection measurement in the modelled locality.


Klimes (2011a) derived that, for a given source, the migrated section is the convolution of the reflectivity function with the corresponding local resolution function, or the convolution of the spatial distribution of the weak-contrast displacement reflection-transmission coefficient with the corresponding local resolution function.


Bucha (2011) demonstrated the behaviour of the 3-D Kirchhoff prestack depth migration in the presence of triclinic anisotropy. He then studied the effects of simplified or incorrectly estimated anisotropy upon the migrated image.


CD-ROM with SW3D software, data and papers

Compact disk SW3D-CD-15 (Bucha and Bulant, 2011) 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-15 also contains over 400 complete papers from journals and from the SW3D consortium research reports, and 3 books by V. Cerveny. The software and papers from compact disk SW3D-CD-15 can be found at "".





(reported by Ondřej Čadek and Zdeněk Martinec)



In the field of planetology, attention has been paid to development of a fully self-consistent model of Saturn's moon Enceladus. Behounkova et al. presented a 3d spherical model of its thermal evolution in which episodes of enormous tidal heating are triggered by an increase of excentricity. The enhanced tidal heating leads to intense melting of the ice mantle, growth of the sub-surface ocean and destruction of the uppermost stiff lid, which further increases tidal dissipation. The research of terrestrial planets concentrated on methodical aspects, namely on including elastic behavior of the lithosphere in models of thermal convection. The analysis of the problem and first numerical results were presented in the paper by Golle et al., prepared in cooperation with our colleagues at the University of Nantes, France. Another study devoted to the same problem is the paper by Kalousova et al., which addresses the role of variable elastic lithosphere thickness. In the framework of her PhD project, Kalousova also studied the equations governing the behavior of a mixture of solid and liquid phases. The aim of this study is to better understand the thermal evolution of and mass transport in the ice shells on Europa and Titan. The project is supervised by Ondrej Cadek (Charles University) and Gael Choblet (University of Nantes) in cooperation with Ondrej Soucek from the Mathematical Institute of the Charles University.


Tectonophysical modeling

Maierova et al. completed the numerical model of Variscan evolution of the Bohemian Massive. The model suscessfully simulates three basic stages of the orogene evolution predicted on the basis of geological data and can also be generalized to modern orogenes where geological data on internal structure and conditions are missing.


Electromagnetic induction research

1. Vector potential formulation of a quasi-static EM induction problem: existence, uniqueness and stability of the weak solution.

We are dealing with the global-scale electromagnetic induction in a spherical conductor with 3-D distribution of electric conductivity in quasi-static approximation. We focus on theoretical aspects related to the solvability and posedness of this problem. We formulate the initial, boundary-value problem for a magnetic vector potential only, first in differential and then in integral forms. After cumbersome mathematical considerations, we prove that the problem of electromagnetic induction with 3-D distribution of electric conductivity which is formulated only in terms of magnetic vector potential, is well posed in Hadamard's sense. This means that a solution exists, is unique and continuously depends on boundary data, that is, a solution is stable. An important feature of the formulation is that no electric scalar potential is employed in the formulation and no gauge condition is imposed on the magnetic vector potential. The three attributes, that is the well posedness, only magnetic vector potential is employed and no gauge condition is imposed on a solution, make the formulation attractive for numerical implementation.



2. The modelling of toroidal magnetic field induced by tidal ocean circulation

We are dealing with the magnetic field induced by ocean circulation. The study is motivated by the fact that observations of the ocean-induced magnetic field by the CHAMP magnetic space mission, and most likely by the SWARM magnetic mission (its launch has been postponed to 2012) have the potential to be used as constraints on ocean dynamics. This has already initiated theoretical studies on the poloidal magnetic field induced by the horizontal ocean-circulation flow where the ocean layer is approximated by a single conductive sheet. Since the toroidal magnetic field cannot be modelled by this approximate model, we treat the ocean as a layer of finite thickness and model the toroidal magnetic field by a matrix-propagator technique with a source of electrical currents in the ocean layer. Although the primary toroidal magnetic field is not observable outside the oceans, it couples with a strong conductivity contrast between the oceans and continents and generates a secondary poloidal magnetic field which is observable by magnetic satellite missions and ground-based magnetic observatories situated close to the shoreline.



3. The adjoint sensitivity method of global electromagnetic induction for downward continuation of secular magnetic variations

Martinec (1999) developed a time-domain spectral-finite element approach for the forward modelling of vector electromagnetic induction data as measured on ground-based magnetic observatories or by the CHAMP satellite. We used the same time and space representations of magnetic induction vector and designed a new method of computing the sensitivity of the ground-based or CHAMP magnetic induction data to a magnetic field prescribed at the core-mantle boundary, which we term the adjoint sensitivity method. The forward and adjoint initial boundary-value problems, both solved in the time domain, are identical, except for the specification of prescribed boundary conditions. The respective boundary-value data on the ground or at the satellite‘s altitude are the X magnetic component measured by a vector magnetometer for the forward method and the difference between the measured and predicted Z magnetic component for the adjoint method. The squares of the differences in Z magnetic component summed up over the time of observation and all spatial positions of observations determine the misfit. Then the sensitivities of observed data, that is the partial derivatives of the misfit with respect to the parameters characterizing the magnetic field at the core-mantle boundary, are then obtained by the surface integral over the core-mantle boundary of the product of the adjoint solution multiplied by the time-dependent functions describing the time variability of magnetic field at the core-mantle boundary, and integrated over the time of observation. We tested and benchmarking numerical code for 1-D electric conductivity model of the Earth's mantle.


We have also contributed to the development of the Swarm Level 2 processing facility for the upcoming Swarm satellite geomagnetic mission commisioned by ESA. Working in close cooperation with ETH Zurich, the time-domain approach to recover the 3-D electrical conductivity of the Earth's mantle has reached maturity, passed initial tests, and shall be implemented as "Category 1" chain within the facility framework.


Jakub Velímský has also participated at the construction of a new geomagnetic observatory at Gan airport, Maldives, as part of the international effort to improve the geomagnetic data coverage of the southern hemisphere prior Swarm launch.



Ice sheat dynamics

1. On the long-term memory of the Greenland Ice Sheet

The memory of the Greenland Ice Sheet (GIS) with respect to its past states is analyzed. According to ice core reconstructions, the present day GIS reflects former climatic conditions dating back to at least 250 thousand years before the present (kyr BP). This fact must be considered when initializing an ice sheet model. The common initialization techniques are paleoclimatic simulations driven by atmospheric forcing inferred from ice core records and steady state simulations driven by the present day or past climatic conditions. When paleoclimatic simulations are used, the information about the past climatic conditions is partly reflected in the resulting present day state of the GIS. However, there are several important questions that need to be clarified. First, for how long does the model remember its initial state? Second, it is generally acknowledged that, prior to 100 kyr BP, the longest Greenland ice core record (GRIP) is distorted by ice flow irregularities. The question arises as to what extent do the uncertainties inherent in the GRIP based forcing influence the resulting GIS? Finally, how is the modeled thermodynamic state affected by the choice of initialization technique (paleo or steady state)? To answer these questions, a series of paleoclimatic and steady state simulations is carried out. We conclude that (1) the choice of an ice covered initial configuration shortens the initialization simulation time to 100 kyr, (2) the uncertainties in the GRIP based forcing affect present day modeled ice surface topographies and temperatures only slightly, and (3) the GIS forced by present day climatic conditions is overall warmer than that resulting from a paleoclimatic simulation.


2. ISMIP-HEINO experiment revisited: effect of higher-order approximation and sensitivity study

We revisit the results of the ISMIP-HEINO benchmark by first analyzing the differences in various model outputs using a wavelet-based spectral technique. Second, the ISMIP-HEINO benchmark experiments are recomputed with a novel numerical ice-sheet model based on the SIA-I algorithm that enables both the shallow-ice and a higher-order approximation of the ice-flow equations to be performed. To assess the significance of the higher-order approximation in the ISMIP-HEINO experiment, a numerical sensitivity study for the shallow-ice approximation (SIA) simulations is also carried out. A high sensitivity of the SIA model response to surface temperature perturbations is found. We conclude that the variations in ISMIP-HEINO results are due to the differences in (1) simulated basal temperatures and (2) numerical treatment of the basal sliding condition.


Earth’s lower mantle dynamics

1. Iron spin transition and the viscosity variations in the lower mantle

Matyska et al. (2011) parameterised variations of the activation energy caused by the Fe++ spin transition and incorporated them into the convection models. Using a parameterized model capturing the basic physics, they studied the dynamical consequences within the framework of a 2-D Cartesian convection model. Their models are characterized first by a low-viscosity channel below the 670 km boundary caused by an increase of temperature and decrease in grain size reduction with the post-spinel transition. Further, the variability of the rheological parameters leads to the prevalent formation of viscosity minimum at a depth of about 1600 km caused by spin crossovers followed by a ’viscosity hill’ in the bottom half of the lower mantle, the latter due to the increase of the activation energy in some low spin irons at the bottom of the mantle. Numerical simulations reveal a tendency to the formation of bigger stabilized plumes and/or plume clustering in the deep lower mantle. Increased vigor of convection in the lower mantle due to a potential increase of thermal conductivity and/or decrease of viscosity in the D”-layer due to the presence of ppv - results in appearance of small-scale convection right below the transition zone. These small-scale flow patterns in the lower mantle may influence the local mixing of geochemical anomalies and appearance of secondary plumes in the upper mantle.


2. Slab sinking speed and the viscosity of the lower mantle

Čížková et al. (2012) proposed a new method for inferring the mantle viscosity using the sinking speed of lithosphere subduction remnants. The sinking speed of paleoslabs recently reported by Van der Meer at al. (2010) was used as a constraint on the viscosity of the lower mantle. They performed a series of elaborate dynamic subduction calculations spanning a range of viscosity profiles from which we select profiles that predict the inferred sinking speed of 12 ± 3 mm/yr. Their modeling shows that sinking speed is very sensitive for lower mantle viscosity. Good predictions of sinking speed are obtained for nearly constant lower mantle viscosity of about 3-4x1022 Pas.  Viscosity profiles incorporating a viscosity maximum in the deep lower mantle, as proposed in numerous studies, only lead to a good prediction in combination with a weak post-perovskite layer at the bottom of the lower mantle, and only for a depth average viscosity of 5x1022 Pas.