Research overview 2009
Department of Geophysics,
belonging to the Faculty of Mathematics and Physics,
In 2009, the Department participated in two EC projects: C2C and IMAGES. C2C (2007-2011), Crust to Core: the Fate of Subducted Material, is a pan-European Marie Curie Research Training Network specialized in the fate of subducting material. IMAGES (2005-2009), Induced Microseismics Applications from Global Earthquake Studies, aims at transfer of knowledge between seismologists and applied geophysicists (Schlumberger Cambridge Research), studying microearthquakes induced by drilling for oil. Preparations were also carried to join one new European project, QUEST, Quantitative Estimation of Earth’s Seismic Sources and Structure (coordinated by H. Igel).
Research project SW3D, Seismic Waves in Complex 3-D Structures, coordinated at the Department of Geophysics since 1993, has continued successfully in 2009. 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.
Seismic station PRA (created
The theses defended at the Department in 2009 included the following:
B.C.: Helena Munzarová, Libor Šachl
M.S.: Ľubica Valentová
The department was happy to host the long-term visitors Olivia Golle, Alexandra Guy, Nicola Tosi within the framework of the C2C project and Ernestina Lucca.
As in 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. Červený and
Paraxial ray methods and Gaussian beams in anisotropic media
The paper by Červený and Moser (2009) is devoted to paraxial ray methods in the computation of seismic wavefields propagating in elastic media. The basic procedure in paraxial ray methods consists in the solution of dynamic ray tracing equations. The authors derived the initial conditions for dynamic ray tracing equations for the three basic situations:
a) initial point situated at an initial curved surface,
b) initial point situated at a planar or nonplanar curved line,
c) isolated initial point (point source).
The derived expressions for initial conditions are very general, valid for homogeneous or inhomogeneous, isotropic or anisotropic, media. Along an initial surface and initial line, the travel time may be constant or variable. The results presented in the paper extend the possible applications of the paraxial ray methods. For example, they will play an important role also in perturbation methods for weakly dissipative media, in which the dynamic ray tracing is needed.
Klimeš (2009a) derived the relations between the propagator matrix of geodesic deviation (i.e., paraxial-ray propagator matrix) and the second-order derivatives of characteristic function (i.e., second-order derivatives of two-point travel time) in general coordinates. The equations are applicable to Finsler geometry, to Riemann geometry, and to their various applications like the general relativity or the high-frequency approximations of wave propagation.
The paper by Červený and Klimeš (2009) is devoted to the second derivatives of travel-time fields in anisotropic media, computed by dynamic ray tracing (DRT) along a reference ray. Two forms of DRT system have been broadly used in applications: the DRT system in global Cartesian coordinates and the DRT system in ray-centred coordinates. Simple transformation relations between the second derivatives of the travel-time field in global Cartesian coordinates and in ray-centred coordinates are derived. These transformation relations can be used in many applications, including the computation of complex-valued paraxial travel times which are necessary in the evaluation of Gaussian beams.
The paper by Červený and Pšenčík (2009a) is devoted to the theory of Gaussian beams concentrated close to rays of high-frequency seismic body waves propagating in inhomogeneous anisotropic layered structures. The basic role in the computation of Gaussian beams is played by a 2x2 complex-valued matrix of second derivatives of the travel time field along the reference ray. This matrix can be computed in several ways by dynamic ray tracing. Different possibilities are discussed in detail and possible simplifications are outlined.
Seismic waves in anisotropic attenuating media
Considerable attention has been devoted to perturbation methods from perfectly elastic media to anisotropic attenuating media. In the paper by Červený and Pšenčík (2009b), the complex-valued travel time and the complex-valued travel-time gradient of seismic body waves propagating in heterogeneous weakly dissipative anisotropic media are studied using the travel-time perturbation theory. The final expressions require simple quadratures along the real-valued reference ray traced in the reference perfectly elastic anisotropic medium. The results are applicable to any type of seismic source, including a point source.
Calculation of complex-valued travel time
In real space, the eikonal equation for complex-valued travel time represents the system of two Hamilton-Jacobi equations for the real and imaginary parts of the complex-valued travel time. Klimeš (2009b) suggested to solve these equations numerically using the wavefronts obtained by wavefront tracing.
The numerical solution of the eikonal equation for complex-valued travel time led Klimeš (2009c) to the idea of incorporating the amplitude into the complex-valued travel time by moving the second-order travel-time derivatives from the transport equation to the eikonal equation.
Gaussian packet prestack depth migration
Bucha (2009) continued the work on the Gaussian packet prestack depth migration, and tested the Gaussian packet prestack depth migration in two simple models composed of homogeneous blocks separated by planar or curved interfaces.
Bulant (2009) smoothed the existing 1-D regional velocity model of the Bohemian Massif for ray tracing. He minimized the Sobolev norm composed of the second-order velocity derivatives in order to prevent ray chaos characterized by the corresponding Lyapunov exponents.
The medium correlation function
is of principal importance in refraction travel-time tomographic inversion.
Travel times measured during refraction seismic experiments contain information
on the correlation function primarily in the horizontal directions, whereas
sonic-log measurements and vertical seismic profiling in deep wells contain
information on the correlation function in the vertical direction. Bulant and
Klimeš (2009) inform about the available data from sonic-log measurements in
the very deep KTB well, which could be used for estimating the correlation
function for the western
CD-ROM with SW3D software, data and papers
Compact disk SW3D-CD-13 (Bucha and Bulant, 2009) 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-13 also contains over 350 complete papers from journals and from the SW3D consortium research reports.
Earthquake and structural studies (reported by
Earthquake studies were focused
onto source processes and strong ground motions, including site effects. The
following earthquakes were most intensively investigated in 2009: the damaging
M6.3 earthquake of June 8,
A new strategy to
reveal the spatio-temporal evolution of the earthquake rupture process from
near-regional data, without assuming a constant rupture velocity was proposed
by Gallovič et al. (2009). The approach
is based on the conjugate gradient method, for which the required
waveform-misfit derivative with respect to slip on the fault is analytically
expressed. The derivative is given by the back-propagation of residual
seismograms towards the source; hence similarity with the adjoint tomographic
inversions. A good initial source approximation is necessary, being obtained from
hypocenter location and centroid-moment tensor solution (see ISOLA software
discussed below). The iterative approach then gradually reveals major
characteristics of the source process. The method was applied to investigate a
line-source model of the damaging Mw6.3 Movri Mountain, June 8, 2008,
The so-called H-C geometrical method for quick
identification of the fault plane, developed at the department in 2008, has
been applied to a M6.3 event south of
Software ISOLA for moment-tensor calculations from regional waveforms has been further upgraded (http://seismo.geology.upatras.gr/isola/). For example, an option of multiple crustal models was launched. The method is routinely applied by E. Sokos, University of Patras, for earthquakes well recorded in Western Greece by satellite network PSLNET, most M>3, and the results are reported to EMSC. Consultations were provided to many external users of the software, too.
The weak earthquakes of the
The study of M5 Amfilochia earthquake, 2002, performed in previous years, was finally published by Adamová et al.(2009). Large non-shear components reported by agencies were explained as an artifact due to complexity of the event.
A new method to measure the size of the non-elastic source volume, based on spectral magnitude of the P- and
S-waves, was suggested (Duda and Janský, 2009). The technique was successfully
applied also to the catastrophic earthquake of 12 May 2008,
Modern topic of the rotational components of seismic
motions was dealt by Brokešová et al. (2009); two instruments and a methodology
were developed and patented: the rotational sensor, the generator of rotational
waves, and technique to apply these instruments in the seismic exploration. The
instruments have been under further intensive testing, including, for example,
the swarm region of the
Wang et al. (2009) investigated combined source and basin effects upon the rotational seismic ground motions. The latter have a great practical significance for earthquake effects upon highway bridges and other long structures.
Determination of the depth of very shallow events, of the order of
Similar problems, complicated also by instrumentation and site effects,
may accompany also the very active region of the western
Strong ground motions
Many scenarios for engineering use were calculated by the previously
developed hybrid integral-composite method of F. Gallovič. For example, the
Zollo et al. (2009) used synthetic ground motion
scenarios calculated by F. Gallovič to test their Early Warning system in the
Irpinia region of
A 3D hybrid method to
simultaneously address the seismic-wave propagation and site effects, upgraded
with a new scheme, was better described from the point of view of its
mathematical formulation (Opršal et al., 2009). The main innovation of the
paper is the generalization of the boundary condition acting between the two
complementary sub-regions, describing the regional and local wave propagation,
respectively. The 3D code has been applied to several sites in
In the territory of the
Geodynamics (reported by C. Matyska)
In the two
papers, Tosi et al. (2009a,b) addressed the role of
lateral viscosity variations (LVV) in predicting the long-wavelength geoid. The
effect of LVV on the geoid has been a puzzling problem since early 90s. Until
recently, no really conclusive answer was available, mostly because the
numerical codes were not able to treat LVV with a sufficiently high resolution
and the prediction of the geoid was only limited to cases with rather small
lateral viscosity contrasts. Using a high-resolution spectral code developed by
Tosi and Martinec (2007), Tosi et al. (2009) investigated the sensitivity of
the geoid to highly viscous slabs in the mantle in 2D axisymmetric spherical geometry.
They found strikingly large impacts of LVV on the geoid and demonstrated a
paradoxical effect, namely that the amplitudes of the geoid above subduction
zones would be unrealistically large (
The recently discovered iron spin transition in major
mantle minerals at high pressures should exert dramatic influences on the
transport properties, such as the activation creep parameters in the deep
mantle. As a consequence of the spin transition, there is a significant
softening of the elastic parameters and an attendant strong reduction of the
activation energy in the creep law, leading to its non-monotonic behavior to
develop in the middle of the mantle. Matyska et al. (submitted) studied the
dynamical consequences of the activation energy reduction within the framework
of a 2-D Cartesian convection model. They used a parameterized model of the
rheological parameters capturing the basic physics and found that the
non-monotonicity in the rheological parameters led to the formation of
viscosity minimum at a depth of about
followed by a ``viscosity hill'' in the bottom half of the lower mantle. The numerical simulations revealed a tendency to the formation of superplumes and/or plume clustering in the deep lower mantle with the dynamical behavior, which was influenced by the changes of the activation parameters throughout the lower mantle and other material properties at the bottom of the mantle (decrease of thermal expansivity, viscosity reduction due to the perovskite-postperovskite phase change, radiative thermal conductivity).
The results of the modelling study about the dynamical consequenses of the low viscosity postperovskite in the lowermost mantle were accepted for publication (Čížková et al., in press).
Important though indirect information about the internal structure of Venus is provided by its topography and geoid. This information has been used in last decades to constrain the Venus mantle viscosity structure and its dynamic regime. Recently, Pauer et al. (2006) performed the geodynamic inversion of the venusian geoid and topography and proposed a group of best fitting viscosity profiles. As there is no information available about the Venus mantle density structure, they assumed a simplified 2D depth-independent approximation of the density anomalies. Benešová and Čížková (to be submitted) use their viscosity model together with several other simple viscosity stratifications (linear increase, constant, constant with a stiff lithosphere) as an input for the mantle convection code. They carry out the simulations of the mantle evolution in both 3D spherical shell and a 2D axisymmetric model and check, whether the character of the dynamic topography and the geoid represented by their power spectra fits the observed quantities. They estimate the effects of the thermally-dependent viscosity, internal heating, adiabatic and viscous heating and varying Rayleigh number. Though the geoid and topography do not vary too much among their considered viscosity profiles, they generally conclude, that the best fit between the observed and model quantities is reached for the viscosity profile by Pauer et al. (2006) with the upper mantle viscosity of 1021 Pa s. Thermally dependent viscosity improves the fit. The observed geoid and topography over several selected plumes successfully fit the model results reached for the preferred viscosity profile.
Zábranová et al. (2009) have studied equations and boundary
conditions describing tidal deformations and the gravitational potential of
prestressed elastic spherically symmetric bodies that were decomposed into a
series of boundary value problems for ordinary differential equations of the
second order and that were converted into a set of linear algebraic equations
with an almost block diagonal matrix by using pseudospectral schemes on Chebyshev
grid. Both accuracy and efficiency of the method were shown by evaluating the
tidal Love numbers for the Earth's model PREM. The method was applied to
evaluation of the tidal Love numbers for models of Mars and Venus. The Love
numbers of the two simplified Martian models - the former optimized to
cosmochemical data and the latter to the moment of inertia - were h2=0.172
(0.212) and k2=0.093 (0.113). For Venus, the value of k2=0.295, well-known from
the gravity-field analysis, was consistent with the results for the model with
the liquid-core radius of
In cooperation with researchers from the University of Nantes, France, M. Běhounková developed a new important tool to study convection in planets and icy satellites driven by the heat released due to tidal dissipation. She successfully combined the 3D spherical code of Oedipus, developed by G. Choblet, for simulating thermal convection, and the code of O. Čadek, designed for evaluating tidal dissipation in planetary and moon mantles with a general 3D viscosity structure. The new code, constructed and successfully tested by Běhounková et al. (2009), is unique since, for the first time, it allows the effect of the tides on mantle convection to be incorporated in a self-consistent manner. The code has a large number of applications in planetary research. It is now used to study thermal evolution of Earth-like exo-planets and it is to be soon applied to investigating processes in some icy satellites (Europa, Enceladus) as well.
EM induction modelling
Using the technique of adjoint modelling, we have inverted one year of CHAMP satellite data in terms of 1-D and 2-D axisymmetric conductivity models of the upper and mid-mantle (Martinec and Velímský, 2009). The results suggest slightly larger average conductivity of the northern hemisphere in the upper mantle.
Seven years of CHAMP satellite data were also inverted in terms of 1-D conductivity with emphasis on the deepest part of the Earth's mantle, the D'' layer (Velímský, in review). No significant conductivity increase was observed. This can be explained (as confirmed by 3-D modelling) as a lack of interconnection of the highly conductive post-perovskite phase in the longitudinal direction, i.e., in the dominant direction of polarization of the external currents.
The oceans play a special role in electromagnetic induction due to thein relatively high conductivity and the dynamo effect of ocean currents. The magnetic field by ocean circulation motion can be divided into toroidal and poloidal parts. The toroidal magnetic field is generated by electric currents closing in vertical planes and is estimated to reach 100 nT in amplitude. The much weaker poloidal field, with amplitudes up to 10 nT, results from electric currents closing horizontally. It has a significant vertical component and reaches remote land and satellite locations. Much attention has been given to the periodic magnetic signals of ocean flow which is driven by the lunar tides, but no attention has been devoted to the induced toroidal field. Dostal and Martinec (submitted) developed the matrix propagator technique to compute the toroidal magnetic field inside the Earth and oceans with the aim to generate the secondary poloidal field due to electrical conductivity inhomogeneities in the Earth‘s crust and lithosphere. The first estimate of its strength is up to 10 nT that may be detectable, summed up with the primary poloidal field by future SWARM satellite mission.
The 3-D time-domain method designed for inversion of Swarm multisatellite data was further developed. An extension study comissioned by the European Space Agency, "Level 2 products and performances for mantle studies with Swarm", was successfuly completed. We have joined the consortium of leading European institutions in the proposal to ESA to run the Swarm Level 2 processing facility.
Glacial isostatic adjustment
The improvement in understanding of dynamic processes in the Earth‘s mantle demands to consider a non-linear rheology of mantle material. Whereas this rheology is accepted in the studies of mantle convection, the need of a non-linear material behaviour in modelling of glacial isostatic adjustment (GIA) is still under discussion. Almost all the predictions of ongoing present-day processes induced by GIA are based on the assumption of a linear Maxwell viscoelastic rheology. To study the influence of non-linear rheology on the GIA-induced motion, Martinec and Klemann (submitted) implemented a non-linearly stress-dependent rheology in the spectral finite-element formulation of a viscoelastic self-gravitating sphere. The main effect of a non-linear rheology on the GIA-induced motion is for times when a surface ice-mass load is changing most rapidly because of large induced loading stresses.
Sasgen and Martinec (submitted) performed
a joint inversion of gravity fields from the Gravity Recovery and Climate
Experiment (GRACE) for glacial-isostatic adjustment over North America and
present-day ice-mass change in
The Gravity Recovery and Climate
Experiment satellite mission has allowed the resolution of temporal variations
in the Earth‘s gravity field to serve as a new observable for monitoring mass
changes in cryosphere. Sasgen and Martinec (submitted) analysed the GRACE time
series from August 2002 to August 2008 with regards to regional ice-mass
Modern modelling approaches to GIA are based on several techniques ranging from purely analytical formulations to fully numerical methods. Various European teams nowadays are independently working on the post-glacial rebound process in order to constrain the rheological profile of the mantle and the extent and chronology of the late-Pleistocene ice sheets which are prerequisites for the determination of the GIA contribution to geodetic observables. Martinec contributed to the benchmark study performed within the Working Group 4 of the ESF COST Action ES0701. Improved constraints on models of Glacial Isostatic Adjustment focuses on i) load Love numbers and relaxation spectra, ii) the deformation and gravity variations driven by surface loads characterized by simple geometry and time-history, and iii) the rotational fluctuations in response to glacial unloading.
Souček and Martinec have incorporated their newly-proposed numerical algorithm SIA-I (Souček and Martinec, 2008) into a large-scale thermo-mechanical ice-sheet model. The mechanical part of the problem, that is the momentum equation and rheology are solved by the SIA-I technique, which improves the standard Shallow-Ice Approximation in a series of iterative steps. In addition, the numerical model involves the evolution equation for the free surface and the heat-transport equation. Preliminary testing of the model was completed by performing two series of European Ice-Sheet Model INiTiative (EISMINT) benchmarks focused on (i) evaluating the effects of thermo-mechanical coupling and (ii) investigating the dynamics of the Greenland Ice-Sheet (GIS) subject to various climatologic scenarios. A simulation coupling the dynamic ice-sheet model with the GIA model by Martinec et al. (2000) is in preparation. Also a simulation evaluating the effects of higher-order (non-shallow) dynamics to surge occurrence and characteristics for the ISMIP-HEINO experiment is being set up.
Bohemian Massif dynamics
Until recently, the geodynamical research at the Department was only focused on studying global and large-scale problems. Fast development of numerical tools and computer power now allows to model also regional and local evolution of geological systems and, in this way, to test geological hypotheses and theories. The challenging aspect of this research lies in a great amount of data collected by the geological community in the past that can be used to constrain the numerical models. P. Maierová, a PhD student at the Department of Geophysics, modified the finite-element software Elmer in order to address various geological problems related to the evolution of the Bohemian Massif during the Variscan orogeny. In cooperation with geologists from other institutions (S. Ulrich, O. Lexa, K. Schulmann), she developed a model of the eastern part of the Bohemian Massif (collision of the Moldanubian domain with Brunia) and she also studied hypothetical lower crust diapirism in the southern part of the Bohemian Massif. The publication of the results is being prepared.