Research overview 2006

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 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 2006, there were 19 staff members at the department (counting permanent and temporary positions), and 17 PhD students (14 of them supervised by the staff members).

In 2006, the Department of Geophysics participated in three EC projects: SPICE, 3HAZ-CORINTH and IMAGES. SPICE (2004-2007) is the pan-European Marie Curie Research Training Network involving 14 universities and specialized in theory of seismic wave propagation. 3HAZ-CORINTH (2004-2006) is a research project specifically targeted on three main natural hazards in the Gulf of Corinth, Greece, viz earthquakes, landslides and tsunamis. IMAGES (2005-2008) aims at transfer of knowledge between seismologists and applied geophysicists (Schlumberger Cambridge Research), studying microearthquakes induced by oil drilling.

An important contribution related to our participation in the SPICE project was completion of a tutorial book by J. Brokesova on the asymptotic rayámethod in seismology. Publication of the book was sponsored by SPICE, copies are available on request.

Research project " Seismic Waves in Complex 3-D Structures" (SW3D), coordinated at the Department of Geophysics since 1993, has continued successfully in 2006. The project has been supported by 6 oil companies (BP America Inc., U.S.A.; Chevron U.S.A. Inc., U.S.A.; ExxonMobil Upstream Research Company, U.S.A.; Petrobras - CENPES, Brazil; Schlumberger Cambridge Research Limited, U.K.; Shell International Exploration and Production B.V., Netherlands) in the framework of the international SW3D consortium.

