**Research overview 2011**

http://geo.mff.cuni.cz/documents/research_11.pdf

http://geo.mff.cuni.cz/documents/publ_11.htm

Department
of Geophysics, belonging to the Faculty of Mathematics and Physics,

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

Operation of all
these stations is partly supported by the project CzechGeo/EPOS, see http://czechgeo.cz . Within this project the
department operates also the CzechGeo/EPOS Seismological Software Centre, see http://epos-eu.cz/ssc.

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

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,

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 (

The
ISOLA software for multiple point-source moment tensor inversion of local and
regional seismograms, http://seismo.geology.upatras.gr/isola/, 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
"http://sw3d.cz".

**Geodynamics **

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

**Planetology**

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

*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-4x10^{22} 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 5x10^{22}
Pas.