RESEARCH INTERESTS

Alik Ismail-Zadeh


What do I do in science?

A principal objective of my research is to understand interplay between crustal and mantle dynamics and surface processes by means of theoretical and experimental studies. My approach is three-pronged. 1: I study geophysical phenomena on the basis of available observations. 2: I develop mathematical and numerical methodologies to model geodynamic processes (such as deformations and stresses in and around the descending lithosphere) or fluid flow (such as lava flow); thermal convection in the mantle; sedimentary basin evolution and salt diapirism; and gravitational, thermal and buckling instabilities of geological structures. 3: I test model predictions against observations and explain specific features of the processes and phenomena.

Why should the general public be interested in what I do?

When nature unleashes itself (resulting in earthquakes, volcanic eruptions, lava flows, landslides or tsunami), it strikes where and when it decides. Extreme events strike without warning and cause devastation in loss of human life and environmental damage. There is no doubt that man will never be able to prevent these occurrences entirely. However, geoscientists can gain a better understanding of the complex mechanisms behind the geohazards, dynamics of the lithosphere and its surface manifestations. Convolving the knowledge of physical processes with the knowledge of vulnerability and exposure we can construct a map of risks alerting policymakers and society on disaster risk reduction.

Why does it interest me?

As a geophysicist, I interest it to understand a complex dynamic behavior of the Earth. As an applied mathematician, I interest it to bring formalism to Earth sciences. As a habitant of our planet, I interest it to know how to cope with disasters.

 

CURRENT to PAST RESEARCH

Sedimentary basins long heated by a fading mantle upwelling

Using mantle plume modelling and global plate motion reconstructions, we have demonstrated that the East European lithosphere is likely to have been located above the weakening (due to a diffusive decay) Perm Anomaly upwelling for approximately 150-200 million years ago, as indicated in the figure (below). As the East European platform moved over the Perm Anomaly during the post-Jurassic period, the vertical tectonic movements recorded in sedimentary hydrocarbon-rich basins indicate either uplift (no sedimentation) or insignificant subsidence. The analysis of heat conduction across the lithosphere reveals that the basins have experienced prolonged heating due to the Perm Anomaly upwelling, resulting in favourable circumstances for the maturation of hydrocarbons. The research findings indicate a robust correlation between a mantle plume that heats the thick continental lithosphere over an extended geological period, the sedimentary basins heated by this plume becoming thermally weak with time, and the subsequent development of hydrocarbon provinces within these basins. The generation of oil and gas fields can occur on oceanic passive margins along the path of mantle upwellings. Hydrocarbon resources discovered in the North Sea and the Norwegian margins may be associated with the Eifel hotspot activity during the Middle-Late Jurassic times. Our model has the potential to provide an explanation for the recent findings of significant resources offshore Brazil and Africa along the track of the Tristan da Cunha plume. This plume led to the Parana-Etendeka large igneous province development approximately 132 million years ago, and subsequently triggered the opening of the South Atlantic Ocean. Our research indicates that shallow or deep hydrocarbon fields (and potentially hydrogen fields) can be concentrated in sedimentary basins where deep-source heat flow enhancements occur. This is primarily owing to the protracted presence of mantle upwelling underneath these basins. Published in Nature Communications, 15, 3915, 2024.

A schematic evolution of the mantle plume (panels a-c) associated with the Perm Anomaly (panel d) since the Early Jurassic times. Tracks of the projection of the Perm Anomaly plume’s centre on the surface of western Eurasia from the present position of this centre (51°N, 52°E) to its position in the Early Jurassic (~200 Ma); the tracks are derived from three models (marked by brown crosses, blue circles, and red diamonds) of global plate models (panel e). Vertical tectonic movements in East European sedimentary basins (panel f) have been derived from the analysis of sediments in many drilled boreholes (yellow dots in panel e).

Natural hazards, climate change and drivers of disasters

