Research overview 2009


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 2009, there were 20 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 2009, the Department participated in two EC projects: C2C and IMAGES. C2C (2007-2011), 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. Preparations were also carried to join one new European project, QUEST, Quantitative Estimation of Earth’s Seismic Sources and Structure (coordinated by H. Igel).


Research project SW3D, Seismic Waves in Complex 3-D Structures, coordinated at the Department of Geophysics since 1993, has continued successfully in 2009. The project has been supported by 5 companies or research institutes (BP America Inc., U.S.A.; Chevron U.S.A. Inc., U.S.A.; NORSAR, Norway; Petrobras, Brazil; Schlumberger Cambridge Research Limited, U.K.) in the framework of the international SW3D consortium.


Seismic station PRA (created in Prague in 1924), equipped with CMG-3ESP, was linked with the international data center ORFEUS to provide on-line data transfer to Utrecht. At the end of 2009, a new broadband seismograph CMG-3T substituted the former instrument. One more seismic station in Greece, Prodromos PDO, joined the family of the stations operated by the department in cooperation with theUniversity of Patras. Each of the stations SER, MAM, LTK, PYL and PDO is equipped with CMG-3T (broadband) and CMG-5T (strong motion) instruments. PYL, LTK and PDO, belonging to the satellite network PSLNET, transmit the signal continuously to the unified Greek network HUSN; LTK data are also online provided to ORFEUS. The department also participated in installation and operation of the borehole (60 m deep) broadband CMG-3T station at Geodetic observatory Pecny, GOPC.


The theses defended at the Department in 2009 included the following:

B.C.: Helena Munzarová, Libor Šachl

M.S.: Ľubica Valentová


The department was happy to host the long-term visitors Olivia Golle, Alexandra Guy, Nicola Tosi within the framework of the C2C project and Ernestina Lucca.


As in previous years, research at the Department was carried out in three directions: Theory of Seismic Waves, Earthquake and Structural Studies, and Geodynamics.



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

L. Klimeš)


Paraxial ray methods and Gaussian beams in anisotropic media


The paper by Červený and Moser (2009) is 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.


Klimeš (2009a) derived the relations between the propagator matrix of geodesic deviation (i.e., paraxial-ray propagator matrix) and the second-order derivatives of characteristic function (i.e., second-order derivatives of two-point travel time) in general coordinates. The equations are applicable to Finsler geometry, to Riemann geometry, and to their various applications like the general relativity or the high-frequency approximations of wave propagation.


The paper by Červený and Klimeš (2009) is devoted to the second derivatives of travel-time fields in anisotropic media, computed by dynamic ray tracing (DRT) along a reference ray. Two forms of DRT system have been broadly used in applications: the DRT system in global Cartesian coordinates and the DRT system in ray-centred coordinates. Simple transformation relations between the second derivatives of the travel-time field in global Cartesian coordinates and in ray-centred coordinates are derived. These transformation relations can be used in many applications, including the computation of complex-valued paraxial travel times which are necessary in the evaluation of Gaussian beams.


The paper by Červený and Pšenčík (2009a) is devoted to the theory of Gaussian beams concentrated close to rays of high-frequency seismic body waves propagating in inhomogeneous anisotropic layered structures. The basic role in the computation of Gaussian beams is played by a 2x2 complex-valued matrix of second derivatives of the travel time field along the reference ray. This matrix can be computed in several ways by dynamic ray tracing. Different possibilities are discussed in detail and possible simplifications are outlined.


Seismic waves in anisotropic attenuating media


Considerable attention has been devoted to perturbation methods from perfectly elastic media to anisotropic attenuating media. In the paper by Červený and Pšenčík (2009b), 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.


Calculation of complex-valued travel time


In real space, the eikonal equation for complex-valued travel time represents the system of two Hamilton-Jacobi equations for the real and imaginary parts of the complex-valued travel time. Klimeš (2009b) suggested to solve these equations numerically using the wavefronts obtained by wavefront tracing.


