Papers by Yury Podladchikov
PROCEEDINGS OF THE INTERNATIONAL CONFERENCE “PHYSICAL MESOMECHANICS. MATERIALS WITH MULTILEVEL HIERARCHICAL STRUCTURE AND INTELLIGENT MANUFACTURING TECHNOLOGY”
A paper considers mathematical model, numerical method and its software implementation on GPU for... more A paper considers mathematical model, numerical method and its software implementation on GPU for the simulation of the development and localization of plastic shear bands in a poroelastoplastic solid. A mathematical model of the process is described by a system of nonlinear partial differential equations of poroelastoplasticity, generalizing the Biot model for a two-phase fluid-saturated poroelastoplastic medium. A nonlinear dependence of material parameters on a dynamic porosity, which, in turn, depends on the volumetric deformation of the skeleton, is taken into account. For the numerical solution of the posed problem an isoparametric spectral element method is used to discretize a geometric model and poroelastoplastic equations in space on curvilinear unstructured meshes. A dynamic relaxation method is used to obtain a stationary solution of the boundary value problem using an explicit time scheme, which leads to a high scalability of the computational process when parallelizing the algorithm on massively parallel GPU. The algorithm is implemented using CUDA technology. A spectral element mesh is naturally mapped onto the CUDA Grid, and accordingly, each spectral element is mapped onto a CUDA Block, so local nodes inside a spectral element are processed by the threads in the Block. This approach allows efficient use of the shared memory for data caching in the Block and significantly increase the throughput of the parallel algorithm, which is essentially memory bound. Numerical results of the simulation of the development and localization of plastic shear bands nearby a borehole drilled in a porous fluid saturated solid are presented. Evolution of porosity and permeability as a result of the accumulation of plastic deformations is analyzed.
Fluid flow processes below the Earth's surface play an important role in many geological settings... more Fluid flow processes below the Earth's surface play an important role in many geological settings, from shallow depths in aquifers and hydrocarbon reservoirs down to tens of kilometers, where subduction of oceanic crust brings fluids to large depths. Industrial applications, environmental studies as well as fundamental geological research depends on accurate description and understanding of fluid flow, mass transport and reaction in rocks. Numerical modeling of fluid flow and reaction through rocks is a powerful tool to understand and interpret geological field observations as well as for predicting large scale geodynamic processes, hydrocarbon migration, CO 2 storage reliability, formation of economic ore deposits, and subsurface contamination, to name just a few. In all these cases it is important to have powerful predictive modeling tools that can run in sufficiently high resolution to capture the physical processes on a large variety of scales, both in time and space. Chemical reactions and fluid flow are intimately linked and have to be solved in a coupled way to accurately and consistently model the geological processes. To date, the usual approach in most reactive transport modelling consists of solving non-linear system of chemical reactions at each grid node in the model domain. Here, we explore the use of Gibbs minimization for solving chemical reactions and fluid flow in a coupled system. We apply our method to dehydration of serpentinite during subduction processes. Dynamic porosity, density, and fluid pressure evolution are coupled to reactions by precomputing tens of thousands of phase diagrams covering P-T and compositional space. We compare the results of two approaches: solving individual reactions versus Gibbs minimization using precomputed look-up tables. Both approaches give similar results, however, the Gibbs minimization approach allows for including most complex solid and fluid solution models, whereas in the reaction approach this is very cumbersome. Additionally, we explore the relation between activity/chemical potential diagrams, and isochemical phase diagram sections and evaluate their usage in reactive transport models.
Journal of Geophysical Research: Solid Earth, 2021
Majority of the most powerful supercomputers on the world host hardware accelerators to sustain c... more Majority of the most powerful supercomputers on the world host hardware accelerators to sustain calculations at the petascale level and beyond. Graphical processing units (GPUs) are amongst widely employed hardware accelerators, initiating a revolution in high-performance computing (HPC) in the last decade. The three-dimensional calculations targeting billions of grid cells-technically impossible resolutions decades ago-became reality. This major breakthrough in HPC and supercomputing comes however with the cost of developing and reengineering scientific codes to efficiently utilize the available computing power. Increasing the low-level parallelism is the key. In Earth sciences, HPC and GPU-accelerated applications target in particular forward and inverse seismic modeling and geodynamics-fields where high spatial and temporal resolutions as well as large spatial domains are required. We here develop a multi-GPU implementation for applications in seismic modeling in porous media. Understanding seismic wave propagation in fluid-saturated porous media enables more accurate interpretation of seismic signals in Earth sciences. The two phase medium is represented by an elastic solid matrix (skeleton) saturated with a compressible viscous fluid. The dynamic response of such an isotropic two phase medium results in two longitudinal waves and one shear wave, as predicted by Frenkel (1944) (see also Pride and Garambois (2005)). The wave of the first kind (fast wave) is a true longitudinal wave where the solid matrix motion and the fluid velocity vector fields are in-phase. The wave of the second kind (slow wave) is a highly attenuated wave where the solid matrix motion and the fluid velocity vector fields are
Poromechanics VI, 2017
Experimental techniques have been developed to measure two poroelastic parameters: unjacketed bul... more Experimental techniques have been developed to measure two poroelastic parameters: unjacketed bulk modulus K s ′ and unjacketed pore modulus K s ". K s ′ is measured on samples with no jacket that are instrumented with strain gages, confining fluid (oil) penetrates inside the material and hydrostatic mean stress is equal to pore pressure, so the overal unjacketed bulk modulus can be calculated. K s " is measured on water-saturated rock by determining the change in pore volume of the specimen from accurate measurements on pore pressure controller, while the difference between mean stress and pore pressure is preserved constant. K s ′ and K s " are measured to be pressure-independent for Berea sandstone and Opalinus clay (shale). K s ′ ≈ K s " for the clay-rich shale, but K s " is smaller than K s ′ and bulk modulus of quartz for quartz-rich Berea sandstone. This result is explained by the presence of non-connected porosity.
Earth and Planetary Science Letters, 2021
Solid volume increases during the hydration of mantle rocks. The fluid pathways necessary to feed... more Solid volume increases during the hydration of mantle rocks. The fluid pathways necessary to feed the reaction front with water can be filled by the low density reaction products. As a result, the reaction front dries out and the reaction stops at low reaction progress. This process of porosity clogging is generally predicted to dominate in reactive transport models, even when processes such as reaction-induced fracturing are considered. These predictions are not consistent with observations at mid-ocean ridges where dense mantle rocks can be completely replaced by low density serpentine minerals. To solve this issue, we develop a numerical model coupling reaction, fluid flow and deformation. High extents of reaction can only be achieved when considering that the increase in solid volume during reaction is accommodated through deformation rather than porosity clogging. The model can generate an overpressure that depends on the extent of reaction and on the boundary conditions. This overpressure induces viscoelastic compaction that limits the extent of the reaction. The serpentinisation rate is therefore controlled by the accommodation of volume change during reaction, and thus by deformation, either induced by the reaction itself or by tectonic processes.
non-ideality may lead to a different sign in the molar formulation at the specific conditions. Th... more non-ideality may lead to a different sign in the molar formulation at the specific conditions. These effects need to be explored in future applications. Data and code availability Details on modelling and the approach are included in the main text. Additional material, as referred in the main text, is available in the Supplementary material. Further correspondence and requests for materials and codes should be addressed to L.T.
Geology, 2020
The densification of the lower crust in collision and subduction zones plays a key role in shapin... more The densification of the lower crust in collision and subduction zones plays a key role in shaping the Earth by modifying the buoyancy forces acting at convergent boundaries. It takes place through mineralogical reactions, which are kinetically favored by the presence of fluids. Earthquakes may generate faults serving as fluid pathways, but the influence of reactions on the generation of seismicity at depth is still poorly constrained. Here we present new petrological data and numerical models to show that in the presence of fluids, densification reactions can occur very fast, on the order of weeks, and consume fluids injected during an earthquake, which leads to porosity formation and fluid pressure drop by several hundreds of megapascals. This generates a mechanically highly unstable system subject to collapse and further seismic-wave emission during aftershocks. This mechanism creates new pathways for subsequently arriving fluids, and thus provides a route for self-sustained dens...
The generation of a tsunami by a landslide is a complex phenomenon that involves landslide dynami... more The generation of a tsunami by a landslide is a complex phenomenon that involves landslide dynamics, wave dynamics and their interaction. Numerous lives and infrastructures around the world are threatened by this phenomenon. Predictive numerical models are a suitable tool to assess this natural hazard. However, the complexity of this phenomenon causes such models to be either computationally inefficient or unable to handle the overall process. Our model, which is based on shallow-water equations, has been developed to address these two problems. In our model, the two materials are treated as two different layers, and their interaction is resolved by momentum transfer inspired by elastic collision principles. The goal of this study is to demonstrate the validity of our model through benchmark tests based on physical experiments performed by Miller et al. (2017). A dry case is reproduced to validate the behaviour of the landslide propagation model using different rheological laws and to determine which law performs best. In addition, a wet case is reproduced to investigate the influence of different still-water levels on both the landslide deposit and the generated waves. The numerical results are in good agreement with the physical experiments, thereby confirming the validity of our model, particularly concerning the novel momentum transfer approach.
Geophysical Research Letters, 2015
ABSTRACT The migration of gases from deep to shallow reservoirs can cause damageable events. For ... more ABSTRACT The migration of gases from deep to shallow reservoirs can cause damageable events. For instance, some gases can pollute the biosphere or trigger explosions and eruptions. Seismic tomography may be employed to map the accumulation of subsurface bubble-bearing fluids to help mitigating such hazards. Nevertheless, how gas bubbles modify seismic waves is still unclear. We show that saturated rocks strongly attenuate seismic waves when gas bubbles occupy part of the pore space. Laboratory measurements of elastic wave attenuation at frequencies <100 Hz are modeled with a dynamic gas dissolution theory demonstrating that the observed frequency-dependent attenuation is caused by wave-induced-gas-exsolution-dissolution (WIGED). This result is incorporated into a numerical model simulating the propagation of seismic waves in a subsurface domain containing CO2-gas bubbles. This simulation shows that WIGED can significantly modify the wavefield and illustrates how accounting for this physical mechanism can potentially improve the monitoring and surveying of gas bubble-bearing fluids in the subsurface.
Tectonophysics, 1993
Experiments to study the behavior of various materials point to the relation that exists between ... more Experiments to study the behavior of various materials point to the relation that exists between elastic properties and the type of stress. The influence of the state of stress on the elasticity of a fractured material will be discussed for a physically non-linear model of an elastic solid. The strain-dependent moduli model of material, presented in this paper, makes it possible to describe this feature of a solid. ft also permits to simulate a dilatancy of rocks. A damage parameter, introduced into the model using a the~o~amical approach, allows to describe a rheologicai transition from the ductile regime to the brittle one, and to simulate the rock's memory, narrow fracture zone creation and strain rate localization. Additionally, the model enables the investigation of the final geometry of fracture zones, and also to simulate their creation process, taking into account pre-existing fracture zones. The process of narrow fracture zone creation and strain rate localization was simulated numerically for single axis compression and shear flow.
Tectonophysics, 2000
In a rheologically layered crust, compositional layers have an upper, elasto-plastic part and a l... more In a rheologically layered crust, compositional layers have an upper, elasto-plastic part and a lower, viscous one. When broken, the upper elastic part undergoes flexure, which is upward for the foot-wall and downward for the hanging wall. As a consequence of bending, stresses will develop locally that can overcome the strength of the plate and, therefore, impose the migration of active fault. In the lower, viscous part of each compositional layer, rocks can potentially flow. Numerical modelling of the behaviour of a crust made up of two compositional layers, during and following extension, shows that flow can take place not only in the lower crust but also, and more importantly, in the lower part of the upper crust. The ability of crustal rocks to flow influences the style and kinematics of rifted regions. When no flow occurs, subsidence will affect the extending areas, both hanging wall and foot-wall will subside with respect to an absolute reference frame such as sea level, and there will be a strict proportionality between extension and thinning. In addition, the downward movement of fault blocks will decrease the local stresses created in the foot-wall and increase those of the hanging wall, thereby imposing a migration of the active fault towards the hanging wall. This is the behaviour of extensional settings developed on stabilised crust and which evolved in a passive margin. When flow does take place, middle crustal rocks will move towards the rifting zone causing isostatically driven upward movements that will be superimposed on movements associated with crustal and lithospheric thinning. Consequently, fault blocks will move upwards and the crust will show more extension than thinning. The upward movements will decrease the stresses developed in the hanging walls and increase those of the foot-wall. Faults will then migrate towards the foot-wall. Such a mode of deformation is expected in regions with thickened crust and has its most apparent expression in core complexes.
Петрология, 2013
We review published evidence that rocks can develop, sustain and record significant pressure devi... more We review published evidence that rocks can develop, sustain and record significant pressure devi ations from lithostatic values. Spectroscopic studies at room pressure and temperature (P-T) reveal that in situ pressure variations in minerals can reach GPa levels. Rise of confined pressure leads to higher amplitude of these variations documented by the preservation of α quartz incipiently amorphized under pressure (IAUP quartz), which requires over 12 GPa pressure variations at the grain scale. Formation of coesite in rock deformation experiments at lower than expected confined pressures confirmed the presence of GPa level pressure variations at elevated temperatures and pressures within deforming and reacting multi mineral and polycrystalline rock samples. Whiteschists containing garnet porphyroblasts formed during prograde metamorphism that host quartz inclusions in their cores and coesite inclusions in their rims imply preserva tion of large differences in pressure at elevated pressure and temperature. Formation and preservation of coherent cryptoperthite exsolution lamellae in natural alkali feldspar provides direct evidence for grain scale, GPa level stress variations at 680°C at geologic time scales from peak to ambient P T conditions. Similarly, but in a more indirect way, the universally accepted 'pressure-vessel' model to explain preservation of coesite, diamond and other ultra high pressure indicators requires GPa level pressure differences between the inclu sion and the host during decompression at temperatures sufficiently high for these minerals to transform into their lower pressure polymorphs even at laboratory time scales. A variety of mechanisms can explain the for mation and preservation of pressure variations at various length scales. These mechanisms may double the pressure value compared to the lithostatic in compressional settings, and pressures up to two times the litho static value were estimated under special mechanical conditions. We conclude, based on these consider ations, that geodynamic scenarios involving very deep subduction processes with subsequent very rapid exhu mation from a great depth must be viewed with due caution when one seeks to explain the presence of micro scopic ultrahigh pressure mineralogical indicators in rocks. Non lithostatic interpretation of high pressure indicators may potentially resolve long lasting geological conundrums.
Nature Geoscience, 2009
Intermediate-depth (50-300 km) earthquakes commonly occur along convergent plate margins but thei... more Intermediate-depth (50-300 km) earthquakes commonly occur along convergent plate margins but their causes remain unclear. In the absence of pore-fluid pressures that are sufficiently high to counter the confining pressure in such settings, brittle failure is unlikely. In such conditions, the rocks could fail by the mechanism of progressively self-localizing thermal runaway 1 , whereby ductile deformation in shear zones leads to heating, thermal softening and weakening of rock 1-3. Here we test this hypothesis by focusing on fault veins of glassy rock (pseudotachylyte) formed by fast melting during a seismic event, as well as associated ductile shear zones that occur in a Precambrian terrane in Norway. Our field observations suggest that the pseudotachylytes as well as shear zones have a single-event deformation history, and we also document mineralogical evidence for interaction of the rocks with external fluids. Using fully coupled thermal and viscoelastic models, we demonstrate that the simultaneous occurrence of brittle and ductile deformation patterns observed in the field can be explained by self-localizing thermal runaway at differential stresses lower than those required for brittle failure. Our results suggest that by perturbing rock properties, weakening by hydration also plays a key role in shear zone formation and seismic failure; however, thermal runaway enables the rocks to fail in the absence of a free fluid phase. The intense seismicity in subduction zones at depths greater than 50 km has puzzled geoscientists for decades because high ambient pressures should inhibit brittle failure at such depths. Two main failure mechanisms for intermediate-depth earthquakes have been hypothesized: (1) dehydration embrittlement and (2) shear instabilities. Metamorphic dehydration reactions have the potential to raise the pore fluid pressure and thereby lower the effective pressure to values that permit brittle failure 4-6. Elastic energy stored in a viscoelastic material may be spontaneously released at seismic strain rates by the formation of very high-temperature self-localizing shear instabilities 1,3,7 , and the most unknown critical factor for this failure mode is the conditions under which an initial perturbation of the system either decays or self-amplifies, leading to extreme localization of deformation and shear heating. Both mechanisms are theoretically viable but it is difficult to discriminate between them. Pseudotachylyte fault veins (quenched melts) provide evidence for paleo-earthquakes 8 , and they can be used to study the details of seismic failure. It has been suggested that pseudotachylytes associated with shallow earthquakes (10-15 km depth) are usually formed by brittle failure followed by frictional melting 9 , but at low background temperatures and high stresses ductile deformation
Mineralogy and Petrology, 2009
Magmatic enclaves from the Rudolfov quarry near Liberec (Czech Republic) are interpreted to repre... more Magmatic enclaves from the Rudolfov quarry near Liberec (Czech Republic) are interpreted to represent remnants of lamprophyric melt that intruded the Karkonosze granite at a stage at which the granite was not fully solidified. Based on the observation that larger enclaves are generally more circular than the smaller ones, we conclude that bigger blobs of mafic magma became more spherical during flow in the gravity field (sink or float). This flow is also interpreted to be responsible for the incorporation of mineral grains into the enclaves and may have facilitated the assimilation of granitic melt. Linear mixing trends on Harker diagrams for most network-forming and mainly slow-diffusing or fluidimmobile elements indicate such an assimilation process between granite and lamprophyre. In contrast, all fastdiffusing or fluid-mobile elements display scattered element distributions, implying that chemical diffusion also played a role. Pressure and temperature for this late magmatic stage are estimated at around 1 kbar and 500°C. These results suggest that two processes modified the composition of the enclaves in the Karkonosze granite: (1) assimilation (mechanical mixing) of granitic melt during the injection of the lamprophyric melt and the subsequent flow of the forming enclaves in the gravity field (responsible for the linear mixing trends) and (2) diffusion-controlled redistribution of elements between both solidifying rock types at the magmatic stage and/or due to late-stage magmatic fluids (responsible for the scattering and deviation from the linear mixing trends).
Journal of Metamorphic Geology, 2009
Chemical reactions and phase changes generally involve volume changes. In confined settings this ... more Chemical reactions and phase changes generally involve volume changes. In confined settings this will cause a mechanical deformation of the matrix that surrounds the reaction sites where the volume change takes place. Consequently, mineral reactions and the mechanical response of the rock matrix are coupled. A companion paper in this issue illustrates this coupling with experiments where quartz and olivine react to form enstatite reaction rims under ambient conditions of 1 GPa and 1000°C. It has been demonstrated that for identical run conditions, the thickness of the reaction rims depends on whether quartz grains are embedded in an olivine matrix or olivine grains are included in a quartz matrix. The experimental conditions, the nature of the results, and the large volume change of the reaction ()6%) leave only viscous creep as a viable matrix response to the reaction progress. A model is developed for this reaction, which combines diffusion of chemical components through the growing rim and viscous creep of the matrix. The resulting rate law for reaction rim growth in spherical geometry shows that the progress rate is proportional to the reaction overstepping and controlled by the slower of the two competing processes; either diffusion or creep. If diffusion is rate limiting the usual linear proportionality between rim growth and ffiffi t p results. However, if viscous creep is rate limiting, then the reaction rates are reduced and may become effectively creep controlled. With respect to the experiments in the companion paper it is inferred that the effective viscosity of the two matrix materials, i.e. polycrystalline quartz and olivine, differ by approximately one order of magnitude with the quartz being the stronger one. The absolute values of the inferred viscosities correspond well to published flow laws. The rheological properties of natural rocks are well within the parameter range for which significant mechanical control on reaction rim growth is expected. This implies that for the interpretation of natural reaction rims and corona structures both diffusion and mechanical control need to be considered. In addition the mechanical effect also needs to be considered when interdiffusion coefficients are retrieved from rim growth experiments. This should also be considered for geospeedometry analyses. Furthermore, the control on reaction rate because of slow creep of the matrix is expected to be even more important, compared to the experiments, under colder crustal conditions and may contribute substantially to the frequent observation of only partially completed reactions. We suggest that this phenomenon is referred to as Ômechanical closureÕ, which may be an important mechanism in the kinetic displacement of the boundaries between the stability fields of phase assemblages.
Journal of Geophysical Research, 2010
Piercement structures such as mud volcanoes, hydrothermal vent complexes, pockmarks and kimberlit... more Piercement structures such as mud volcanoes, hydrothermal vent complexes, pockmarks and kimberlite pipes, form during the release of pressurized fluids. The goal of this work is to predict under which conditions piercement structures form from the insights gained by sand box experiments injecting compressed air through an inlet of width w at the base of a bed of glass beads of height h. At an imposed critical velocity v f , a fluidized zone consisting of a diverging cone-like structure formed with morphological similarities to those observed in nature. Dimensional analysis showed that v f is correlated to the ratio of h over w. In addition, we derived an analytical model for v f which is compared to the experimental data. The model consists of a force balance between the weight and the seepage forces imparted to the bed by the flowing gas. The analytic model reproduces the observed correlation between v f and h/w, although a slight underestimate was obtained. The results suggest that the gas-particle seepage force is the main triggering factor for fluidization and that the commonly used proxy, which the fluid pressure must equal or exceed the lithostatic weight, needs to be reconsidered. By combining the experiments and the model, we derived critical pressure estimates which were employed to a variety of geological environments. Comparing the estimated and measured pressures prior to the Lusi mud volcano shows that the presented model overestimates the critical pressures. The model paves the way for further investigations of the critical conditions for fluidization in Earth systems.
Journal of Geophysical Research, 2009
Water-, mud-, gas-, and petroleum-bearing seeps are part of the Salton Sea geothermal system (SSG... more Water-, mud-, gas-, and petroleum-bearing seeps are part of the Salton Sea geothermal system (SSGS) in southern California. Carbon dioxide is the main component behind the seeps in the Davis-Schrimpf seep field ($20,000 m 2). In order to understand the mechanisms driving the system, we have investigated the seep dynamics of the field by monitoring the temperature of two pools and two gryphons for 2180 h (90.8 days) in the period from December 2006 to March 2007, with a total of 32,700 measurements per station. The time series have been analyzed by statistical methods using cross correlation, autocorrelation and spectral analysis, and autoregressive modeling. The water-rich pools never exceed 34.0°C and are characterized by low-amplitude temperature variations controlled by the diurnal cycles in air temperature. The long-term validity of these results is evident from a second period of temperature monitoring of one of the pools from December 2007 to April 2008 (120 days). In contrast to the pools, the mud-rich gryphons have a strikingly different behavior. The gryphons are hotter (maximum 69.7°C) and have large amplitude variations (standard deviation of 6.4) that overprint any signal from external diurnal forcing. Autoregressive modeling shows the presence of distinct hot and cold pulses in the gryphon temperature time series, with amplitudes up to 3°C. These pulses likely reflect a combination of hydrothermal flux variations from the SSGS and the local temporal changes in bubbling activity within the gryphons.
Uploads
Papers by Yury Podladchikov