**Research
overview 2008**

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

In 2008, the Department
participated in two EC projects: C2C and IMAGES. C2C (2007-2010), *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. We were invited to join also three new European
project proposals, QUEST, *Quantitative
Estimation of Earth’s Seismic Sources and Structure* (coordinated by H. Igel), INSIGHT, *In-situ
observatories for understanding earthquake generation* (M. Cocco) and CRLab3000, *The
3000m deep observatory in the * (F. Cornet). The Department also participated in the PASSEQ
project,

Research project SW3D, *Seismic Waves in Complex 3-D Structures*,
coordinated at the Department of Geophysics since 1993, has continued
successfully in 2008. The project has been supported by 8 companies or research
institutes (BP America Inc., U.S.A.; Chevron U.S.A. Inc., U.S.A.; NORSAR,
Norway; Observatorio Nacional,
Brazil; OMV Exploration & Production GmbH, Austria; Petrobras,
Brazil; Rock Solid Images, U.S.A.; Schlumberger Cambridge Research Limited,
U.K.) in the framework of the international SW3D consortium.

Seismic station PRA (created
in

The 2008 was rich in the
number of theses defended at the Department in all three programs, bachelor,
master and doctoral (the theses can be found at
http://geo.mff.cuni.cz/prace.htm):

Bc.: Adela Androvičová, Martin Čochner, Soňa Formánková, Klára Kalousová, Miroslav Kuchta, Radim Kusák, Pavla
Procházková, Jan Vrba

Mgr.: Dana Červinková,

Ph.D.: Martin Kukačka, Nicola Tosi

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

**Theory of seismic waves ****(reported by V. Červený and**

**L. Klimeš)**

Theoretical investigation of
seismic wave propagation in complex structures and of seismic data inversion is
mostly being performed within the scope of the consortium project "Seismic
Waves in Complex 3-D Structures". The research is focused primarily on the
fundamental issues of high-frequency seismic wave propagation in complex 3-D
isotropic and anisotropic structures, which go beyond the traditional
approaches. The emphasis is put on new, stable, more efficient and flexible
algorithms for both forward numerical modeling and inversion of seismic wave
fields in 3-D inhomogeneous, isotropic or anisotropic, elastic or attenuating
structures.

**Stochastic
inversion of travel times**

The paper by Klimeš
(2008a), based on research conducted in period 1996-2002, lists the algorithms
we need for numerical realization of the stochastic inversion of travel times.
The results of the inversion are represented by the velocity model and the
model covariance function which specifies the deviation of the velocity model
from the geological structure. Whereas the velocity model is composed of smooth
functions of 3 coordinates, the model covariance function is a function of 6
coordinates with pronounced singularities. The computer representation of the
model covariance function in terms of two matrices and the rays used for
inversion has been designed.

The extensive paper by Klimeš (2008b) presents new equations for calculation of
the geometrical covariances of travel times in a
self-affine random medium. This calculation represents one of the cardinal
steps in the stochastic inversion of travel times.

**Inversion of
seismograms from a reflection seismic survey**

In 2007, it was determined how
the perturbations of a generally heterogeneous isotropic or anisotropic
structure manifest themselves in the wavefield, and
which perturbations can be detected within a limited aperture and a limited
frequency band. Perturbations of elastic moduli and
density can be decomposed into Gabor functions. The wavefield scattered by the perturbations is then composed
of waves scattered by individual Gabor functions. The
derived approximate solution of the forward scattering problem can be utilized
in designing the algorithm of the linearized
inversion of the complete set of seismograms recorded for all shots during a
reflection seismic survey. This challenge of replacing migration by true
inversion has represented the main research topic in 2008.

The paper by Klimeš
(2008c) describes the whole inversion method. Other 4 papers present the
solutions of individual problems related to individual steps in developing the
complete inversion algorithm.

The paper by Klimeš
(2008d) contains the equations required for coding the calculation of the
scalar products of the Gabor functions.

The paper by Klimeš
(2008e) presents a new theoretical result: the approximate analytic expressions
for the frame bounds corresponding to the standard discrete Gabor
transform. These expressions can be used to analytically study both the discretization error of the continuous Gabor
transform and the stability of the discrete Gabor
transform in dependence on the phase-space lattice intervals. Since this
theoretical result cannot be practically applied to the designed wavefield inversion, its considerable but unproved
generalization has been conjectured in the paper by Klimeš
(

The paper by Klimeš
(2008g) represents a simple attempt to specify the finite set of structural Gabor functions for the wavefield
inversion. Specification of the finite set of structural Gabor
functions constitutes the first step in the wavefield
inversion.

**Gaussian
packet prestack depth migration**

Bucha (2008) has continued the work
on the Gaussian packet prestack depth migration
commenced by Karel Žáček in
the period 1999-2005. The algorithm of smoothing the parameters describing the
shape of optimized Gaussian packets has been modified. The code for
optimization has been generalized from the zero offset to the common shot. The
code for decomposition of recorded seismograms into Gaussian packets and the
code for migration have considerably been sped up. Václav
Bucha then proceeded to the migration with optimized
Gaussian packets.

**Seismic waves
in anisotropic attenuating media**

Considerable attention has been
devoted to homogeneous and inhomogeneous waves propagating in anisotropic
attenuating media, especially to perturbation methods from perfectly elastic
media to anisotropic attenuating media.

In the paper by Červený, Klimeš & Pšenčík (2008), the attenuation vector of seismic body waves
propagating in heterogeneous weakly dissipative anisotropic media is studied
using the travel-time perturbation theory proposed by Klimeš
in 2002. The final expressions for the attenuation vector require simple quadratures along the reference ray, constructed in a
relevant reference perfectly elastic anisotropic medium. The results are
applicable to any type of seismic source, including a point source. It is shown
that the waves in heterogeneous weakly dissipative anisotropic media are, in
general, inhomogeneous. The exception are waves
generated by a point source in a homogeneous isotropic medium, which are always
homogeneous.

In the paper by Červený & Pšenčík (2008b),
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. This paper differs from the related paper by Červený,
Klimeš & Pšenčík (2008)
in the first place by considerably more general forms of the perturbation
Hamiltonian, which allow for more accurate perturbation expansions. Attention
has also been devoted to the special case of an isotropic reference medium.

The papers by Červený & Pšenčík (2008a,
2008d) are devoted to the exact and approximate expressions for the quality
factor Q of time-harmonic, homogeneous and inhomogeneous plane waves
propagating in unbounded homogeneous dissipative anisotropic medium. The
direction-dependent quality factor Q is defined as the ratio of the time
averaged complete stored energy and the dissipative energy, per unit volume. For
weakly inhomogeneous plane waves, propagating in weakly dissipative media, the
exact expression for Q can be simplified to the approximate one using the
perturbation method. Relation of the approximate Q to the attenuation
coefficient, measured along an arbitrary straight-line profile, has been
discussed. It has been shown that the approximate Q and the attenuation
coefficient do not depend on the inhomogeneity of the
plane wave if the profile is parallel to the ray-velocity vector of the plane
wave under consideration. Consequently, the approximate quality factor is a
convenient measure of the intrinsic dissipative properties of rock in the ray
direction.

In the paper by Červený & Pšenčík (2008c),
weakly inhomogeneous plane waves propagating in weakly dissipative anisotropic
media are discussed. The first-order perturbations are used to derive simple
expressions for the propagation and attenuation vectors and for the
polarization vectors. It is shown that the scalar product of the attenuation
vector with the energy-velocity vector represents an intrinsic attenuation
factor, which does not depend on the inhomogeneity of
the plane wave under consideration. Actually, the intrinsic attenuation factor represents
1/Q, where Q is the generalization of the quality factor for anisotropic media.
The accuracy of the perturbation method is studied by comparison of approximate
and exact results.

**Anisotropic
ray theory**

The papers by Červený & Moser (

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.

**Coupling ray
theory**

The paper by Bulant
& Klimeš (2008a) is the reprinted version of the
paper incorrectly printed in 2007. The authors numerically tested the equations
for the anisotropic-common-ray

approximation of the coupling ray theory in
three smooth 1-D velocity models of differing degree of anisotropy. They
calculated the synthetic seismograms obtained by both the isotropic-common-ray
and anisotropic-common-ray approximations of the coupling ray theory, and
compared these seismograms with the more accurate coupling-ray-theory synthetic
seismograms simulated by the second-order perturbation expansion of travel time
from the anisotropic common rays.

**Velocity
macro models**

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

**CD-ROM with
SW3D software, data and papers**

Compact disk SW3D-CD-12 (Bucha & Bulant, 2008)
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-12 also contains over 330 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 in Greece were investigated: M5 Amfilochia
2002, M5.5 Zakynthos 2006, M5.2 Trichonis
2007, six events M3-

A geometrical method for quick identification of
fault plane has been published (Zahradník et al.,
2008). The method is based on the strike and dip of the moment tensor solution,
and the relative position of two points, hypocenter (H) and centroid
(C), hence also H-C method. The application included M6.2 Leonidio
earthquake, January 6, 2008. The fault plane has been verified by the Coulomb
stress change, calculated for the previously published regional stress filed of
the region (after A. Kiratzi and C.B. Papazachos).

The H-C method has been also used during the year
for the other four significant events (M6) in

One of the main possible applications of the quick
fault identification is the Coulomb stress analysis, to assess places of
increased aftershock probabilities after strong events. The problem of stress
loading was also addressed in a theoretical study by Gallovič
(2008). Using ‘rate and state’ model of fault rheology,
regular repetition of events on the fault plane can be modeled. The paper
answers the question how the 3D fault behavior is modified under homogeneous or
heterogeneous Coulomb stress load. Advanced or even delayed earthquakes might
appear on the fault. Moreover, if the earthquake is triggered relatively soon
after the loading, position of the nucleation point seems to be spatially
related to the loaded zone. Knowledge of the nucleation point (or region) plays
an important role in the calculations of seismic scenarios as well as in the
time-dependent hazard analysis due to aftershocks (Gallovič,
Brokešová, 2008a,b).

Method ISOLA for moment-tensor calculations from
regional waveforms has been published (Sokos, Zahradník, 2008) and new version of the fortran-matlab
software released on web (http://seismo.geology.upatras.gr/isola/). Of about 40
earthquakes in

The ISOLA software was also used for two largest
(M>3.5) earthquakes of the 2008 West-Bohemia swarm, the solution also
reported to EMSC. The study was interesting mainly because of a substantial
contribution of near-field terms in the waveforms at the nearest stations, as
well as due to presence of long-period disturbances at the Novy
Kostel (NKC) station of the Czech Regional
Seismological Network. The study will continue in 2009 under cooperation with
A. Plešinger.

Two papers were devoted to
non-shear components. Zahradník et al. (2008b)
studied three main events (M5.5) of the sequence near Zakynthos,
2006, one of which had double-couple part as low as 20% only. An equivalent
interpretation in terms of a double event with significantly different
mechanism was suggested. The M5 Amfilochia earthquake, 2002, was analyzed by Adamová et
al. (in press). In the latter case, most likely, the large non-shear component
reported by agencies was an artifact due to time shifts introduced in the
agency solutions to compensate the unknown crustal
structure at long epicentral distances (>

An interesting attempt to define a seismic volume using estimates of the
highly deformed zone (>1e-4) was made by
Duda and Janský (2008); see
also Duda and Janský
(submitted). The size of the non-elastic region depends on period and wave
type, hence the difference between P and S magnitude is a key to decipher
the source size.

During the strong
earthquake swarm in the West-Bohemia/Vogtland region
at the turn of the years 1985 and 1986, significant co-seismic changes were
observed in the discharge of many mineral springs at the nearby Františkovy Lázně spa. These
changes were studied soon after the earthquake swarm. Some discharge changes
were very distinct, amounting up to 40 %. Moreover, recent detailed
analyses of the springs parameters have also revealed
certain variations in the temperature of some mineral springs (Stejskal et al., 2008). The temperature variations were
less noticeable than the discharge ones, but preceded the beginning of the
swarm by several months. The highest temperature increase was about

**Strong-ground
motion studies**

Many scenarios for engineering use were calculated by the previously
developed hybrid integral-composite method of F. Gallovič.
For example, the M5.7

A 3D method to simultaneously
address the propagation and site effects, continuously developed by I. Opršal, found its interesting application in explaining
archeological data due to historical earthquake near

A puzzling feature of the strong ground motions is
the effect of radiation pattern, by some authors considered
unimportant or even non-existent at high frequencies, say > 1 Hz
(e.g. A. Pitarka). Indeed, if using radiation pattern
in the codes to model near-fault records we face large amplitude ratios between
the three geometrical components of the motions. On the other hand, data
exhibit all three components rather similar to each other as far as their peak
amplitudes are concerned. The working hypothesis is that vanishing radiation
pattern at > 1 Hz is just a phenomenological effect. The source posses its
radiation pattern, but it is canceled, or reduced, by other effects, such as
the fault complexity, e.g. non-planar geometry, 3D heterogeneity of the crust,
including the fault region itself, surface topography, etc. Kaeser
and Gallovič (2008) investigated these effects. For a
specific event, Parkfiled 2004, data and modeling
indicate that just the 3D heterogeneity of the crust (partly reflecting also
site effects) is a very important candidate to explain the puzzling feature (Gallovič et al., submitted).

**Structural
studies**

The above discussed source and ground-motion
studies in the *v _{P}/v_{S}* velocity ratio and by higher
velocities to a depth of

In the territory of the

Feasibility to construct a 1D structural model by
means of the Neighbourhood algorithm, and jointly
locate weak events, was studied in synthetic tests simulating real conditions
of microearthquakes induced by hydrofracturing
at oil fields (Janský et al., submitted). The
attention was focused on receivers in wells. Two wells are needed to avoid
trade-off between location and structural model.

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

**Planetary research**

In this field attention has been
paid to the effects of tidal deformation on rotational dynamics and mantle
convection in icy satellites (Tobie et al., 2008; Robuchon et al., submitted) and terrestrial planets. The
research presented in the paper by Tobie et al.
(2008) concentrated on Saturn’s tiny moon Enceladus
which is one of four solar system solid objects that are sufficiently
geologically active for their internal heat to be detected by remote sensing.
The endogenic activity on Enceladus
is only located on a specific region at the south pole.
Using a model of tidal deformation of a viscoelastic
body, Tobie and co-authors demonstrated that only
interior models with a liquid water layer at the base of the icy mantle can
explain the observed magnitude of the heat production and its particular
location at the south pole. Another moon of Saturn, Iapetus, was investigated by Robuchon
et al. (submitted). The Cassini mission revealed that
Iapetus (i) is geologically
old and has a high equatorial ridge, and (ii) has a large flattening which does
not correspond to its very slow present rotation. The authors studied the
latter phenomenon by simulating the evolution of the interior thermal structure
and its spin rate and global shape for a simplified viscoelastic
model. They showed that the present shape of the satellite provides a strong
constraint on the concentration in the short-lived radiogenic elements, namely ^{26}Al.

Important though indirect
information about the structure of Venus provides the topography and the venoid. We have developed a numerical code to simulate the
thermal convection in the mantle of the planets and carried out several
simulations of the Venus mantle convection. We concentrated on the question how
the geoid spectra computed from density distribution
obtained in our convection models correspond with those of the observed data.
We found that the models with the viscosity increasing with the depth produce
relatively good fit with the observed geoid spectra.
The thermally dependent viscosity and extended Boussinesq
approximation improve the fit. On the other hand, in the model with a constant
viscosity the predicted spectra deviate from the observed one considerably,
especially its wavelengths part.

Using synthetic gravity data, Pauer et al. attempted to answer the question of the detectability of oceanic floor topography on Europa (Pauer et al., submitted).
Gravity field data were also used to estimate crustal
density on Mars (Pauer et al., 2008).

Šrámek and collaborators in France
studied separation of metal from the silicates leading to core formation in
growing proto-planets. They developed a multiphase numerical code that
self-consistently accounts for the presence of solid silicates, and solid and
liquid iron (Ricard et al., submitted; Šrámek et al., to be submitted). At each point an average
velocity describes the incompressible flow due to the lateral density
variations in the mixture, while a segregation (irrotational)
velocity field allows the sinking of the denser molten metal across the
silicate matrix. The energy balance accounts for changes in potential energy
associated with segregation. It is shown that the core formation starts before
a significant melting of the silicates, as soon as impact heating is large
enough to reach the melting temperature of the metallic component. The first
metallic diapirs leave behind a trailing conduit. The
segregation proceeds in a runaway process as the gravitational energy released
as heat upon segregation brings up the temperature
above the metal melting temperature in adjacent regions. This results in the
whole core-mantle segregation and inevitably occurs for undifferentiated
embryos of Moon to Mars size.

**Dynamics of
the lower mantle**

Recent evidence on perovskite to post-perovskite
phase change in the lowermost mantle suggests that post-perovskite
piles or lens should be present in the relatively cold downwelling
areas, while the roots of the hot upwelling plumes consist of perovskite. Post-perovskite is
generally believed to be deformed predominantly by the dislocation creep and
there are some indications that the activation parameters of dislocation creep
in post-perovskite induce lower viscosity than the
surrounding perovskite, which may have important
consequences for the slab deformation in the D''. Čížková et al.
(to be submitted) performed the simulations of
the thermal convection in a 2D cartesian model of the
lower mantle with composite rheology including
diffusion creep and dislocation creep. Different creep parameters are used for
the perovskite lower mantle and for post-perovskite lens (or layer) respectively. The presence of
the rheologically distinct post-perovskite
strongly influences not only the slab deformation above the CMB but results in
different dynamic regimes of the CMB region. Further, the presence of a very
low viscosity lens or layer above the core-mantle boundary has a strong
influence to the heat flux from the core.

The mechanical properties of lithospheric
slabs in the lower mantle are still a controversial issue. Since the
temperature of subducting slabs is lower than the
average temperature at a given depth, it is reasonable to assume that the
effective viscosity of slabs in the lower mantle remains significantly higher
than that of the surrounding mantle. Using a simply parameterized mechanical
model, Tosi et al. (to be submitted) investigated the
effect of viscous slab thickened in the lower mantle on the gravitational
response of the Earth. They found that the lower mantle slabs generate too big geoid signal (~200 m) if their viscosity is significantly
(more than 10 times) higher than the mean viscosity of the lower mantle. Hence,
to obtain reasonable amplitudes of the geoid, the
slabs must have either similar viscosity as the surrounding mantle or the same
thickness as in the upper mantle. However, recent tomographic
studies have indicated that the shape of slabs changes as they pass through
upper/lower mantle boundary and the subducted
lithosphere has a blob-like rather than plate-like structure in the lower
mantle.

As shown by Tosi et al.
[2008], another way to decrease amplitudes of the geoid
consists in incorporating a low-viscosity postperovskite
phase, presumably existing in the lowermost mantle, into numerical models.
Although the viscosity of postperovskite is unknown
and extremely difficult to measure, there are indications (mostly based on
analogy between electrical conductivity and diffusion creep viscosity) that it
could be significantly lower than that of perovskite.
If this is the case, reasonable amplitudes of the geoid
can be obtained even for a relatively high slab viscosities
above the D” layer.

**Subduction**** modeling**

Kukačka and Matyska (2008) modeled
mantle wedge flow driven by subducting plates and its
consequences to surface heat flow in back-arc regions, where
an increase in surface heat flow
has been observed in almost all subduction
zones. They showed that such a robust behavior can be reproduced
if a strong pressure-dependence of viscosity is taken
into account together with its
temperature-dependence. The
models resulted in a nearly isothermal mantle wedge, where temperatures
exceeded

Corner flow and
temperature in the mantle wedge caused by viscous coupling between the subducting
slab and the overriding mantle wedge was studied by Kukačka and Matyska (submitted). They developed a numerical model of a generic pseudo-plastic subduction plate, in which the steady-state thermal structure
of the mantle wedge is calculated
without any *a priori *assumptions on the geometry of the subducting
slab. The results showed that both
the temperature field and the
geometry of the slab are strongly controlled by weakness zones in the topmost part of the subducting
plate and in the cold tip of the
mantle wedge. They showed that a channel
of backward flow may evolve
along the plate-wedge interface if the cold portion
of the mantle wedge is weak,
which likely is true as this
region is hypothesized to be strongly serpentinized.

**Convection
modeling**

Van Hunen
and Čadek (submitted) studied the effect
of a small-scale convection on the magnitude and orientation of anisotropy in
the upper mantle. Using a fully 3d convection model, they demonstrated that the
strikingly small amplitudes of seismic anisotropy in comparison with a
theoretical prediction result from an interaction between large-scale plate
motion, which tends to establish anisotropy parallel to the plate motion, and
the small-scale convection, which tends to destroy a uniform anisotropy
pattern.

**EM induction
modeling**

Velímský et al. (submitted) have
developed and tested the adjoint method to the global
3-D EM induction problem. This allows efficient evaluation of misfit gradient
in the model space and thus makes feasible solution of 3-D inverse problem by
gradient methods. This approach was also reported as part of the SWARM Science
Study comissioned by the European Space Agency and is
considered to become part of the Level 2 processing chain for the upcoming
SWARM multisatellite mission. Similar approach was
applied to the CHAMP satellite data in a simplified 2-D axisymmetric
case (Martinec and Velímský,
submitted). 1-D conductivity model was recovered and sensitivity to possible
conductivity differences between northern and southern hemisphere was studied.

**Ice-sheet
modeling**

Souček and Martinec
(2008) have continued testing their algorithm for a fast iterative improvement
of the extensively used Shallow-Ice Approximation. Their algorithm enables for
a wide range of realistic conditions to attain the full-Stokes solution in
significantly shorter times (up to 2 orders of magnitude) compared to standard
modeling approaches such as based on the finite-element method. The
computational performance and accuracy of the algorithm were also successfully
tested in the ISMIP-HOM international benchmark. The algorithm was implemented
to a new numerical code for ice-sheet evolution, a
novel approach in glaciology was applied to track the boundary of a glacier by
means of the level-set method.

**Glacial isostasy and plate motion **

The influence of glacial-isostatic adjustment (GIA) on the motion of tectonic plates
is usually neglected. Employing recently developed numerical approach, the
effect of glacial loading on the motion of the Earth's tectonic plates has been
examined by Klemann et al. (2008), where an elastic
lithosphere of laterally variable strength and the plates losely
connected by low viscous zones has been considered. The aim was to elucidate
the physical processes which control the GIA-induced horizontal motion and to
assess the impact of finite plate-boundary zones. The present-day motion of
tectonic plates induced by GIA was obtained at, or above, the order of accuracy
of the plate motions determined by very precise GPS observations. Therefore,
its contribution should be considered when interpreting the mechanism
controlling plate motion.

**Postseismic**** deformation
in a viscoelastic self-gravitating spherical Earth **

Theoretical models of the viscoelastic relaxation of a spherical earth are derived to
model large-scale post-seismic deformation resulting from great earthquakes
(M>7) over decadal time scales. Most existing models of post-seismic
deformation do not consider strong lateral heterogeneities in mantle viscosity,
in particular in the subducting slab, where such
events occur. In addition, the self-gravitation effect is often treated only
approximately. Both effects become important when observations from space
geodetic techniques such as GPS and GRACE are interpreted. Tanaka et al. (2008)
has used a spectral finite-element approach, that
allows these two effects to be considered in a rigorous way. In this way, much
larger lateral viscosity variations can be handled than by perturbation
techniques. Interface conditions for an arbitrary shear fault in the form of
double-couple forces that are equivalent to a prescribed dislocation, have been
derived, and a relaxation process for an incompressible Maxwell earth with a
three-dimensional viscoelastic structure has been
simulated. As an example, this approach has been applied to the 2004
Sumatra-Andaman earthquake, and a large-scale post-seismic gravity potential
variation has been simulated by a forward calculation. In the presence of a
slab, a secular variation in the geoid height change
decreases by 30% for wavelengths longer than ^{19} Pa s, which are
larger than observation errors of GRACE. For a displacement field, a relative
decrease in deformation rates can amount to 70% due to the inclusion of a slab,
which is detectable with geodetic observations such as GPS. The effect of the
slab is attenuated in the gravity field for such longer wavelengths, since
horizontal scales of the slab are smaller than its spatial resolution. Lateral
heterogeneities in viscosity due to a slab should be considered for large
thrust events, in order to reveal rheology in a subduction zone through comparisons between modeled and
observed post-seismic deformations.

**Semi-analytical
solution for viscoelastic relaxation in eccentrically
nested spheres induced by surface toroidal traction**

A semi-analytical solution to
the 2-D forward modeling of viscoelastic relaxation
in a heterogeneous sphere induced by a surface toroidal
force has been derived by Martinec (2008). The model
consists of a concentrically nested elastic lithosphere, a viscoelastic
mantle, and an eccentrically nested viscoelastic
core. Since numerical codes based on finite-element or
spectral-finite-difference techniques for modeling viscoelastic
relaxation in a spherical geometry in the presence of lateral viscosity
variations are becoming more popular, reliable examples for testing and
validating such codes are essential. The eccentrically-nested sphere solution
has been tested by comparing it with two distinct results: the analytical
solution for viscoelastic relaxation in
concentrically nested spheres and the time-domain, spectral finite-element
numerical solution for viscoelastic relaxation in
eccentrically nested spheres, with excellent agreement being obtained.