Many nations face challenges in assessing, understanding, and responding to the time-dependent nature of disaster risk. Changes in the intensity of occurrences of extreme events coupled with changes in vulnerability and exposure alter the impacts of natural hazards on society in mostly negative ways. An interrelationship between natural hazard (NH), climate change (CC), vulnerability (V), exposure (E), and decisionmaking (DM) is considered. While NHs trigger disasters and CC is likely to intensify occurrences of disasters, V and E present major drivers of disasters. Informed DM on disaster risk reduction should be based on scientific evidence from NH and CC, knowledge of V and E, and relevant options for actions on preventive disaster measures as a part of preparedness and public awareness. Published in Natural Hazards, 111, 2147–2154, 2022.

Lava dome morphology and viscosity inferred from data-driven numerical modeling of dome growth at Volcan de Colima, Mexico during 2007-2009

Magma extrusion, lava dome growth, collapse of domes, and associated pyroclastic flow hazards are among important volcanological studies. We analyze the influence of the magma viscosity and discharge rates on the lava dome morphology at Volcan de Colima in Mexico during a long dome-building episode lasting from early 2007 to fall 2009 without explosive dome destruction. Camera images of the lava dome growth together with recorded volumes of the erupted lava have been used to constrain numerical modeling and hence to match the history of the dome growth by nudging model forecasts to observations. Our viscosity model incorporates crystal growth kinetics and depends on the characteristic time of crystal content growth (or CCGT) and the crystal-free magma viscosity. Initially, we analyze how this viscosity, CCGT, and the rate of lava extrusion influence the morphology of the growing dome. Several model scenarios of lava dome growth are then considered depending on the crater geometry, the conduit location, the effective viscosity of dome carapace, and the extrusion rates. These rates are determined either empirically by optimizing the fit between the morphological shape of modeled domes and that of the observed dome or from the recorded lava dome volumes. The maximum height of the modeled lava dome and its horizontal extent are in a good agreement with observations in the case of the empirically-derived extrusion rates. It is shown that the topography of the crater at Volcan de Colima is likely to be inclined toward the west. The viscosity of the modeled lava dome (~10^12 Pa s) is in a good agreement with the effective viscosity estimated experimentally from lavas of Volcan de Colima. Due to the interplay between the lava extrusion and the gravity forces, the dome reaches a height threshold, and after that a horizontal gravity spreading starts to play an essential role in the lava dome evolution. The model forecasts that the dome carapace of higher viscosity (~10^14 Pa s) influences the dome growth and its morphology during long dome-building episodes by retarding horizontal advancement and developing steep-sided eastern edge of the dome at the volcano. The developed model can be used in assessments of future effusive eruptions and lava dome growth at Volcan de Colima or elsewhere. History matching modeling of lava dome growth sheds a light on dynamic processes inside the dome and may assist in assessing stress state in the dome carapace and in forecasting the dome failures. The results are published in Frontiers in Earth Science, 9, 735914, 2021.

Inferring lava viscosity from lava dome morphology

The viscosity of magma nonlinearly depends on the volume fraction of crystals and temperature. We estimate magma viscosity based on a comparison of observed and simulated morphological forms of lava domes. We consider a two-dimensional axisymmetric model of magma extrusion on the surface and lava dome evolution, and assume that the lava viscosity depends only on the volume fraction of crystals. The crystallization is associated with a growth of the liquidus temperature due to the volatile loss from the magma, and it is determined by the characteristic time of crystal content growth (CCGT) and the discharge rate. Lava domes are modeled using a finite-volume method implemented in Ansys Fluent software for various CCGTs and volcanic vent sizes. For a selected eruption duration a set of morphological shapes of domes (shapes of the interface between lava dome and air) is obtained. Lava dome shapes modeled this way are compared with the observed shape of the lava dome (synthesized in the study by a random modification of one of the calculated shapes). To estimate magma viscosity, the deviation between the observed dome shape and the simulated dome shapes is assessed by three functionals: the symmetric difference, the peak signal-to-noise ratio, and the structural similarity index measure. These functionals are often used in the computer vision and in image processing (AI methods). Although each functional allows to determine the best fit between the modeled and observed shapes of lava dome, the functional based on the structural similarity index measure performs it better. The viscosity of the observed dome can be then approximated by the viscosity of the modeled dome, which shape fits best the shape of the observed dome. This approach can be extended to three-dimensional case studies to restore the conditions of natural lava dome growth. The results are published in J. Volcanol. Seismol., 15(3), 1159–168, 2021.