The numerical solution of the eikonal equation for complex-valued travel time led Klimeš (2009c) to the idea of incorporating the amplitude into the complex-valued travel time by moving the second-order travel-time derivatives from the transport equation to the eikonal equation.


Gaussian packet prestack depth migration


Bucha (2009) continued the work on the Gaussian packet prestack depth migration, and tested the Gaussian packet prestack depth migration in two simple models composed of homogeneous blocks separated by planar or curved interfaces.


Structural studies


Bulant (2009) smoothed the existing 1-D regional velocity model of the Bohemian Massif for ray tracing. He minimized the Sobolev norm composed of the second-order velocity derivatives in order to prevent ray chaos characterized by the corresponding Lyapunov exponents.


The medium correlation function is of principal importance in refraction travel-time tomographic inversion. Travel times measured during refraction seismic experiments contain information on the correlation function primarily in the horizontal directions, whereas sonic-log measurements and vertical seismic profiling in deep wells contain information on the correlation function in the vertical direction. Bulant and Klimeš (2009) inform about the available data from sonic-log measurements in the very deep KTB well, which could be used for estimating the correlation function for the western Bohemia region in the vertical direction.


CD-ROM with SW3D software, data and papers


Compact disk SW3D-CD-13 (Bucha and Bulant, 2009) 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-13 also contains over 350 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 were most intensively investigated in 2009: the damaging M6.3 earthquake of June 8, 2008 in Western Greece, called Andravida or Movri Mountain earthquake, the disastrous L’Aquila event M6.2 of April 6, 2009, Italy, and two M>3.5 earthquakes belonging to the 2008 swarm in West Bohemia (Czech Republic).


Source studies


A new strategy to reveal the spatio-temporal evolution of the earthquake rupture process from near-regional data, without assuming a constant rupture velocity was proposed by Gallovič et al. (2009). The approach is based on the conjugate gradient method, for which the required waveform-misfit derivative with respect to slip on the fault is analytically expressed. The derivative is given by the back-propagation of residual seismograms towards the source; hence similarity with the adjoint tomographic inversions. A good initial source approximation is necessary, being obtained from hypocenter location and centroid-moment tensor solution (see ISOLA software discussed below). The iterative approach then gradually reveals major characteristics of the source process. The method was applied to investigate a line-source model of the damaging Mw6.3 Movri Mountain, June 8, 2008, earthquake in Greece. It revealed the predominantly unilateral rupture propagation toward north-east direction, with two or three main slip patches. One of them was significantly delayed, thus indicating a temporary rupture arrest. The region of the largest slip coincided with the region of least abundant aftershocks between hypocenter and centroid. The method is very fast, so may find its application for shake-maps, emergency response, and/or aftershock hazard assessment.


The so-called H-C geometrical method for quick identification of the fault plane, developed at the department in 2008, has been applied to a M6.3 event south of Crete (July 1, 2009). The same technique was also used for the L’Aquila event M6.2 of April 6, 2009, Italy. Quick reports were sent to EMSC, see Earthquake News and Highlights on the web page The L’Aquila event was processed (upon invitation) in co-operation with the University of Naples, together with E. Sokos of the University of Patras. However, due to normal mechanism of the event, the H-C method experienced many problems. In fact, the idea was finally used in its ‚reverted‘ form: when the fault plane was unequivocally determined by the satellite interferometry, the assumption about joint plane comprising both the hypocenter and centroid was used to precise the centroid position.


