**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
in

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**

**L. Klimeš)**

**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.

**Structural
studies**

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**

**J. Zahradník)**

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,
earthquake in

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 ^{-5}m. The paper (Zahradník and
Plešinger) was accepted in BSSA and will appear in 2010.

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,

**Rotational
motions**

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.

**Location
**

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
M5.7

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

**Structural
studies**

In the territory of the

**Geodynamics ****(reported by C. Matyska)**

**Mantle dynamics**

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).

**Planetary research**

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 10^{21} 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
variability in

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.

**Ice-Sheet dynamics **

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.