Lava dome morphology inferred from numerical modelling

Lava domes form when highly viscous magmas erupt on the surface. Several types of lava dome morphology can be distinguished depending on the flow rate and the rheology of magma. We develop a two-dimensional axisymmetric model of magma extrusion on the surface and lava dome evolution and analyse the dome morphology using a finite-volume method implemented in Ansys Fluent software. The magma/lava viscosity depends on the volume fraction of crystals and temperature. We show that the morphology of domes is influenced by two parameters: the characteristic time of crystal content growth (CCGT) and the discharge rate (DR). At smaller values of the CCGTs, that is, at rapid lava crystallisation, obelisk-shaped structures develop at low DRs and pancake-shaped structures at high DRs; at longer CCGTs, lava domes feature lobe-shaped to pancake-shaped structures. A thick carapace of about 70% crystal content evolves at smaller CCGTs. We demonstrate that cooling does not play the essential role during a lava dome emplacement, because the thermal thickness of the evolving carapace remains small in comparison with the dome’s height. A transition from the endogenic to exogenic regime of the lava dome growth occurs after a rapid increase in the DR. A strain-rate-dependent lava viscosity leads to a more confined dome, but the influence of this viscosity on the dome morphology is not well pronounced. The model results can be used in assessments of the rates of magma extrusion, the lava viscosity, and the morphology of active lava domes. The results are published in Geophys. J. Int., 223, 1597–1609, 2020.

Examples of morphological shapes of lava domes: on the left are images of domes of the Soufriere Hills volcanic complex (according to Watts et al., 2002); on the right are lava domes according to the studied models (colors indicate the distribution of viscosity in lava domes; the most viscous part of lava in red, the least viscous part in blue).

Nonlinear dynamics of crustal blocks and faults and earthquake occurrences in the Transcaucasian region

The Caucasus is a part of ongoing collision of the Arabian and Eurasian plates, where moderate to strong earthquakes caused significant losses of lives and livelihood in the past. To better understand seismic hazard in the region, we develop a model of block-and-fault dynamics for Transcaucasia, the largest part of the Caucasus to the south of the Greater Caucasus Mountains, to simulate regional earthquakes. The model structure is developed by employing the results of the morphostructural analysis to delineate crustal blocks and the geodetic observations on crustal movements in the region. The model incorporates a nonlinear rate-dependent slip of the faults separating the blocks. A set of numerical experiments has been performed to address the following questions: (i) where strong earthquakes occur and what their reoccurrence time is; (ii) how rigid crustal blocks react to the Arabian plate push and to movements of the ductile part of the crust in Transcaucasia; and (iii) whether the fault slip rates and the block displacements in the model correlate with observed GPS-velocities. The model results confirm that the contemporary crustal dynamics and seismicity pattern in Transcaucasia are determined by the north-northeastern motion of the Arabian plate relative to Eurasia and by the movement of the ductile crust underlying the rigid crustal blocks. Variations in the rheological properties of the fault zones and/or of the ductile crust influence the pattern of seismicity. The number and maximum magnitude of synthetic earthquakes change with the variations in the movements of the crustal blocks and in the rheological properties of the lower crust and the fault zones. The model results can be used in comprehensive seismic hazard assessment of the Caucasus region based on instrumentally observed, historical and synthetic seismicity. The results are published in Physics of the Earth and Planetary Interiors, 297, 106320, 2019.

Distribution of synthetic seismicity (with earthquake magnitudes greater than 6) imposed on the map of observed seismicity: blue stars mark earthquakes of magnitude Ì7+, blue circles Ì6.5+ events, and light blue dots Ì6+ events.