Software ISOLA for moment-tensor calculations from regional waveforms has been further upgraded ( For example, an option of multiple crustal models was launched. The method is routinely applied by E. Sokos, University of Patras, for earthquakes well recorded in Western Greece by satellite network PSLNET, most M>3, and the results are reported to EMSC. Consultations were provided to many external users of the software, too.


The weak earthquakes of the West Bohemia swarm in 2008 were studied from the viewpoint of the near-field source terms. Data of two collocated instruments CMG-40T and STS-2 of the Geophysical Institute, Academy of Sciences, were used; the station Nový Kostel (NKC), belonging to the Czech Regional Seismological Network, was situated at a very small epicentral distance, of about 3 km. The records revealed special long-period disturbances. After removing these artifacts, the records enabled assessment of the very small static displacements of the order of 10-5m. The paper (Zahradník and Plešinger) was accepted in BSSA and will appear in 2010.


The study of M5 Amfilochia earthquake, 2002, performed in previous years, was finally published by Adamová et al.(2009). Large non-shear components reported by agencies were explained as an artifact due to complexity of the event.


A new method to measure the size of the non-elastic source volume, based on spectral magnitude of the P- and S-waves, was suggested (Duda and Janský, 2009). The technique was successfully applied also to the catastrophic earthquake of 12 May 2008, Sichuan, China.


Rotational motions


Modern topic of the rotational components of seismic motions was dealt by Brokešová et al. (2009); two instruments and a methodology were developed and patented: the rotational sensor, the generator of rotational waves, and technique to apply these instruments in the seismic exploration. The instruments have been under further intensive testing, including, for example, the swarm region of the Western Bohemia.


Wang et al. (2009) investigated combined source and basin effects upon the rotational seismic ground motions. The latter have a great practical significance for earthquake effects upon highway bridges and other long structures.




Determination of the depth of very shallow events, of the order of 5 km or less, is a particularly difficult problem. It is often hampered by missing details of the structural models (poor knowledge of the low-velocity surface layers), and missing stations at epicentral distances comparable to the source depth. Erroneous locations (i.e. apparently larger depths) might significantly bias seismo-tectonic deductions about so-called seismogenic layers. Also the tsunami computationally-based scenarios might be biased due to the same effect. The non-reviewed discussion by Janský et al. (2009a) issues an important warning of this kind. The paper is based both on observations and synthetic tests. Usefulness of the combined location and moment-tensor studies is an important finding of this paper. Success of this technique follows from the fact that it makes use of both the high- and low-frequency seismic information. Stable depth estimate from the relatively low-frequency moment tensors, independent of the inaccurately known subsurface layers, has not been broadly recognized yet as a tool to improve the accuracy of the shallow seismicity pattern.


Similar problems, complicated also by instrumentation and site effects, may accompany also the very active region of the western Gulf of Corinth; the efficiency of the stations locating events in that region was quantified by Janský et al. (2009b).


Strong ground motions


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., 2009) was investigated by this technique, comparing results with independent method of F. Pacor.


Zollo et al. (2009) used synthetic ground motion scenarios calculated by F. Gallovič to test their Early Warning system in the Irpinia region of Southern Italy. Two new parameters (so-called Effective Lead-Time and Probability of Prediction Error) to measure the system performance were introduced and their regional maps drawn. The latter show a significant variability around the fault up to large distances, thus indicating that the early warning system’s capability to accurately predict the observed peak ground motion strongly depends on distance and azimuth from the fault.


A 3D hybrid method to simultaneously address the seismic-wave propagation and site effects, upgraded with a new scheme, was better described from the point of view of its mathematical formulation (Opršal et al., 2009). The main innovation of the paper is the generalization of the boundary condition acting between the two complementary sub-regions, describing the regional and local wave propagation, respectively. The 3D code has been applied to several sites in Italy as well as to data of an international test site. Parallelization of the code, to enhance its application range, was performed and published by Caserta et al. (2009).


Structural studies


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., 2009). 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-Herglotz method. The model is useful for location of induced seismic events in the area of the Příbram-Háje underground gas storage, and also for weak earthquakes occurring occasionally in the Orlík reservoir region on Vltava river.



Geodynamics (reported by C. Matyska)


Mantle dynamics


