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.



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.