Forging a paradigm shift in disaster science

Despite major advancements in knowledge on disaster risks and disasters caused by natural hazards, the number and severity of disasters are increasing. Convolving natural, engineering, social and behavioral sciences and practices with policymaking should significantly reduce disaster risks caused by natural hazards. To this end, a fundamental change in scientific approaches to disaster risk reduction is needed by shifting the current emphasis on individual hazard and risk assessment dominant in the geoscientific community to a transdisciplinary system analysis with action-oriented research on disaster risk reduction co-produced with multiple stakeholders, including policymakers. This paradigm shift will allow for acquisition of policy-relevant knowledge and its immediate application to evidence-based policy and decision making for disaster risk reduction. The need for the paradigm shift is more critical now than ever before because of the increasing vulnerability and exposure of society to disaster risk and the need for cross-cutting actions in policy and practice related to climate change and sustainability. The results are published in Natural Hazards, 86, 969–988, 2017.

Reconstruction of thermal and dynamic characteristics of volcanic lava from surface thermal measurements

We study a model of volcanic lava flow to determine its thermal and dynamic characteristics from thermal measurements of the lava at its surface. Mathematically this problem is reduced to solving an inverse boundary problem. Namely, using known conditions at one part of the model boundary we determine the missing condition at the remaining part of the boundary. We develop a numerical approach to the mathematical problem in the case of steady-state flow. Assuming that the temperature and the heat flow are prescribed at the upper surface of the model domain, we determine the flow characteristics in the entire model domain using a variational (adjoint) method. We have performed computations of model examples and showed that in the case of smooth input data the lava temperature and the flow velocity can be reconstructed with a high accuracy. As expected, the noise imposed on the smooth input data results in a less accurate solution, but still acceptable below some noise level. Also we analyse the influence of optimization methods on the solution convergence rate. The proposed method for reconstruction of physical parameters of volcanic lava can also be applied to other problems in geophysical fluid flows. The first results are published in Geophys. J. Int., 205, 1767–1779, 2016.

Surface temperature of lava can be extracted from the satellite measurements, e.g. from LANDSAT 7 ETM+ thermal and infrared bands.

Reconstruction of the lava temperature (a) and the flow velocity (c) after 20 and 80 iterations. The relevant residuals of the temperature (b) and the velocity (d) indicate the quality of the reconstruction.

 

Numerical modelling of lava flow with a ruptured crust

Although volcanic lava flows do not significantly affect the life of people, its hazard is not negligible as hot lava kills vegetation, destroys infrastructure, and may trigger a flood due to melting of snow/ice. The lava flow hazard can be reduced if the flow patterns are known, and the complexity of the flow with a debris is analyzed to assist in disaster risk mitigation. In this paper we develop three-dimensional numerical models of a gravitational flow of multi-phase fluid with rafts (mimicking rigid lava-crust fragments) on a horizontal and topographic surfaces to explore the dynamics and the interaction of lava flows. We have obtained various flow patterns and spatial distribution of rafts depending on conditions at the surface of fluid spreading, obstacles on the way of a fluid flow, raft landing scenarios, and the size of rafts. Furthermore, we analyze two numerical models related to specific lava flows: (i) a model of fluid flow with rafts inside an inclined channel, and (ii) a model of fluid flow from a single vent on an artificial topography, when the fluid density, its viscosity, and the effusion rate vary with time. Although the studied models do not account for lava solidification, crust formation, and its rupture, the results of the modeling may be used for understanding of flows with breccias before a significant lava cooling. The first results are published in J. Geodyn., 97, 31-41, 2016.

Spreading lava with a ruptured crust.

A model of lava flow in a channel with ruptured crust (moving rafts are blue, and landed are red).

 

On the use of multiple-site estimations in probabilistic seismic hazard assessment