In the two papers, Tosi et al. (2009a,b) addressed the role of lateral viscosity variations (LVV) in predicting the long-wavelength geoid. The effect of LVV on the geoid has been a puzzling problem since early 90s. Until recently, no really conclusive answer was available, mostly because the numerical codes were not able to treat LVV with a sufficiently high resolution and the prediction of the geoid was only limited to cases with rather small lateral viscosity contrasts. Using a high-resolution spectral code developed by Tosi and Martinec (2007), Tosi et al. (2009) investigated the sensitivity of the geoid to highly viscous slabs in the mantle in 2D axisymmetric spherical geometry. They found strikingly large impacts of LVV on the geoid and demonstrated a paradoxical effect, namely that the amplitudes of the geoid above subduction zones would be unrealistically large (100 m and more) if the viscosity of slabs in the lower mantle were significantly more viscous than that of the surrounding mantle. They proposed several possible scenarios of behavior of slabs in the lower mantle and discussed them in connection with thermal convection models (e.g. Běhounková and Čížková, 2008).


The recently discovered iron spin transition in major mantle minerals at high pressures should exert dramatic influences on the transport properties, such as the activation creep parameters in the deep mantle. As a consequence of the spin transition, there is a significant softening of the elastic parameters and an attendant strong reduction of the activation energy in the creep law, leading to its non-monotonic behavior to develop in the middle of the mantle. Matyska et al. (submitted) studied the dynamical consequences of the activation energy reduction within the framework of a 2-D Cartesian convection model. They used a parameterized model of the rheological parameters capturing the basic physics and found that the non-monotonicity in the rheological parameters led to the formation of viscosity minimum at a depth of about 1600 km

followed by a ``viscosity hill'' in the bottom half of the lower mantle. The numerical simulations revealed a tendency to the formation of superplumes and/or plume clustering in the deep lower mantle with the dynamical behavior, which was influenced by the changes of the activation parameters throughout the lower mantle and other material properties at the bottom of the mantle (decrease of thermal expansivity, viscosity reduction due to the perovskite-postperovskite phase change, radiative thermal conductivity).


The results of the modelling study about the dynamical consequenses of the low viscosity postperovskite in the lowermost mantle were accepted for publication (Čížková et al., in press).


Planetary research


Important though indirect information about the internal structure of Venus is provided by its topography and geoid. This information has been used in last decades to constrain the Venus mantle viscosity structure and its dynamic regime. Recently, Pauer et al. (2006) performed the geodynamic inversion of the venusian geoid and topography and proposed a group of best fitting viscosity profiles. As there is no information available about the Venus mantle density structure, they assumed a simplified 2D depth-independent approximation of the density anomalies. Benešová and Čížková (to be submitted) use their viscosity model together with several other simple viscosity stratifications (linear increase, constant, constant with a stiff lithosphere) as an input for the mantle convection code. They carry out the simulations of the mantle evolution in both 3D spherical shell and a 2D axisymmetric model and check, whether the character of the dynamic topography and the geoid represented by their power spectra fits the observed quantities. They estimate the effects of the thermally-dependent viscosity, internal heating, adiabatic and viscous heating and varying Rayleigh number. Though the geoid and topography do not vary too much among their considered viscosity profiles, they generally conclude, that the best fit between the observed and model quantities is reached for the viscosity profile by Pauer et al. (2006) with the upper mantle viscosity of 1021 Pa s. Thermally dependent viscosity improves the fit. The observed geoid and topography over several selected plumes successfully fit the model results reached for the preferred viscosity profile.


Zábranová et al. (2009) have studied equations and boundary conditions describing tidal deformations and the gravitational potential of prestressed elastic spherically symmetric bodies that were decomposed into a series of boundary value problems for ordinary differential equations of the second order and that were converted into a set of linear algebraic equations with an almost block diagonal matrix by using pseudospectral schemes on Chebyshev grid. Both accuracy and efficiency of the method were shown by evaluating the tidal Love numbers for the Earth's model PREM. The method was applied to evaluation of the tidal Love numbers for models of Mars and Venus. The Love numbers of the two simplified Martian models - the former optimized to cosmochemical data and the latter to the moment of inertia - were h2=0.172 (0.212) and k2=0.093 (0.113). For Venus, the value of k2=0.295, well-known from the gravity-field analysis, was consistent with the results for the model with the liquid-core radius of 3110 km (Zábranová et al., 2009).


