Research overview 2008



Department of Geophysics, belonging to the Faculty of Mathematics and Physics, Charles University in Prague, has its roots seated as deep as in the 20's of the last century. Its structure and priorities have evolved from nearly pure seismology and geomagnetism to the present-day broad scope covering nearly all main branches of the physics of the Earth and, most recently, comprising also planetary aspects. Research is tightly coupled with education at the bachelor, master and doctoral levels. In 2008, there were 19 staff members at the department (counting permanent, temporary and part-time positions), and 20 PhD students (15 of them supervised by the staff members).


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 Corinth Rift Laboratory (F. Cornet). The Department also participated in the PASSEQ project, Passive Seismic Experiment in Trans-European Suture Zone (M. Wilde-Piorko) and the SWARM Science Study organized by the European Space Agency.


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 Prague in 1924), recently equipped with CMG-3ESP, was linked with the international data center ORFEUS to provide on-line data transfer to Utrecht. Four seismic stations in Greece have been operated by the department in cooperation with the Patras University, which started in 1997 ( Each of the stations SER, MAM, LTK and PYL is equipped with CMG-3T (broadband) and CMG-5T (strong motion) instruments. PYL and LTK, belonging to the satellite network PSLNET, are transmitted to the unified Greek network HUSN; LTK data is also transferred to ORFEUS. Equipment for a new station was purchased by our Department in 2008.


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


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á, Petra Maierová, Martina Ulvrová, Eliška Zábranová

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š (2008f).


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 (2008, in press 2008) are 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.


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-4.5 in Corinth Gulf, five events M>6 of 2008 (Leonidio, 2x Methoni, Andravida, Rhodos); then also M6 Parkfield 2004, M5.7 Gubbio, Italy 1984, two M<4 earthquakes in West Bohemia (Czech Republic), as well as hypothetical strong earthquakes M>6 near Basel and Los Angeles. The structural investigations concerned the Gulf of Corinth, Central Bohemia and simulated oil-field sites.


Source studies


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 Greece. For all of them, the identified fault plane was reported to EMSC, see Earthquake News and Highlights (“full list”):


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 ( Of about 40 earthquakes in Greece were processed by E. Sokos with ISOLA software at the University of Patras during 2008, using satellite broadband network PSLNET, and routinely reporting to EMSC. The same software was employed by Červinková (2008) in her diploma thesis to study the shallow M5.2 Trichonis 2007 earthquake in Greece. As shown in the thesis, the same focal mechanism can be obtained from 8 regional stations as well as for each of the stations individually. On the other hand, any attempt to apply a similar single-station inversion for the intermediate-depth M6.2 earthquake Leonidio 2008 has failed, probably due to missing Lg waves.


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 (> 1000 km). The near-regional records (< 400 km) can be explained without any non-DC component.


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.


Corinth Gulf is a typical rift zone, hence most expected focal mechanism is a normal fault. However, using weak events, hence going into smaller spatial scale, diversity of the mechanisms starts to appear, thus reflecting complexity of the tectonics. At the same time, for small events, severe methodical problems accompany the moment tensor calculation from the near-regional waveform inversion. The methodical aspects have been analyzed by Zahradnik et al. (2008a), who simulate the situation that four agencies determine the moment tensor, with slightly different methods, and obtain different results. Warning about several pitfalls is presented. They include, for example, erroneous mechanisms coming from inaccurate hypocenter location. Use of very near stations (<20km) is also discouraged, as the large amplitudes at very near stations might bias the inversion due to small errors in the location or small data problems. By data problem we mean clipping cutting a few waveform peaks, or a low-frequency instrumental disturbance, often not well apparent in band-pass waveforms, but having potentially a damaging effect on the inversion.


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 3 °C. The temperature variations represent a weak, but significant precursory phenomenon of the earthquake swarm.


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 Gubbio, Italy 1984 earthquake (Ameri et al., in press) was investigated by this technique, comparing results with independent method of F. Pacor. In Wang et al. (2008) the finite source model was implemented in conjunction with a 3D finite difference method to investigate jointly the source and site effect for a hypothetical strong (M7) earthquake at the Newport-Inglewood fault in the Los Angeles sedimentary basin. Again, an important effect of the nucleation point position (unknown during predictions) was emphasized. A new aspect is the source and site effect upon rotational components of strong motions, important for highway bridges and other long structures (Wang et al., submitted).


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 Basel (Faeh et al., 2008). The method itself, based on Alterman and Karal, often criticized for its purely algorithmic formulation, has received its exact mathematical description (Opršal et al., 2008 and Opršal et al., in press). The same paper also discusses various implementations of this useful technique. To encourage more applications, the method of I. Opršal has been parallelized (Caserta et al., submitted).


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 Corinth Gulf region have a tight link with the crustal structure investigations. Attention has been paid to the vicinity of the town of Aigion where a sequence of smaller earthquakes occurred in the first half of the year 2001. Arrival times from this sequence have been used to derive a local structural model of the Aigion region, composed of homogeneous layers (Novotný et al., 2008). In comparison with the present model for the western part of the Corinth Gulf, the new model is characterized by a higher vP/vS velocity ratio and by higher velocities to a depth of 7 km. The new model was derived with the aim to get more accurate locations of future events in the vicinity of the town of Aigion.


In the territory of the Czech Republic, shallow crustal structure has been studied along a refraction profile in the Central Bohemian Pluton (Málek et al., submitted). The profile passed approximately from the Orlík water reservoir to the town of Příbram. A 1-D velocity model was derived using the Wiechert-Herlotz method. The model will be used to locate the induced seismic events in the area of the Příbram-Háje underground gas storage, and weak earthquakes that sometimes occur in the Orlík reservoir region.


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 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. [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 1200 °C even in relatively shallow depths. Circulation in the wedge was obtained since the employed viscosity law resulted in creation of a zone of substantially reduced viscosity in shallow depths below the continental crust.


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 500 km, with respect to the case excluding the slab. The effect of the slab can exceed 0.3 mm/yr on short-term variations when the asthenosphere viscosity is 1019 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.