We analyze specific features of multiple-site (MS) probabilistic seismic hazard assessment (PSHA), i.e. the annual rate of ground-motion level exceedance in at least one site of several sites of interest located within in an area or along a linear extended object. The relation between MS hazard estimates and strong ground-motion data obtained during large earthquakes is discussed in the cases of the 1999 Chi-Chi MW 7.6 earthquake and the 2008 Wenchuan MW 7.9 earthquake. The strong-motion records obtained in the epicentral zones may be considered as an example of the ground motion exceeding the design level estimated using the conventional point-wise (PW) PSHA. We show that the MS-PSHA, when being performed for the standard return period 475 years, provides reasonable estimations of the intensity level that may occur during the earthquakes, parameters of which are close to the parameters of events with maximum possible magnitude accepted in PSHA for the regions. The MS-PSHA efficiency is discussed with respect to (a) evaluation of the performance of probabilistic seismic hazard maps and (b) application of the MS hazard estimates as a basis for design loads. Based on the results of this work, we propose a multi-level approach to PSHA considering fixed reference probability of exceedance (e.g. 10% in 50 years): (i) a standard PW-PSHA to be performed in a seismic-prone region (first level), and (ii) this analysis should be supplemented by a MS-PSHA for urban and industrial areas, or zones of a particular economic and social importance (second level). The results are published in Bull. Seismol. Soc. Am., 106(5), 2233–2243, 2016.

Relation between the multiple-site PGAMLT and point-wise PGAPNT hazard estimations for the 475-yr (a) and 2475-yr (b) return periods. The dashed line shows ratio between the maximum PGA recorded during the Wenchuan earthquake (about 810 cm s-2, geometric mean of two horizontal components) and the estimated design ground motion PGA475 ~ 300 cm s-2 (a) and PGA2475 ~ 600 cm s-2 (b). 1: area 600 km2; 2: 400 km2; 3: 100 km2.

Seismic hazard from instrumentally recorded, historical and simulated earthquakes

We study regional seismic hazard accounting for observed (instrumentally recorded and historic) earthquakes, as well as for seismic events simulated for a significantly longer period of time than that of observations. We use a probabilistic seismic hazards analysis (PSHA) with elements of scenario-based analysis for the Tibet-Himalayan region. The large magnitude synthetic events, which are consistent with the geophysical and geodetic data, together with the observed earthquakes are employed for the Monte-Carlo PSHA. Earthquake scenarios for hazard assessment are generated stochastically to sample the magnitude and spatial distribution of seismicity, as well as the distribution of ground motion for each seismic event. The peak ground acceleration values, which are estimated for the return period of 475 yr, show that the hazard level associated with large events in the Tibet-Himalayan region significantly increases if the long record of simulated seismicity is considered in the PSHA. The magnitude and the source location of the 2008 Wenchuan M=7.9 earthquake are among the range of those described by the seismic source model accepted in our analysis. We analyze the relationship between the ground motion data obtained in the earthquake’s epicentral area and the obtained PSHA estimations using a deaggregation technique. The proposed approach provides a better understanding of ground shaking due to possible large-magnitude events and could be useful for risk assessment, earthquake engineering purposes, and emergency planning. The results were published in Tectonophysics, 657, 187-204, 2015.

Comparison of instrumental intensity (USGS shake map for the 2008 Wenchuan MW 7.9 earthquake) with the PGAs (in cm s-2; for rock site conditions and the return period of 475 yrs) taken from (a) the national Seismic Ground Motion Zonation Map GB 18306-2001, values in parentheses are given for the medium-stiff soil; (b) the GSHAP map; and (c) from the results of this study.

Mantle/lithosphere dynamics beneath the Japanese islands

Recent seismic tomography studies image a low velocity zone (interpreted as a high temperature anomaly) in the mantle beneath the subducting Pacific plate near the Japanese islands at the depth of about 400 km. This thermal feature is rather peculiar in terms of the conventional view of mantle convection and subduction zones. Here we present a dynamic restoration of the thermal state of the mantle beneath this region assimilating geophysical, geodetic, and geological data up to 40 million years. We hypothesise that the hot mantle upwelling beneath the Pacific plate partly penetrated through the subducting plate into the mantle wedge and generated two smaller hot upwellings, which contributed to the rapid subsidence in the basins of the Japan Sea and to back-arc spreading. Another part of the hot mantle migrated upward beneath the Pacific lithosphere, and the presently observed hot anomaly is a remnant part of this mantle upwelling. The results were published in Nat. Sci. Rep., 3, 1137, 2013.

 