In cooperation with researchers from the University of Nantes, France, M. Běhounková developed a new important tool to study convection in planets and icy satellites driven by the heat released due to tidal dissipation. She successfully combined the 3D spherical code of Oedipus, developed by G. Choblet, for simulating thermal convection, and the code of O. Čadek, designed for evaluating tidal dissipation in planetary and moon mantles with a general 3D viscosity structure. The new code, constructed and successfully tested by Běhounková et al. (2009), is unique since, for the first time, it allows the effect of the tides on mantle convection to be incorporated in a self-consistent manner. The code has a large number of applications in planetary research. It is now used to study thermal evolution of Earth-like exo-planets and it is to be soon applied to investigating processes in some icy satellites (Europa, Enceladus) as well.


EM induction modelling


Using the technique of adjoint modelling, we have inverted one year of CHAMP satellite data in terms of 1-D and 2-D axisymmetric conductivity models of the upper and mid-mantle (Martinec and Velímský, 2009). The results suggest slightly larger average conductivity of the northern hemisphere in the upper mantle.


Seven years of CHAMP satellite data were also inverted in terms of 1-D conductivity with emphasis on the deepest part of the Earth's mantle, the D'' layer (Velímský, in review). No significant conductivity increase was observed. This can be explained (as confirmed by 3-D modelling) as a lack of interconnection of the highly conductive post-perovskite phase in the longitudinal direction, i.e., in the dominant direction of polarization of the external currents.


The oceans play a special role in electromagnetic induction due to thein relatively high conductivity and the dynamo effect of ocean currents. The magnetic field by ocean circulation motion can be divided into toroidal and poloidal parts. The toroidal magnetic field is generated by electric currents closing in vertical planes and is estimated to reach 100 nT in amplitude. The much weaker poloidal field, with amplitudes up to 10 nT, results from electric currents closing horizontally. It has a significant vertical component and reaches remote land and satellite locations. Much attention has been given to the periodic magnetic signals of ocean flow which is driven by the lunar tides, but no attention has been devoted to the induced toroidal field. Dostal and Martinec (submitted) developed the matrix propagator technique to compute the toroidal magnetic field inside the Earth and oceans with the aim to generate the secondary poloidal field due to electrical conductivity inhomogeneities in the Earth‘s crust and lithosphere. The first estimate of its strength is up to 10 nT that may be detectable, summed up with the primary poloidal field by future SWARM satellite mission.


The 3-D time-domain method designed for inversion of Swarm multisatellite data was further developed. An extension study comissioned by the European Space Agency, "Level 2 products and performances for mantle studies with Swarm", was successfuly completed. We have joined the consortium of leading European institutions in the proposal to ESA to run the Swarm Level 2 processing facility.


Glacial isostatic adjustment


The improvement in understanding of dynamic processes in the Earth‘s mantle demands to consider a non-linear rheology of mantle material. Whereas this rheology is accepted in the studies of mantle convection, the need of a non-linear material behaviour in modelling of glacial isostatic adjustment (GIA) is still under discussion. Almost all the predictions of ongoing present-day processes induced by GIA are based on the assumption of a linear Maxwell viscoelastic rheology. To study the influence of non-linear rheology on the GIA-induced motion, Martinec and Klemann (submitted) implemented a non-linearly stress-dependent rheology in the spectral finite-element formulation of a viscoelastic self-gravitating sphere. The main effect of a non-linear rheology on the GIA-induced motion is for times when a surface ice-mass load is changing most rapidly because of large induced loading stresses.


