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
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
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
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
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
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).
The above discussed source and ground-motion
studies in the
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)
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 26Al.
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. , 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.
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
models resulted in a nearly isothermal mantle wedge, where temperatures
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.
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.
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
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.