Crustal block-and-fault dynamics, fault slip rates and earthquake flow in the Tibet-Himalayan region

The Tibetan plateau and Himalayans have resulted from the continuous Indian and Eurasian plate convergence following their initial collision at about 55 million years ago. Earthquakes in the region occur mainly in response to the crustal motion and stress localization associated with this convergence. To understand the basic features of the motion and seismicity in the Tibet-Himalayan region, I develop a model of its block-and-fault dynamics. The model structure is composed of twelve interacting upper crustal blocks separated by fault planes. I develop several sets of numerical experiments constrained by the regional seismic observations, geodetic and geological data. Synthetic large events in the numerical experiments are clustered mainly on the fault segments associated with the Himalayan Frontal Thrust as well as at some internal faults of the Tibetan plateau. The clustering of earthquakes at some fault is a consequence of dynamics of the regional fault system rather than that of the fault only. I show that variations in the relationship of magnitude to frequency of the events depend on changes in the motion of the upper crustal blocks and on the rheological properties of the lower crust and fault zones. I demonstrate that the predicted crustal motion in the region is characterized by the north-northeastern movement of India toward Eurasia. Fluctuations in rheological properties of fault plane zones and/or the lower crust influence displacement rates of the crustal blocks and hence slip rates at the faults separating the blocks. This can explain the discrepancies in estimates of slip rates at major faults in the region (e.g., Altyn Tagh, Karakorum) over short and long time scales.

Observed and modeled large earthquakes in the Tibet-Himalayan region. (a) Seismicity from 1902 to 2000. White bold lines (model faults) delineate the structural geological elements, and the white arrow indicates the motion of India relative to Eurasia. Earthquake epicenters are marked by colored (depending on the depth of the hypocenters) circles. (b)-(d) Simulated earthquakes for three BAFD catalogs (modified after Ismail-Zadeh et al., 2007).

Quantitative approach to thermo-mechanical restoration of salt-bearing sedimentary basins (SBSBs).

The research activity will address most important processes occurring in SBSBs: thermal and structural evolutions of the basins. In applications they are referred to as a maturity of hydrocarbons, their migration paths, deformations of sedimentary layers, and stresses regime. I intend to study thermo-mechanical restoration models for Pricaspian salt basin (East European platform), which plays an important role in hydrocarbon explorations and exploitations. The goal of the intended research is to understand the interplay between geodynamic, geothermal, and tectonic processes involved in the evolution of the SBSB. The specific objectives of the research are (i) to develop a 3D numerical approach to structural restoration of SBSBs; (ii) to develop a 3D numerical approach to thermo-mechanical restoration of SBSBs, and hence to determine temperatures within the basins in the geological past based on their present-day estimations; and (iii) to analyze the evolution of the Pricaspian basin using thermo-mechanical restoration models and geological and geophysical observations.

Inverse Problem of the Thermal Convection

To restore quantitatively the seismically detected mantle structures and temperature field, we need a numerical tool for solving an inverse problem of thermal convection at infinite Prandtl number. I develop a variational approach to three-dimensional numerical restoration of thermoconvective mantle flow. The approach is based on a search for the mantle temperature and flow in the geological past by minimizing differences between mantle temperature derived from seismic velocities (or their anomalies) and temperature predicted by forward models of mantle flow. The mantle temperatures and flow fields in the past so obtained could be employed as constraints on forward models of mantle dynamics.

Salt Tectonics and Basin Evolution