Sasgen and Martinec (submitted) performed a joint inversion of gravity fields from the Gravity Recovery and Climate Experiment (GRACE) for glacial-isostatic adjustment over North America and present-day ice-mass change in Alaska and Greenland. We determine the Alaskan and Greenlandic contribution to sea-level change from the adjusted ice-mass change models. The residual misfit over the GIA-dominated region around the Hudson Bay is interpreted with regard to the mantle viscosities beneath North America by applying forward model calculations of the GIA signal in these regions. We compare our results based on satellite gravimetry with constraints derived from sea-level indicators, absolute gravimetry, tide-gauge stations and GPS, and show their sensitivity to the GRACE release considered, as well as to the glacial history underlying the GIA forward model.


The Gravity Recovery and Climate Experiment satellite mission has allowed the resolution of temporal variations in the Earth‘s gravity field to serve as a new observable for monitoring mass changes in cryosphere. Sasgen and Martinec (submitted) analysed the GRACE time series from August 2002 to August 2008 with regards to regional ice-mass variability in Greenland. We found that the mass change of the Greenland Ice Sheet amounts to 0.5 mm/a equivalent sea-level change and significantly accelerated during the observation period. The comparison with Interferometric Synthetic Aperture Radar (InSAR) data and output from surface-mass balance modelling indicates that mass-loss acceleration is mainly caused by increasing discharge in the Northwest starting in the year 2005. In the year 2007, mass loss additionally accelerates in the Southwest caused by a reduced surface-mass balance. We conclude that GRACE allows for the detection of regional scale mass variations, including accelerations, and may further contribute to the understanding the processes governing the current changes of Greenland Ice Sheet.


Modern modelling approaches to GIA are based on several techniques ranging from purely analytical formulations to fully numerical methods. Various European teams nowadays are independently working on the post-glacial rebound process in order to constrain the rheological profile of the mantle and the extent and chronology of the late-Pleistocene ice sheets which are prerequisites for the determination of the GIA contribution to geodetic observables. Martinec contributed to the benchmark study performed within the Working Group 4 of the ESF COST Action ES0701. Improved constraints on models of Glacial Isostatic Adjustment focuses on i) load Love numbers and relaxation spectra, ii) the deformation and gravity variations driven by surface loads characterized by simple geometry and time-history, and iii) the rotational fluctuations in response to glacial unloading.


Ice-Sheet dynamics


Souček and Martinec have incorporated their newly-proposed numerical algorithm SIA-I (Souček and Martinec, 2008) into a large-scale thermo-mechanical ice-sheet model. The mechanical part of the problem, that is the momentum equation and rheology are solved by the SIA-I technique, which improves the standard Shallow-Ice Approximation in a series of iterative steps. In addition, the numerical model involves the evolution equation for the free surface and the heat-transport equation. Preliminary testing of the model was completed by performing two series of European Ice-Sheet Model INiTiative (EISMINT) benchmarks focused on (i) evaluating the effects of thermo-mechanical coupling and (ii) investigating the dynamics of the Greenland Ice-Sheet (GIS) subject to various climatologic scenarios. A simulation coupling the dynamic ice-sheet model with the GIA model by Martinec et al. (2000) is in preparation. Also a simulation evaluating the effects of higher-order (non-shallow) dynamics to surge occurrence and characteristics for the ISMIP-HEINO experiment is being set up.


Bohemian Massif dynamics


Until recently, the geodynamical research at the Department was only focused on studying global and large-scale problems. Fast development of numerical tools and computer power now allows to model also regional and local evolution of geological systems and, in this way, to test geological hypotheses and theories. The challenging aspect of this research lies in a great amount of data collected by the geological community in the past that can be used to constrain the numerical models. P. Maierová, a PhD student at the Department of Geophysics, modified the finite-element software Elmer in order to address various geological problems related to the evolution of the Bohemian Massif during the Variscan orogeny. In cooperation with geologists from other institutions (S. Ulrich, O. Lexa, K. Schulmann), she developed a model of the eastern part of the Bohemian Massif (collision of the Moldanubian domain with Brunia) and she also studied hypothetical lower crust diapirism in the southern part of the Bohemian Massif. The publication of the results is being prepared.