Apart from the seismic station PRA (created in Prague in 1924), four seismic stations in Greece have been operated by the department in cooperation with the Patras University since 1997 ( http://seis30.karlov.mff.cuni.cz).

Four Bc. theses were defended at the department in 2006 (D. Cervinkova, P. Maierovß, M. Ulvrova and E. Zabranova), two MSc theses (P. Adamova, N. Benesova) and one PhD thesis (F. Gallovic).

Similarly to the previous years, research at the department was carried out in three directions: Geodynamics, Theory of seismic waves, and Earthquake and structural studies

Geodynamics (reported by O. Cadek)

Thermal convection in a mantle with post-perovskite phase change

Matyska and Yuen (2006) investigated the dynamical impact arising from the interaction of temperature-dependent and depth-dependent viscosity with radiative thermal conductivity on ascending and descending flows in the presence of an endothermic phase change at 670-km depth and an exothermic post-perovskite transition at 2650 km. Their results for temperature- and depth-dependent viscosity suggest that the radiative thermal conductivity is important for sustaining superplumes in the lower mantle, once the post-perovskite phase change is brought into play. This aspect is especially emphasised, when the radiative thermal conductivity is restricted only to the post-perovskite phase. In this case, mass and heat-transfer between the upper- and lower mantle deviates substantially from the traditional whole-mantle convection model since an overall complete communication between the top and lower mantle is difficult to be achieved.

Viscosity structure of the lowermost mantle

Cadek and Fleitout (2006) studied the lateral distribution of the viscosity anomalies in the core-mantle boundary (CMB) region and its relationship to the long-wavelength gravitational field and the structure information provided by seismic tomography. By inverting the geoid and free-air gravity data at harmonic degrees 2-8, they derived a large-scale model of lateral viscosity variations in the lowermost mantle. The viscosity pattern found for the CMB region shows a high density of hotspots above the regions of higher-than-average viscosity. This result suggests an important role for petrological heterogeneity, potentially associated with a post-perovskite phase transition.

Conductivity of the Earth from satellite measurements

A time-domain approach to the global electromagnetic induction problem was applied to vector magnetometer data observed by the CHAMP satellite (Velimsky et al., 2006). Data recorded during 11 storm events in 2001-2003 were processed track by track, yielding time series of spherical harmonic coefficients. The data were interpreted in terms of 1-D layered electrical conductivity models. In the upper 50 km the inversion solidly recovered a conductive layer corresponding to averaged surface conductance. The conductivity of the lower mantle was established at 6 S/m assuming the upper-lower mantle interface at the seismic-based 670-km boundary. The resolution of the method in the resistive upper mantle sandwiched between conductive crust and lower mantle was poor. Nevertheless, an upper bound of 0.01 S/m was suggested by the data. A conductivity increase in the transition zone was not observed. Subsequently, J. Velimsky modified the time-domain approach to EM induction using combined boundary conditions prescribed at the Earth's surface and satellite altitude. Preliminary inversions of 3-year time series of observatory and satellite data yield less conductive lower mantle than previous studies.

Ice sheet modelling

Soucek and Martinec (2006) continued investigating the thermodynamic formulation for polythermal ice on the basis of rational thermodynamics of reacting mixtures and have finalised the shape of the governing equations for polythermal regions, such as a binary mixture of ice and water. The boundary and transition conditions were derived consistently from the rational thermodynamic framework. The authors derived the shallow ice scaling approximation for a general setting of an ice sheet in spherical geometry and investigated solution techniques for numerical implementation of the theory.

Earth gravitational field from satellite data

Together with his co-workers from GFZ in Potsdam, Germany, Z. Martinec presented a spatial averaging method for Gravity Recovery and Climate Experiment (GRACE) gravity-field solutions based on the Wiener optimal filtering (Sasgen et al., 2006). Their optimal filter is designed from the least-square minimisation of the difference between the desired and filtered signals. It requires information about the power spectra of the desired gravitational signal and the contaminating noise, which is inferred from the average GRACE degree-power spectrum. Sasgen, Martinec and Fleming (2006) showed that the signal decreases with increasing spherical harmonic degree j with approximately j(-b), where b = 1.5 for GRACE data investigations. The degree power of the noise increases, in the logarithmic scale, linearly with increasing j. The Wiener optimal filter obtained for the signal model with b = 1.5 closely corresponds to a Gaussian filter with a spatial half width of 4 degrees (similar to 440 km). The authors demonstrate that the filtered GRACE gravity signal is relatively insensitive to the exponent b of the signal model, which indicates the robustness of Wiener optimal filtering.

Planetary physics

O. Cadek and co-workers published a paper in which they argued for a viscosity increase in depth in the mantle of Venus (Pauer et al., 2006). Their conclusion is based on analysis of the long- and intermediate wavelength part of the gravity and topography spectrum of the planet, carried out in the framework of internal loading theory. Together with researchers from the University in Nantes, France, O. Cadek studied effects of tidal dissipation on the total heat production in the Saturn moon of Enceladus (Tobie, Cadek and Sotin, in preparation). The preliminary results indicate, that the observed anomalous heat flow can be explained by a model with a low viscosity plume beneath the south pole of the moon and a limited region of liquid water at the interface between the icy mantle and silicate core.

Theory of seismic waves (reported by V. Cerveny and L.áKlimes)

Recent developments in seismic ray method

Important developments in seismic ray method achieved during last 20 years have been described in an extensive invited review paper by Cerveny, Klimes & Psencik (2006). The paper is devoted to the basic features of the seismic ray method, its recent extensions, and future possibilities. The topics include ray histories, two-point ray tracing, controlled initial-value ray tracing, wavefront tracing, interpolation within ray cells, paraxial ray methods, third-order and higher-order spatial derivatives of travel time, second-order and higher-order perturbation derivatives of travel time, optimization of model updates during linearized inversion of travel times, coupling ray theory for S waves, quasi-isotropic approximations of the coupling ray theory, Gaussian beams, Gaussian packets, optimization of the shape of Gaussian beams or packets, asymptotic summation of Gaussian beams and packets, linear canonical transforms, coherent state transforms, Maslov methods, decomposition of a general wave field into Gaussian packets or beams, sensitivity of waves to heterogeneities, Gaussian packet migrations, higher-order ray-theory approximations, direct computation of first-arrival travel times, ray method with complex eikonal, hybrid methods, ray chaos, Lyapunov exponents and rotation numbers, models suitable for ray tracing, application of Sobolev scalar products to smoothing models.

Seismic waves in anisotropic viscoelastic media

Considerable attention was devoted to homogeneous and inhomogeneous harmonic plane waves propagating in anisotropic viscoelastic media, and to the corresponding reflection/transmission laws.

Cerveny & Psencik (2006b) studied the properties of the energy flux of plane waves propagating in viscoelastic anisotropic media. A great attention was devoted to the energy velocity and to the loss factor. Both P and S waves were investigated, and numerical examples were presented. Cerveny & Psencik (2006c) studied the polarization of plane waves propagating in viscoelastic anisotropic media. They demonstrated that the polarization is, in general, elliptical. For homogeneous plane waves, the polarization is usually nearly linear, with large eccentricity. The eccentricity decreases with increasing inhomogeneity. Many numerical examples were presented. Time-averaged and time-dependent energy-related quantities of harmonic waves in inhomogeneous viscoelastic anisotropic media were studied by Cerveny & Psencik (2006a, submitted 2006). Several forms of the energy-balance equation were derived.

Brokesova & Cerveny (submitted 2006) studied the viscoelastic reflection/trans-mission problem for general orientation of propagation and attenuation vectors in isotropic media. Cerveny (submitted 2006) proposed the reflection/transmission laws for homogeneous and inhomogeneous plane waves at a plane interface between two homogeneous anisotropic viscoelastic media. The proposed laws are valid for unrestricted anisotropy and viscoelasticity.

Anisotropic ray theory

Great attention was devoted to dynamic ray tracing and ray-propagator matrices in layered inhomogeneous anisotropic media. Moser & Cerveny (2006, in press 2006) used the dynamic ray tracing system in Cartesian coordinates to construct surface-to-surface paraxial matrices. These 4x4 matrices define the transformation of paraxial ray quantities from one surface crossed by the reference ray, to another. They also used the surface-to-surface paraxial matrices in paraxial ray methods, and discussed certain important applications of the proposed formalism in seismology and seismic exploration. Cerveny & Moser (2006, in press 2006) studied the relations between the 6x6 ray-propagator matrices in Cartesian and ray-centered coordinates, and found useful transformations between them. They also derived the 6x6 and 4x4 interface propagator matrices in ray-centered coordinates. Cerveny (2006) compared the dynamic ray tracing systems, consisting of four linear ordinary differential equations of the first order, in ray-centered coordinates and in local wavefront orthonormal coordinates. The systems, known form literature and derived independently, consist of equations expressed in different forms. Cerveny showed that both systems are fully equivalent and may be used alternatively. Klimes (2006b) discussed six different particular choices of ray-centered coordinates in an anisotropic medium. The equations have been derived for a general homogeneous Hamiltonian of an arbitrary degree and are thus applicable both to the anisotropic-ray-theory rays and anisotropic common S-wave rays.

Klimes (2006a) derived explicit equations for the perturbations and spatial derivatives of amplitude in isotropic and anisotropic media. The perturbations and spatial derivatives of the amplitude exponent can be calculated by numerical quadratures along an unperturbed ray in the reference medium, analogously as the perturbations and spatial derivatives of travel time.

The caustic identification algorithm for isotropic media has been generalized to anisotropic media, and the rules for the phase shift of the anisotropic-ray-theory wave field due to caustics have been derived (Klimes, submitted 2006).

Coupling ray theory

Coupling ray theory is required for modelling propagation of S waves in heterogeneous weakly anisotropic media. Equations for the numerical common S-wave ray tracing and for the corresponding dynamic ray tracing in a smooth elastic anisotropic medium have been proposed (Klimes, 2006b, 2006c).

Bulant & Klimes (2006, in press 2006) 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.

Klimes & Bulant (2006a, 2006b) derived the equations for calculating the second-order perturbation expansion of travel time along the anisotropic common S-wave ray. The second-order terms in the perturbation expansion from the anisotropic common S-wave ray to the anisotropic-ray-theory rays can be used to estimate the errors due to the anisotropic-common-ray approximation of the coupling ray theory. The authors took advantage of their experience with calculating the analogous second-order perturbation expansions along the reference isotropic-ray-theory rays.

Velocity macro models and numerical ray tracing

Construction of velocity models suitable for ray tracing from field VSP measurements and from sonic velocity logs was tested. A possibility to determine the medium correlation function from sonic velocity logs was also studied.

Numerical ray tracing has been tested on various 3-D synthetic structures, including the smoothed SEG/EAGE Salt Model (Bucha, 2006c). The finite-difference seismograms in the elastic SEG/EAGE Salt Model have been compared with the ray-theory seismograms calculated using the SW3D software (Bucha, 2006b, 2006e). Bucha (2006a, 2006d) also applied the SW3D software to the construction of reflection maps on the bottom flat interface in the smoothed SEG/EAGE Salt Model with interfaces, and discussed the illumination of the large shadow area below the trunk of the salt body.

Eisner & Bulant (2006a, 2006b), Eisner, Bulant & Le Calvez (2006), and Bulant, Eisner, Psencik & Le Calvez (submitted 2006) paid attention to the borehole deviation surveys. Not performing accurate borehole deviation surveys for hydraulic fracture monitoring and neglecting the effects of the borehole trajectory results in significant errors in the calculated fracture azimuth and other parameters. For common geometries, a 5-degree deviation uncertainty of monitoring or treatment boreholes can cause more than 40-degree uncertainty in inverted fracture azimuths. The errors caused by the deviation uncertainty should not be compensated by artificially adjusting the velocity model.

Gaussian-packet prestack depth migration

Optimization of the shape of Gaussian beams (Zacek, 2006a) enables the Gaussian-packet prestack depth migration to be applied to more complex velocity models, and its accuracy to be improved.

The decomposition of the time sections into optimized Gaussian packets is of key importance in the Gaussian packet migration. The equations for the decomposition have been derived and the decomposition was numerically tested (Zacek, 2006b).

The theory used for the Gaussian packet common-shot migration was described and a numerical test was performed in the Marmousi model (Zacek, submitted 2006). The numerical test includes examples of both a single common-shot migrated section and a stacked common-shot migrated section.

CD-ROM with SW3D software, data and papers

Compact disk SW3D-CD-10 (Bucha & Bulant, 2006) 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-10 also contains over 240 complete papers from journals and from SW3D consortium research reports.

Earthquake and structural studies (reported by J.áZahradnik)

Studies of earthquakes in Greece

The program package ISOLA designed for inversion of regional waveforms into multiple point-source model of the earthquake source was upgraded. The Fortran version 06a was released on internet with numerous examples and synthetic tests ( http://geo.mff.cuni.cz/~jz/tmp/isola06a_release/). In cooperation with the Seismological Laboratory, the University of Patras, the package was also developed into a user-friendly Fortran-Matlab tool (Sokos and Zahradnik, submitted). A joint web page was created (http://seismo.geology.upatras.gr/isola/) from where the software is freely available. Numerous applications of the code in 2006 included, for example, the M6.7 Kythira in 2006, the M4.5 Amfilochia in 2002, three earthquakes M5.5 of Zakynthos in 2006, and six earthquakes of M3.4 to 4.6 in the Corinth Gulf in 2001-2005. Most of the applications were posted at http://seis30.karlov.mff.cuni.cz, presented at ESC 2006 in Geneva, and are being prepared for publication. Since 2006, ISOLA has been routinely used by the Patras University to process data from their new broad-band satellite network, and the moment tensors are reported to EMSC.

Ronnie Quintero of the OVSICORI-UNA, Heredia, Costa Rica, visited the department in December 2006 to learn the ISOLA software and to apply it in their new broad-band seismograph network.

The new software based on the stabilized EGF method and the slip patch method was developed by V. Plicka and applied to investigation of the source process of the strong Kythira earthquake (Mw=6.7), January 8, 2006, Aegean Sea. The relative source time functions of this earthquake at the regional broad-band stations were obtained and clearly indicated the source complexity. The Neighborhood Algorithm was used to invert the relative source time functions into the basic kinematic source parameters (fault plane orientation, slip value, position and size of the asperities and rupture velocity). The model provided preference of one of the nodal planes to be the fault plane (Plicka, V., submitted).

Inverse problem of the slip distribution on a fault

A new kinematic finite extent source inversion scheme, introducing an innovative parametrization of the problem, was developed by J. Burjanek and presented at ESC 2006, Geneva. The spatial slip distribution is represented by overlapping 2D Gaussian functions on regular grid. The linear inversion for static slip is done with positivity constraint, so called "Quadratic programming" is applied to solve such problem. The method benefits from favourable spectral properties of Gaussian function. The latter prevent from employment of artificial smoothing operators. Thus a direct insight into spectral properties of the static slip distribution of real earthquakes is obtained.

The method was successfully validated in the framework of the SPICE blind test (organized by M. Mai, ETH Zurich) on the slip retrieval on a fault, an important experiment to study uncertainties of this common seismological inverse problem. The other submitted solution of our department was that of J. Zahradnik, who combined inversion by ISOLA software (above) with the forward simulation of a finite extent source. Both submitted reports are available on the web.: http://geo.mff.cuni.cz/~burjanek/report_blind_test_jb.pdf http://geo.mff.cuni.cz/~jz/tmp/presentations/kg2006/jz_report.doc

Modeling the dynamic stress field

An important first step to enrich the scope of the department by studies of the source dynamics has been made by J. Burjanek, the PhD student. It is based on the idea that the kinematic source models reproduce source complexity observed in the earthquake slip inversions. Therefore, the analysis of the stress field associated with the slip history prescribed in the advanced kinematic models is studied. J. Burjanek performed a detail numerical study of of the full dynamic stress history on a fault by the boundary element method, formulated in spectral domain. Basic characteristics of the resulting stress field were discussed in relation to recent dynamic source models of real earthquakes. The study included the static stress drop, the strength excess and the dynamic stress decrease distributions corresponding to the kinematic k-squared stochastic source model with k-dependent rise time. A new quantity, the stress delay, was introduced. The paper passed the first stage of the review process in the Journal of Geophysical Research, and is expected to appear in 2007.

Modeling strong ground motions

PhD thesis on this exciting topic was defended in 2006 by F. Gallovic. In Gallovic and Burjanek (2005; accepted forápublication in Annals of Geophysics), two common approaches to strong ground motion simulations due to finite-extent sources are discussed: the integral and composite approach. It is shown that the integral model overestimates empirical scatter of peak ground acceleration, which is explained by overestimated directivity effect. In Gallovic and Brokesova (2006; accepted forápublication ináPEPI), a new hybrid approach is proposed and applied to modeling of the 1997 M6.1 Kagoshima earthquake, Japan. In cooperation with Italian colleagues the model has been also applied to the 2002 M5.8 Molise earthquake (Franceschina et al., 2006).

The hybrid modeling has been further applied to the probabilistic aftershock hazard assessment (PAHA), see Gallovic and Brokesova, submitted to J. Seismology. PAHA intrinsically combines: 1) probability of the aftershock occurrence (usually described by the Poissonian process with mean rate governed by the generalized Omori's law), and 2) statistical description of strong ground motions (in terms of occurrence of a given strong motion characteristic at a given receiver depending on earthquake location and magnitude). The application of PAHA, modified for finite-extent sources, was illustrated on an example of the A25 aftershock (M5.8) of the Izmit earthquake. In this case the empirical attenuation relations have been substituted by synthetic ones, which in principle allows taking into account structural models (e.g. sedimentary basins) and the activated fault with likely aftershock and its mechanism. This in turn needs utilization of strong motion model tuned for given area.

Modeling 3D local site effects (finite-difference method)

Local 3D site effects due to surface topography and complex geometry of the bedrock interface were investigated by hybrid method developed in recent years by I. Oprsal. His software, able to manage realistic source-path-site configurations and wavefields of the engineering interest (up to 10 Hz), has attracted a number of scientists and extended the collaboration of the department quite significantly.

The Roman city of Augusta Raurica, Switzerland probably suffered from a strong earthquake in the year 250 AD. In order to explain the damage, modeling of site effects at alluvial sediment sites with pronounced topography was performed for this site. Resulting strong 3D effects and spectral amplifications reach a factor 9 in the frequency band 0-10Hz; correlation of high spectral amplification spots and archaeological findings were indicated (Oprsal and Faeh, 2006, Faeh et al., 2006).

Together with S.M. Richwalski and R.Wang (both GFZ Potsdam) we combined the finite-extent source modeling and local site effects. The target zone of this 3D study was the Tricarico town area (Southern Italy). The town is situated on a calcarenite outcrop, which is underlain by argillaceous material with lower shear wave velocity. Analysis of the microtremor data showed several features related to site effects, too. A detailed 3D model of the sub-surface underneath the town was constructed, including topography. The respective earthquake hazard is moderate to high according to the Italian official seismic classification (0.25 g for a return period of 475 years). Damage was reported during the 1694, 1857, and 1980 earthquakes. The results were presented at EGU Vienna and ESC, Geneva. A paper is under preparation.

Application of the 3D code to the southern part of the city of Rome (S. Paolo Cathedral) is still under tests because it comprises parallelization (MPI) of the code. So far the 35-processor computation shows a promissing speedup of 25 (cooperation of A. Caserta, INGV, Rome). This team-demanding work will continue and the results are due to be published within 2007.

A simple source modeling based on the published slip distribution was applied in the 2005-2006 blind prediction experiment. The experiment focused on the M6 (Sep 28, 2004) Parkfield earthquake, with special emphasis to site effects at the Turkey Flat test site. Our contribution was the only one among those applied by the participants which focused on the finite-extent source effect. The result were presented and published in Oprsal and Zahradnik (2006).

Structural studies in the Gulf of Corinth, Greece

Based on the Corinth Rift Laboratory local seismic network, CRL, quality of the routine seismic location can be increased. Station corrections for that purpose were derived and published by Jansky et al. (2006). The abundant data of the CRL network allow also non-standard assessment of the velocity structure based on non-linear inversion of the P and S-wave travel times. Application of the Neighbourhood Algorithm for this purpose was studied in synthetic tests (paper by Jansky et al. Submitted to Acta Geodynamic et Geomaterialia) and the practical results were accepted for publication (Jansky et al., J. of Seismology).

Structural studies in the Czech Republic

In the paper by K. Holub et al. (Acta Geodynamica et Geomaterialia), co-authored by O. Novotny of our department , the structure of the uppermost part of the Earth┤s crust in the territory of northern Moravia and Silesia was investigated. Quarry blasts and mining induced seismic events served as seismic sources. Permanent, temporary and portable seismic stations were used for the monitoring of these seismic events. Data of the CELEBRATION 2000 and SUDETES 2003 refraction experiments were incorporated, as well. The velocity-depth dependence of body waves was searched by joint inversions of travel times of Pg/Sg phases. In this way, velocity structure was determined down to a depth of 5-6 kilometers.

Seismic refraction measurements along several profiles in the western part of the Ohre (Eger) Rift and in the Moravo-Silesian region were interpreted by Holub et al. (paper in transations VSB Ostrava). The travel-time curves indicate a significant azimuthal dependence, the wave propagation along the rift being faster than that in the perpendicular direction. Remarkable lateral inhomogeneities are also documented. The mean travel times were inverted in terms of vertically inhomogeneous models (1-D models). For this purpose, O. Novotny suggested to approximate the travel-time curves by rational functions (quotients of two polynomials), the parameters of which were computed by the conjugate-gradient method. The classical Wiechert-Herglotz method was then used to interpret the smoothed travel-time curves. The velocity distributions were obtained down to a depth of about 2ákm.

P. Kolinsky, the PhD student supervised by J. Brokesova, successfully applied his sophisticated new software designed for the frequency analysis of seismic waves. The target zone of their Love wave dispersion investigation was the western Bohemia uppermost crust and its shear wave velocity structure. Their paper was submitted to the Journal of Seismology.

Publications: http://geo.mff.cuni.cz/documents/publ_06.htm