Numerical studies of ductile deformations induced by salt movements have, until now, been restricted to two-dimensional modeling (2D) of diapirism. I study several 3D models of salt diapirism and deformation of overlying sediments. We analyze a model of salt diapirs that develop from an initial random perturbation of the interface between salt and its overburden and then restore the evolved salt diapirs to their initial stages. An evolutionary model of a 2D salt wall loaded by a 2D clinoform of sediments predicts a decomposition of the salt wall into 3D diapiric structures when the overburden of salt is supplied by 3D syn-kinematic wedge of sediments. Also we study how lateral flow effects the evolution of salt diapirs. The shape of a salt diapir can be very different, if the rate of horizontal flow is much greater than the initial rate of diapiric growth solely due to gravity.

Geodynamics, Stress and Seismicity

Tectonic stress in the lithospheric slab descending beneath the Vrancea (SE-Carpathians)

To understand processes of stress generation in the descending slab, I analyze tectonic stress in the slab by means of analytical and numerical modeling. I have found that the maximum shear stress migrates from the upper plane of the Benioff-Wadati seismic zone to its lower plane in course of changes in slab dynamics from its active subduction through roll-back movements to passive sinking solely due to gravity. It can explain a location of hypocenters of Vrancea events at the side of slab adjacent to the Eastern European platform. To understand processes of stress generation, I analyze the tectonic stress induced by the sinking Vrancea slab employing 3D Eulerian FEM and using (as input model data) seismic tomographic high-resolution images of the Vrancea lithosphere, results of refraction seismic study, and other geophysical and geological data. On the basis of experimental data on elastic parameters and anelasticity I obtain initially a model of mantle temperature from the P-wave velocity anomalies; crustal temperature is derived from heat flow data. Based on temperature- and pressure-dependent viscosity and density, modeled tectonic stress predicts horizontal compression beneath the Vrancea region, which coincides with the stress regime defined from fault-plane solutions for intermediate-depth earthquakes. The stress reaches maximum at depths of 70 to 110 km and 130 to 180 km and decays below being in a good agreement with the observed seismicity.

Stress and buoyancy-driven deformations of the lithosphere beneath central Apennines

The juxtaposed contraction and extension, a long-standing geological enigma recognized in different geodynamic frameworks world-wide, observed at the surface and active nowadays, hence better studied, only in the Italian Peninsula, has for a long time attracted the attention of geoscientists. Several models, invoking mainly external forces, have been put forward to explain the close association of these two end-member deformation mechanisms clearly observed by geophysical, geodetic and geological investigations. These models appeal to interactions along plate margins or at the base of the lithosphere such as back-arc extension or shear tractions from mantle flow or to subduction processes such as slab roll back, retreat or pull and detachment. On the basis of seismic images of the lithosphere beneath the central Apennines and numerical modeling of lithosphere dynamics and in-situ stress, I find that buoyancy forces can explain solely the present-day stress distribution in the Apennine lithosphere, complex lithospheric deformations, and unusual earthquake distribution.

Dynamic quantitative restoration of profiles across diapiric structures

The backstripping method that is widely used in basin analysis sometimes fails for salt-bearing basins, because the highly mobile and buoyant salt deforms its sedimentary overburden. I developed a new numerical approach for two-dimensional dynamic restoration of cross-sections through successive earlier depositional stages. Our models interpreted basin profiles as multiple ductile layers with various physical parameters (e.g., density, viscosity, Young’s modulus). The evolution of salt structures was modeled backwards in time by removing successively younger layers and restoring older layers and any diapirs to the stage they were likely to have been. The applicability of the technique was demonstrated by reconstructions of upbuilt and downbuilt diapirs. I used the technique to restore a depth-converted seismic cross-section through the south-eastern part of the Pricaspian salt basin [5]. Mature salt diapirs in the section are shown to have been downbuilt from a salt layer with an initially uniform thickness, as a result of differential sedimentary loading until the end of the Triassic before one of the diapirs was buried and actively upbuilt. Our approach is well suited for restoration of cross-sections with ductile overburdens and now is being developed to three-dimensional restorations and other rheologies.

Gravitational and buckling instabilities of rheologically stratified geostructures

I investigated the gravitational and buckling instabilities of a structure consisting of a buoyant layer of viscous fluid overlain by a dense perfectly plastic layer. The structure is subject to either horizontal extension or shortening and models rocksalt under a brittle overburden. Considering the viscosity of the buoyant layer to be much less than the effective viscosity of the overlying layer, we obtained the following results. The instability pattern of the structure is similar to that of a perfectly plastic structure. The characteristic wavelength, corresponding to the most unstable mode, increases initially with the thickness ratio between the lower and upper layers, but then decreases by a series of abrupt jumps. This is in contrast to the result that the characteristic wavelength approaches to a constant value with the increasing thickness ratio in the case of viscous layers. The buckling instability induced by rapid horizontal straining overwhelms the gravitational instability, and the growth rate of the instability depends linearly on the effective viscosity ratio. We studied models of diapirism in the Great Kavir, Iran. We show that a small interdiapir spacing in the salt canopy province and a wide range of the spacings in the salt pillow province of the region can be explained by the perfectly plastic sedimentary overburden and horizontal shortening.

Stress in descending lithospheric slabs

I examined the effects of viscous flow, phase transition, and dehydration on stress in a relic slab to explain the intermediate-depth seismic activity in the Vrancea region. I developed a two-dimensional finite-element model of a slab gravitationally sinking in the mantle which predicts (i) downward extension in the slab as inferred from the stress axes of earthquakes, (ii) the maximum stress occurring in the depth range of 70 km to 160 km, and (iii) a very narrow area of the maximum stress. The depth distribution of the annual average seismic energy released in earthquakes has a shape similar to that of the depth distribution of the stress in the slab. Estimations of the cumulative annual seismic moment observed and associated with the volume change due to the basalt-eclogite phase changes in the oceanic slab indicate that a pure phase-transition model cannot explain solely the intermediate-depth earthquakes in the region. We consider that one of realistic mechanisms for triggering these events in the Vrancea slab can be the dehydration of rocks which makes fluid-assisted faulting possible. The approach can be applied to other regions of intermediate-depth seismicity.

Block-and-fault dynamics of the lithosphere and earthquakes in the SE-Carpathians

I studied a numerical model of block-and-fault dynamics of the lithosphere beneath the earthquake-prone Vrancea region. The model presents a system of absolutely rigid blocks separated by infinitely thin plane faults. The interaction of the blocks along the fault planes and with the surrounding medium is assumed to be a viscoelastic. Motions of boundary blocks cause the displacements of the block system. The velocities of the motions are found from a model of mantle flows induced by a sinking slab beneath the Vrancea region. When a ratio of stress to pressure for some portion of a fault plane exceeds a certain strength level, a stress-drop ('earthquake') occurs. As a result of the numerical simulations catalogues of synthetic earthquakes are produced. Several numerical experiments for various model parameters illustrate that the spatial distribution of synthetic events is significantly sensitive to the directions of the block movements. Small variations in a slab rotation control the pattern of the synthetic seismicity. The results of the analysis indicate that the catalogues obtained by the simulation of the block structure dynamics have certain features similar to those of the real earthquake catalogue of the Vrancea region. The results are in a good agreement with recent seismic tomography data.

Computational Methodology to Study Geodynamic Problems

I developed a new two-dimensional Eulerian finite element methodology to study geodynamical problems where chemical composition changes abruptly across interfaces. The method combines a Galerkin-spline technique with a method of integration over regions bounded by advected interfaces to represent discontinuous variations of material parameters. It allows to directly approximate a natural free surface position, instead of a posteriori calculation of topography from the normal stress at the upper free-slip boundary. The suggested approach avoids the difficulties due to material discontinuities at intermediate boundaries, like Moho or the earth's surface, and is also free from the deficiencies of the Lagrangian approach always resulting in mesh distortion. Also I developed new three-dimensional methods for study of lithospheric deformation and mantle convection. Special computer codes were written for parallel supercomputers with distributed memory.