EMRP – Earth Magnetism & Rock Physics

EMRP1.14 – Multiscale rock damage in geology, geophysics and geo-engineering systems

EGU21-15691 | vPICO presentations | EMRP1.14

Microfracture evolution leading to catastrophic failure observed by hearing and seeing

Maria-Daphne Mangriotis, Andrew Curtis, Alexis Cartwright-Taylor, Edward Andò, Ian Main, Andrew Bell, Ian B. Butler, Florian Fusseis, Roberto Rizzo, Martin Ling, Sina Marti, Derek Doug Vick Leung, Jonathan Singh, and Oxana Magdysyuk

Catastrophic failure is a critical phenomenon present in Earth systems on a variety of scales, and is associated with the evolution of damage leading to system-size failure. Laboratory testing of rock failure permits characterization of fracture network evolution at the micro-scale to understand the interaction of cracks, pores and grain boundaries to an applied stress field, and the relationship between deformation and seismic response. Previous studies have relied on acoustic emissions (hearing) or X-ray imaging (seeing) to study the process of localization, which involves spontaneous self-organization of smaller cracks along faults and fractures on localised zones of deformation. To combine hearing and seeing of the microscopic processes and their control of system-sized failure, a novel x-ray transparent cell was used for deformation experiments of rock samples, which permits integration of acoustic monitoring with fast synchrotron x-ray imaging. To increase temporal characterization of damage beyond the temporal resolution of the fast 3D synchrotron system, acoustic emission (AE) feedback control was used to regulate the applied stress and slow down the deformation processes. As a result, there is increased temporal resolution of the incremental deformation between successive x-ray scanned states allowing synchronized comparison of acoustic emissions to x-ray scans. Here, we present the seismic analysis used to characterize the velocity evolution of the rock samples, and the location and characteristics of individual AE events in relation to microscopic deformation processes. Time-lapse velocity measurements are linked to internal stress changes and structural damage corresponding to seismic and aseismic deformation processes, while acoustic emissions are a direct indication of local cracking.  We show that we can successfully locate AE events in 3D using only two sensors on either end of the sample, based on ellipsoid mapping, and x-ray image to event correlation. We explore temporal and spatial statistics of AE signatures and how those are linked to the strain field in the samples measured with incremental digital volume correlation between pairs of recorded x-ray tomograms. The direct observation of AE and X-ray images enables quantification of relevant seismic (local cracking leading to AE) to aseismic (elastic loading and silent irreversible damage) processes, with information extracted over fine temporal resolution throughout the deformation process through the AE-feedback control.

How to cite: Mangriotis, M.-D., Curtis, A., Cartwright-Taylor, A., Andò, E., Main, I., Bell, A., Butler, I. B., Fusseis, F., Rizzo, R., Ling, M., Marti, S., Leung, D. D. V., Singh, J., and Magdysyuk, O.: Microfracture evolution leading to catastrophic failure observed by hearing and seeing, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15691,, 2021.

EGU21-2629 | vPICO presentations | EMRP1.14

From monodisperse to polydisperse: the influence of grain size distribution on the mechanical behavior of porous synthetic rocks

Lucille Carbillet, Michael Heap, Fabian Wadsworth, Patrick Baud, and Thierry Reuschlé

Sedimentary crustal porous rocks span a wide range of grain size distributions – from monodisperse to highly polydisperse. The distribution of grain size depends on the location and conditions of rock formation, the chemico-physical processes at play, and is influenced by subsequent geological processes. Well-sorted granular rocks, with a grain size distribution close to monodisperse, and granular rocks with a more polydisperse grain size distribution, have repeatedly been subjected to laboratory experiments. And yet the natural variability from sample to sample and structural heterogeneity within single natural samples all conspire to prevent us from constraining the effect of grain size polydispersivity. While a few studies have focused on the influence of grain size, the control of grain size distribution on the mechanical behavior of rocks has scarcely been studied, especially in the laboratory. In this study, we address this knowledge-gap using synthetic samples prepared by sintering glass beads with controlled polydisperse grain size distributions. When heated above the glass transition temperature, the beads act as viscous droplets and sinter together. Throughout viscous sintering, a bead pack evolves from an initial granular discontinuous state into a solid connected porous state, at which the microstructural geometries and final porosity are known. Variably polydisperse individual samples were prepared by mixing glass beads with diameters of 0.2, 0.5, and 1.15 mm in various proportions, which were sintered together to a final porosity of 0.25 or 0.35. Hydrostatic and triaxial compression experiments were performed for each combination of polydispersivity. The samples were water-saturated, deformed at room temperature, and deformed under drained conditions (with a fixed pore pressure of 10 MPa). Triaxial experiments were conducted at a constant strain rate at effective pressure corresponding to the ductile (compactive) regime. Our mechanical data provide evidence that polydispersivity exerts a significant control on the compactive behavior of porous rocks. Insights into the microstructure were gained using scanning electron microscopy on thin sections prepared from samples before and after deformation. These data allow for the observation of the different deformation features, and by extension the deformation micro-mechanisms, promoted by the different type and degree of polydispersivity. Overall, our data show that, at a fixed porosity, increasing polydispersivity decreases the stress required for compactant failure.

How to cite: Carbillet, L., Heap, M., Wadsworth, F., Baud, P., and Reuschlé, T.: From monodisperse to polydisperse: the influence of grain size distribution on the mechanical behavior of porous synthetic rocks, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-2629,, 2021.

EGU21-172 | vPICO presentations | EMRP1.14

Investigation of failure and damage processes of andesitic rocks

Özge Dinç Göğüş, Deniz Yılmaz, Elif Avşar, and Kamil Kayabalı

In this research, failure and deformation processes of andesitic rocks are investigated through laboratory and discrete element modeling (DEM) analysis to reveal the transition of the cracking, namely from microscale to mesoscale (lab scale). For this purpose, the mechanical properties of Ankara andesites were initially investigated by performing uniaxial - triaxial compressive and indirect tensile laboratory tests. Further, these properties were used as reference parameters for the calibration process in a numerical model, generated through a three-dimensional open source code (Yade) based on the discrete element method (DEM). Our results show that during the linear-elastic region of the stress-strain curve, the major mechanism of rock behavior is driven by tensile cracks. When the crack damage threshold is reached, as a result of plastic strain, the strength related to the inherent cohesion significantly decreases and damage in the rock cannot be prevented anymore. At the peak stress of the curve, both tensile and shear cracks accumulate, intensively. Even the mesoscale failure mechanism is controlled by shearing at the residual stage of the yielding, based on the micro-scale process in the DEM model, the number of micro shear cracks is very limited compared to the tensile ones. This finding shows that the friction takes the control of the damage process as the only driving force during residual phase time.  

How to cite: Dinç Göğüş, Ö., Yılmaz, D., Avşar, E., and Kayabalı, K.: Investigation of failure and damage processes of andesitic rocks , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-172,, 2021.

EGU21-9048 | vPICO presentations | EMRP1.14

Seeing and hearing quasi-static shear band localization in a sandstone

Alexis Cartwright-Taylor, Ian G. Main, Ian B. Butler, Florian Fusseis, Maria-Daphne Mangriotis, Andrew Curtis, Andrew Bell, Martin Ling, Edward Andò, Roberto Rizzo, Sina Marti, Derek Leung, and Oxana Magdysyuk

The localisation of structural damage, in the form of faults and fractures, along a distinct and emergent fault plane is the key driving mechanism for catastrophic failure in the brittle Earth. However, due to the speed at which stable crack growth transitions to dynamic rupture, the precise mechanisms involved in localisation as a pathway to fault formation remain unknown. Understanding these mechanisms is critical to understanding and forecasting earthquakes, including induced seismicity, landslides and volcanic eruptions, as well as failure of man-made materials and structures. We used time-resolved synchrotron x-ray microtomography to image in-situ damage localisation at the micron scale and at bulk axial strain rates down to 10-7 s-1. By controlling the rate of micro-fracturing events during a triaxial deformation experiment, we deliberately slowed the strain localisation process from seconds to minutes as failure approached. This approach, originally established to indirectly image fault nucleation and propagation with acoustic emissions, is completely novel in synchrotron x-ray microtomography and has enabled us to image directly processes that are normally too transient even for fast synchrotron imaging methods. Here, we first present the experimental apparatus and control system used to acquire the data, followed by damage localisation and shear zone development in a sample of Clashach sandstone viewed in unprecedented detail. Time-resolved microtomography images demonstrate a strong intrinsic correlation between shear and dilatant strain in the localised zone, with bulk shear strain accomodated by the nucleation and rotation of en-echelon tensile microcracks within a grain-scale shear band. Rotation is accompanied by antithetic to synthetic shear sliding of neighbouring crack surfaces as they rotate. The evolving 4D strain field, measured with incremental digital volume correlation between pairs of recorded x-ray tomographic volumes, independently confirm the correlation between shear and dilatant strain and show how strain localises spontaneously, first through exploration of several competing shear bands at peak stress before transitioning to failure along the optimally-oriented final fault plane. In order to ‘ground-truth’ inferences made from bulk measurements and seismic waves (the primary method of detecting deformation at the field-scale where direct imaging of the subsurface is impossible), we (a) compare rupture energy estimates from local slip measurements with those from bulk slip data, and (b) use AE source location estimates to identify individual cracks and other local changes in the microstucture that may explain the AE source.

How to cite: Cartwright-Taylor, A., Main, I. G., Butler, I. B., Fusseis, F., Mangriotis, M.-D., Curtis, A., Bell, A., Ling, M., Andò, E., Rizzo, R., Marti, S., Leung, D., and Magdysyuk, O.: Seeing and hearing quasi-static shear band localization in a sandstone, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9048,, 2021.

EGU21-14332 | vPICO presentations | EMRP1.14

Deformational characteristics of thermally treated sandstone from an underground coalmine fire region, India

Adarsh Tripathi, Noopur Gupta, Ashok Kumar Singh, Nachiketa Rai, and Anindya Pain

The Jharia region of lower Gondwana in India is one of the largest Underground Coalmine Fires (UCF) affected coalfield in the world. The UCF induced small scale as well as large-scale surface fracturing often creates the life-threatening conditions to coal miners and local surroundings. So, there is a need to understand the thermomechanical behaviour of coal measures rocks to predict the land disturbances in thermo-environmental conditions. It will provide an insight into the UCF induced subsidence mechanism and its preventive measures. The Jharia coal field predominantly consists of sandstone (75-80% by volume) and rest is composed of coal, shale and carbonaceous shale. The present study focuses on thermo-mechanical behaviour of Barakar sandstone (BS) under elevated temperatures. The cores of BS sample were prepared according to the ISRM standards. Further, samples were grouped and thermally treated in temperature range of 25°C, 100°C, 150°C, 300°C, 400°C, 500°C, 600 °C, 700°C and 800°C at a heating rate of 5°C/min for 24 hours in furnace.  Then, these thermally treated BS samples were subjected to laboratory test for stress-strain characteristics. In the process of deformational characteristics evaluation, effect of mineralogical changes and mode of fracture pattern were also studied at the mentioned elevated temperature. Based on the obtained results, the deformational behaviour of thermally treated BS specimens can be grouped into three zones, viz., zone 1 (25-300°C), zone 2 (300-500°C) and zone 3 (500-800°C). In zone 1, the characteristics of the stress-strain curve is similar to those under air dried sandstone specimen. However, small increment in stiffness were observed upto 300°C. The stress-strain curves in this zone shows dominantly brittle fracturing. The increment in stiffness may be related to evaporation of pore water that increases the cohesion between the mineral grains resulting higher stiffness value. In zone 2, the deformation pattern again shows brittle fracturing with continuous decrement in stiffness. The reduction in stiffness may be related to thermally induced porosity and increased microcrack density. In zone 3, the stress strain curve is observed to be concave upward. It indicates the pseudo-ductile behaviour of the thermally treated BS specimens. The observed results suggest a typical behaviour of deformation pattern under UCF induced rock fracturing which may be useful in predicting the land subsidence in UCF affected areas. Present research outcome may be used to design the support measures to reduce the associated hazards.

How to cite: Tripathi, A., Gupta, N., Singh, A. K., Rai, N., and Pain, A.: Deformational characteristics of thermally treated sandstone from an underground coalmine fire region, India, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14332,, 2021.

EGU21-5959 | vPICO presentations | EMRP1.14

Effect of a heterogeneity on tensile failure: interaction between fractures in a limestone

Anne Pluymakers, Auke Barnhoorn, and Richard Bakker

Not all rocks are perfect. Frequently heterogeneities will be present, either in the form of pre-existing fractures, or in the form of sealed fractures. To date, investigation of sample heterogeneity, specifically tensile strength and strength anisotropy has focused on layered rocks, such as shales, sandstones and gneisses. Data is lacking on the effect of single planar heterogeneities, such as pre-existing fractures or stylolites, even though these frequently occur in geo-energy settings.

We have performed Brazilian Disc tests on limestone samples containing planar heterogeneities, investigating Brazilian test Strength (BtS) and the effects of orientation on strength. We used prefractured Indiana limestone to represent a planar heterogeneity without cohesion and Treuchtlinger Marmor samples with central stylolites to represent a planar heterogeneity of unknown strength (as an example of a sealed fracture). The planar discontinuity was set at different rotation angles of approximately 0–20–30–45–60–90⁰, where 90⁰ (steep angle) is parallel to the principal loading direction, and 0⁰ (low angle) to the horizontal axis of the sample. All experiments were filmed, and where possible Particle Image Velocimetry was used to determine internal particle motion. Moreover, we used a 2D Comsol model in which we simplified the stylolite surface as a sinusoid. The model was used to qualitatively determine how i) a different period of the sinusoid and ii) relative strength of sinusoid/matrix affect the results.

Our results show that all imperfect samples are weaker than intact samples. The 2D Comsol model indicates that the qualitative results remain unaffected by changing the period (assumed to be representative of roughness) of the cohesive heterogeneity, nor by the relative strength contrast: the location of the first fracture remains unaffected. For both heterogeneity types, the fracture patterns can be divided into four categories, with two clear endmembers, and a more diffusive subdivision in between.

For a cohesion-less heterogeneity:

  • steep angles lead to frictional sliding along the interface, and only a small hypothesized permeability increase.
  • Intermediate angles lead to a combination of tensile failure of the matrix and sliding along the interface, where for steeper angles more new fractures form which follow the path of the existing fracture.
  • Low angles lead to closure of the old fracture and new tensile failure.

For a cohesive heterogeneity of unknown cohesion:

  • Steep angles lead to intensive failure of the heterogeneous zone, attributed to the presence of a stress concentrator.
  • Intermediate angles lead to partial failure along the heterogeneous zone, and the formation of new fractures in the matrix, potentially instigated by mode II failure to accommodate motion.
  • Low angles lead to the formation of a new fracture plus opening within the heterogeneous zone.

These results imply that hydrofracture (i.e. creating tensile stresses) of a stylolite-rich zone will lead to more fractures than fractures in a homogeneous zone, where the orientation of the stylolites and bedding will control the orientation of the permeable pathways.

How to cite: Pluymakers, A., Barnhoorn, A., and Bakker, R.: Effect of a heterogeneity on tensile failure: interaction between fractures in a limestone, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5959,, 2021.

EGU21-11553 | vPICO presentations | EMRP1.14

A postmortem approach via X-ray computed tomography and thresholding to investigate fracture network evolution in Onagawa shale

Eranga Jayawickrama, Jun Muto, Osamu Sasaki, and Hiroyuki Nagahama

A postmortem technique is introduced to investigate the fracture connectivity evolution under elevated confining pressures via a sensitivity analysis. Three Onagawa shale samples are deformed under brittle, ductile, and transition conditions, by increasing the confining pressures. Brittle deformation is characterized by longitudinal splitting of the sample at 3% axial strain, and the onset of transition from brittle to ductile deformation is between 4% ~ 5% axial strain. The ductile deformation is characterized by a distributed conjugate fracture network and strain hardening. In completion of the deformation, the samples are scanned in a commercially available X-ray CT machine. The grayscale values of the primary 2D images were reversed, stacked, and surface rendered to obtain the 3D volume distribution of the fractures. Reversing and surface rendering allowed the acquisition of volume and surface data of the fractures along with their direct visualization. Further, utilizing a residual analysis, the voxel value density distribution that fabricated the fracture network is extracted (Residual histogram). Thresholding of the residual histogram generated volume segments of the final fracture network demonstrating the sensitivity of the fracture network to the choice of threshold. Voxel volumes of fractures alone are obtained by thresholding post-peak voxel values of the residual histogram and consecutive post-peak thresholding shows that the generated volume segments of the fracture network can be utilized to interpret, possible nucleation sites after strain localization, propagation of fractures, and coalescence. Fracture connectivity is quantified by means of relative entropy from information theory, and the relative entropy of size distribution of fracture volumes showed that it is closer to zero with the fractures being well connected. Moreover, the cumulative fracture volume shows a power-law growth towards the failure after a unique threshold to each sample. These results have been validated by previous acoustic emission studies and a 4D tomographic investigation on strain localization of shale. Therefore, despite the postmortem nature of the investigation, the new technique has opened possibilities to investigate the fracture properties and their evolution under elevated confining pressures.

How to cite: Jayawickrama, E., Muto, J., Sasaki, O., and Nagahama, H.: A postmortem approach via X-ray computed tomography and thresholding to investigate fracture network evolution in Onagawa shale, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-11553,, 2021.

Large rock slopes instabilities form over long timescales through progressive rock mass strength weakening of initially stable slopes. Progressive rock mass damage is driven by environmental loads and is thus strongly dependent on the local setting and environmental conditions of the rock slope, which can vary over time. It is often assumed that the strong variations of the thermal and hydraulic boundary conditions during deglaciation in combination with unloading due to ice downwasting cause enhanced rates of rock mass damage. However, in-situ observations to quantify deformation, damage and the relevance of different drivers in such environments are rare. This presentation is related to the contribution of Oestreicher et al., presenting in the same session, addressing similar questions, but at different scales and based on different field data and analysis.

In this contribution we analyze continuous pore pressure, temperature and micrometer-scale deformation time series from a subsurface monitoring system comprised of three, 50 m deep, highly instrumented boreholes in a crystalline rock slope which is located beside the rapidly retreating glacier tongue of the Great Aletsch Glacier (Switzerland). We compare high-resolution reversible and irreversible deformation signals with potential drivers, including locally measured pore pressure fluctuation, rock temperature variations, and nearby earthquakes. We show that shallow (10 - 15 m deep) deformations in our rock slope are dominated by thermo-mechanical forcing, whereas deformation measured below this depth is mainly driven by hydro-mechanical effects related to pore pressure fluctuations. Both reversible deformation and irreversible damage events occur more frequently during the snow-free summer season, when we observe higher dynamics in thermal and hydraulic boundary conditions. In our 2.5 years long time series, we do not find any significant deformation event coinciding with a nearby earthquake. Additionally, we discuss differences in the deformation signal with respect to the stability state and the rock mass quality at the different monitoring locations. Also, we assess longer term impacts of glacier retreat and ice downwasting on rock slope deformation and damage. Such information is critical for an improved understanding and quantification of factors contributing to the formation of paraglacial rock slope instabilities.

How to cite: Hugentobler, M., Aaron, J., and Loew, S.: Drivers for micrometer-scale deformation and progressive rock mass damage in a deglaciating rock slope – insights from subsurface in-situ monitoring, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5192,, 2021.

EGU21-8872 | vPICO presentations | EMRP1.14

Multiscale fracture density analysis at Stromboli Volcano, Italy: implications to flank stability

Thomas Alcock, Sergio Vinciguerra, Phillip Benson, and Federico Vagnon

Stromboli volcano has experienced four sector collapses over the past 13 thousand years, resulting in the formation of the Sciara del Fuoco (SDF) horseshoe-shaped depression and an inferred NE / SW striking rift zone across the SDF and the western sector of the island. These events have resulted in the formation of steep depressions on the slopes on the volcano where episodes of instability are continuously being observed and recorded. This study aims to quantify the fracture density inside and outside the rift zone to identify potential damaged zones that could reduce the edifice strength and promote fracturing. In order to do so we have carried out a multiscale analysis, by integrating satellite observations, field work and seismic and electrical resistivity analyses on cm scales blocks belonging to 11 lava units from the main volcanic cycles that have built the volcano edifice, ie. Paleostromboli, Nestromboli and Vancori. 0.5 m resolution Pleiades satellite data has been first used to highlight 23635 distinct linear features across the island. Fracture density has been calculated using Fracpaq based on the Mauldon et al (2001) method to determine the average fracture density of a given area on the basis of the average length of drawn segments within a predetermined circular area. 41.8 % of total fracture density is found around intrusions and fissures, with the summit area and the slopes of SDF having the highest average fracture density of 5.279  . Density, porosity, P- wave velocity in dry and wet conditions and electrical resistivity (in wet conditions) were measured  via an ultrasonic pulse generator and acquisition system (Pundit) and an on purpose built measuring quadrupole on cm scale blocks of lavas collected from both within and outside the proposed rift zone to assess the physical state and the crack damage of the different lava units.  Preliminary results show that P-wave velocity between ~ 2.25 km/s < Vp < 5km/s decreases with porosity while there is high variability electrical resistivity with 21.7 < ρ < 590 Ohm * m. This is presumably due to the lavas texture and the variable content of bubble/vesicles porosity and crack damage, that is reflected by an effective overall porosity between 0 and 9 %. Higher porosity is generally mirrored by lower p-wave velocity values. Neostromboli blocks show the most variability in both P-wave velocity and electrical resistivity. Further work will assess crack density throughout optical analyses and systematically investigate the UCS and elastic moduli. This integrated approach is expected to provide a multiscale fracture density and allow to develop further laboratory testing on how slip surfaces can evolve to a flank collapse at Stromboli.

How to cite: Alcock, T., Vinciguerra, S., Benson, P., and Vagnon, F.: Multiscale fracture density analysis at Stromboli Volcano, Italy: implications to flank stability, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8872,, 2021.

The Fracture Induced Electromagnetic Radiation (FEMR) technique has gradually progressed in the past decade as a useful geophysical tool to determine the direction and magnitude of recent crustal stresses, visualize the modification and realignment of stresses inside tunnels thus proving to be an important precursor for geohazards, earthquake forecasting, as well as delineate landslide-prone slip planes in unstable regions. Its working principle is based on the generation of geogenic electromagnetic radiation emanating from the brittle rock bodies that are fractured being subjected to an incremental increase of the differential stress in the near-surface of the Earth’s crust. The “Process zone” at the fractured crack tip contains numerous microcracks which subsequently creates dipoles due to the polarization of charges on such microcrack tips which rapidly oscillates emitting FEMR waves of frequencies between KHz to MHz range. The coalescence of the microcracks eventually leads to a macro failure dampening the amplitude of the FEMR pulses. The attenuation of FEMR pulses is comparatively lesser than seismic waves making it a more efficient precursor to potential tectonic activities indicating an upcoming earthquake a few hours/days before the actual event. In the current study, we have attempted to exploit this technique to identify the locations of the potential active faults across the tectonically active Narmada-Son Lineament (NSL), Central India. Although the first tectonic stage involved rifting and formation of the NSL during the Precambrian time, the rifting continued at least till the time of Gondwana deposition. Later, tectonic inversion took place as a result of the collision between the Indian and the Eurasian plate resulting in reverse reactivation of the faults. Episodic reverse movement along NSL caused recurrent earthquakes and linear disposition of the sediments that were deposited at the foothills of the Satpura Horst. Although the origin of East-West trending NSL dates back to the Precambrian time, it is very much tectonically active as manifested by recent earthquakes. The study has been conducted by taking linear FEMR readings across 3 traverses along the NSL which on analysis provides an idea about the potential active faults, their locations, and frequency of occurrence. The accumulation of strain in the brittle rocks that can eventually lead to a macro failure is demarcated as an anomalous increase in the amplitude of the FEMR pulses indicative of an upcoming tectonic episode in the region. To further corroborate the analysis, we have attempted to determine the neo-tectonic activity in the region by calculating the morphometric parameters across the Khandwa-Itarsi-Jabalpur region, Central India. Finally, we attempt to comment on the tectonic evolution of Central India in the recent past. We also encourage researchers to adapt the novel technique of FEMR which is swift, affordable, and feasible compared to conventional techniques deployed to survey the active tectonics of a region.

How to cite: Das, S. and Mallik, J.: Delineating active faults across the Narmada-Son Lineament (NSL), India using the technique of Fracture Induced Electromagnetic Radiation (FEMR), EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-2370,, 2021.

EGU21-9403 | vPICO presentations | EMRP1.14

Integrated geophysical and geotechnical monitoring for multiscale rock mass damaging investigation at the Acuto Field-Lab (Italy)

Guglielmo Grechi, Danilo D'Angiò, Matteo Fiorucci, Roberto Iannucci, Luca Lenti, Gian Marco Marmoni, and Salvatore Martino

Rock mass damaging has become a topic of great interest in the engineering-geology research community during the last decades as it can significantly influence slopes stability. In this sense, the study of mechanics and dynamics of jointed rock masses represents a challenge because it will allow to better understand how external continuous and transient stressors can influence the short- to long-term stability controlling their pre-failure behavior. Consequently, the detection of permanent changes in physical and mechanical parameters, due to periodic or transient stressors, is an important target to mitigate the related geological risk as it can potentially lead rock masses to failure, especially when infrastructures and natural or cultural heritages are exposed elements. In this framework, the Acuto field laboratory (Central Italy) has been designed and implemented in 2016 within an abandoned quarry by employing an integrated geotechnical and geophysical monitoring system, with the aim of investigating how natural and anthropic conditioning factors could lead fractured rock masses to failure. The integrated monitoring system, which is installed on a potentially unstable 20-m3 jointed rock block, is composed of several strain devices (i.e., strain gauges -SG- and jointmeters -JM-), one fully equipped weather station, one rock thermometer, eight high-sensitivity microseismic uniaxial accelerometers and optical and InfraRed Thermal cameras. The acquisition of long-term monitoring time-series, coupling multimethodological approaches, allowed to establish cause-to-effect relationships among different environmental stressors and induced strain effects, highlighting the continuous action of thermal stresses on rock mass deformations both at the daily and seasonal timescales. In fact, while the analysis of thermal and strain monitoring data allowed to characterize the cyclic contraction and relaxation response of major rock fractures and microcracks to temperature fluctuations, the microseismic monitoring array was able to detect during thermal transient (i.e., freezing conditions) the occurrence of microseismic emissions potentially related to the genesis or progressive growth of pre-existing cracks.

Starting from 2018, experimental activities at the Acuto field lab are supported by the “Dipartimento di Eccellenza” project of the Italian Ministry of Education Universities and Research funds attributed to the Department of Earth Sciences of the University of Rome “Sapienza”.  In this framework, the Acuto filed laboratory will undergo a structural upgrade that will be aimed at the investigation of two new sectors of the abandoned quarry. These new sectors will be instrumented with innovative thermal profiles probe, fiber Brag grating sensors and traditional SG and JM for detailed stress-strain monitoring, acoustic emission sensors and high-frequency and low-frequency geophones for ambient seismic noise monitoring and microseismic events detection as well as accelerometers for evaluating the rock mass response in the case of seismic shaking. The main goal of such an improvement will be both technical and methodological, and will shed light on the application of integrated geophysical and geotechnical monitoring approaches in investigating the multiscale rock mass damaging process as well as the detection of rock mass failure precursors by using non-conventional combinations and configurations of geotechnical and broad-band geophysical devices.

How to cite: Grechi, G., D'Angiò, D., Fiorucci, M., Iannucci, R., Lenti, L., Marmoni, G. M., and Martino, S.: Integrated geophysical and geotechnical monitoring for multiscale rock mass damaging investigation at the Acuto Field-Lab (Italy), EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9403,, 2021.

EGU21-8082 | vPICO presentations | EMRP1.14

Drivers of reversible and irreversible slope deformations in a paraglacial environment

Nicolas Oestreicher, Clément Roques, Marc Hugentobler, Jordan Aaron, and Simon Loew

Retreating glaciers around the world lead to rapid and profound changes in the surrounding landscapes. In the Alps, many glaciers are rapidly retreating and downwasting, substantially modifying stresses and hydro-thermal boundary conditions on the adjacent slopes. There is an increase in observations of bedrock responses and the formation of large-scale instabilities in paraglacial environments, but still a little knowledge about the underlying preparatory factors and drivers. This presentation is linked to the one from Hugentobler et al. in the same session. Both studies take place in the same catchment and address the same questions at different spatial scales, with other techniques and datasets.

We analyse surface deformation data monitored in a crystalline bedrock catchment, on the recently deglaciated slopes of the Great Aletsch Glacier (Valais, Switzerland). Our monitoring system has been in operation for six years and comprises 93 reflectors, 2 robotic TPS, and 4 cGPS stations distributed on both sides of the glacier tongue. This unique dataset allows studying the main processes involved at relevant spatial and temporal scales. The response of potential drivers for reversible and irreversible deformation is evaluated through combined multivariate (vbICA) and cross-correlation statistical analysis. We found that the variability in deformation near the glacier tongue is primarily controlled by glacier unloading through melting and seasonal groundwater fluctuations. At the catchment scale, the later effect is poroelastic and hence reversible, but we argue that it could also induce hydromechanical fatigue. By investigating the deformation's spatial pattern, we observed that the reversible deformation is mostly controlled by discrete structures such as hectometer-scale brittle-ductile shear zones striking subparallel to the valley axes and the main Alpine foliation. Field mapping and pressure monitoring during borehole drilling suggest that infiltration into the fractured rockmass is very heterogeneous and mainly controlled by the presence of interconnected tensile fractures.

How to cite: Oestreicher, N., Roques, C., Hugentobler, M., Aaron, J., and Loew, S.: Drivers of reversible and irreversible slope deformations in a paraglacial environment, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8082,, 2021.

EGU21-5083 | vPICO presentations | EMRP1.14

The Luco dei Marsi deep-seated gravitational deformation: first evidence of a basal shear zone in the central Apennine mountain belt (Italy) 

Emiliano Di Luzio, Marco Emanuele Discenza, Maria Luisa Putignano, Mariacarmela Minnillo, Diego Di Martire, and Carlo Esposito

The nature of the boundary between deforming rock masses and stable bedrock is a significant issue in the scientific debate on Deep-Seated Gravitational Slope Deformations (DSGSDs). In many DSGSDs the deforming masses move on a continuous sliding surface or thick basal shear zone (BSZ) [1-3]. This last feature is due to viscous and plastic deformations and was observed (or inferred) in many worldwide sites [4]. However, no clear evidence has been documented in the geological context of the Apennine belt, despite the several cases of DSGSDs documented in this region [5-6].

This work describes a peculiar case of a BSZ found in the central part of the Apennine belt and observed at the bottom of a DSGSD which affects the Meso-Cenozoic carbonate ridge overhanging the Luco dei Marsi village (Abruzzi region). The NNW-SSE oriented mountain range is a thrust-related Miocene anticline, edged on the east by an intramountain tectonic depression originated by Plio-Quaternary normal faulting. The BSZ appears on the field as a several meters-thick cataclastic breccia with fine matrix developed into Upper Cretaceous, biodetritic limestone and featuring diffuse rock damage.

The gravity-driven process was investigated through field survey, aerial photo interpretation and remote sensing (SAR interferometry) and framed into a geological model which was reconstructed also basing on geophysical evidence from the CROP 11 deep seismic profile. The effects on slope deformation determined by progressive displacements along normal faults and consequent unconfinement at the toe of the slope was analysed by a multiple-step numerical modelling constrained to physical and mechanical properties of rock mass.

The model results outline the tectonic control on DSGSD development at the anticline axial zone and confirm the gravitational origin of the rock mass damage within the BSZ. Gravity-driven deformations were coexistent with Quaternary tectonic processes and the westward (backward) migration of normal faulting from the basin margin to the inner zone of the deforming slope.


[1] Agliardi F., Crosta G.B., Zanchi A., (2001). Structural constraints on deep-seated slope deformation kinematics. Engineering Geology 59(1-2), 83-102.

[2] Madritsch H., Millen B.M.J., (2007). Hydrogeologic evidence for a continuous basal shear zone within a deep-seated gravitational slope deformation (Eastern Alps, Tyrol, Austria). Landslides 4(2), 149-162.

[3] Zangerl C., Eberhardt E., Perzlmaier S., (2010). Kinematic behavior and velocity characteristics of a complex deep-seated crystalline rockslide system in relation to its interaction with a dam reservoir. Engineering Geology 112(1-4), 53-67.

[4] Crosta G.B., Frattini P., Agliardi F., (2013). Deep seated gravitational slope deformations in the European Alps. Tectonophysics 605, 13-33.

[5] Discenza M.E., Esposito C., Martino S., Petitta M., Prestininzi A., Scarascia-Mugnozza G., (2011). The gravitational slope deformation of Mt. Rocchetta ridge (central Apennines, Italy): Geological-evolutionary model and numerical analysis. Bulletin of Engineering Geology and the Environment,70(4), 559-575.

[6] Esposito C., Di Luzio E., Scarascia-Mugnozza G., Bianchi Fasani G., (2014). Mutual interactions between slope-scale gravitational processes and morpho-structural evolution of central Apennines (Italy): review of some selected case histories. Rendiconti Lincei. Scienze Fisiche e Naturali 25, 161-155.

How to cite: Di Luzio, E., Discenza, M. E., Putignano, M. L., Minnillo, M., Di Martire, D., and Esposito, C.: The Luco dei Marsi deep-seated gravitational deformation: first evidence of a basal shear zone in the central Apennine mountain belt (Italy) , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5083,, 2021.

Rainfall-induced landslides pose a substantial risk to people and infrastructure worldwide, but their mechanical behavior is not well understood. As a result, hazard predictions for these landslides, especially for rainfall and slope-failure correlations, remain an active area of research. Many operational rainfall-induced landslide hazard maps still assume a classical Coulomb type failure criterion where slope-failure must occur either before or at peak subsurface pore pressure reached during a precipitation event. Using satellite-derived surface precipitation data and soil infiltration simulations over a 15 day period preceding 121 rainfall-induced landslides across India, we find that these events occurred systematically 2-12 days after the simulated peak pore pressures on the inferred failure slope nucleating between 0.5 and 5 m depth. These observations cannot be explained with the Coulomb failure criterion, since failure on these slopes is significantly delayed behind the occurrence of the inferred strength minimum. Instead, in this study, we investigate whether a slope failure model with time- and slip-variable shear strength, governed by the rate-state friction (RSF) equations widely used in earthquake mechanics, can explain the observed ranges of time-delays between slope failure and inferred peak pore pressure.

To concentrate on the role of the constitutive behavior of the failure surface, we examine spring-slider dynamics under a classical RSF framework driven by variable on-slope and far-field pore pressure and flux time histories. We derive analytical expressions for the time-to-failure of such a spring-slider under simple pore pressure perturbation histories and find that the delay-times can vary significantly depending on the laboratory derivable RSF parameters, soil bulk properties, and particulars of the pressure history. We further examine the roles of dilatancy strengthening and pore compaction in determining the time-lag between peak pore pressure and slope failure. We find that dilatancy can have either a stabilizing or a destabilizing effect on slope failure depending on the hydrological and mechanical properties of the failure plane and the soil column. Finally, we show with numerical simulations that periodic pore pressure or flux oscillations can also drive asynchronous repeated slope failures in both the presence and absence of the coupling of pore pressure and shear deformation. Our results show that the observed rainfall-landslide correlations for these 121 landslides can be explained with inherently time- and slip-dependent shear strength prescriptions like the RSF equations. This, in turn, implies that realistic landslide hazard monitoring might require the examination of soil shear strength under the experimental protocols widely used in rock friction experiments to determine whether the constant friction assumption inherent in the Coulomb criterion needs to be revised in favor of RSF or similar constitutive equations for shallow landslides.

How to cite: Paul, K., Bhattacharya, P., and Misra, S.: Investigation of The Observed Time-lag Between the Simulated Peak Pore Pressures and Slope Failures in Rainfall Induced Landslides: A Numerical Approach, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15352,, 2021.

In 2009, a large-scale landslide was triggered by typhoon rainfall and buried an entire village, which named Hsiaolin and located in Taiwan. 
After that, Soil and Water Conservation Bureau (SWCB) has promoted a national project for the prevention work of large-scale landslide. The national project includes with the investigation of potential area, the design of monitoring system, and the design of warning system, etc. 
The investigation of potential large-sclae landslide was based on the  digital elevation model with 1 meter resolution. However, the investigation of the underground was lack and not clear enough. Therefore, the specific landslide's body is hardly to estimate and it causes difficulty in follow-up works. 
This study applied two methods to investigate the scenario of slope failure. The first method is based on the limited equilibrium method, which proposed by Yoshino and Uchida (2019). The method was used to search the specific region of unstable slope based on a series of high-resolution digital elevation models. After the specific region of unstable slope was confirmed, the landslide can be simulated by a numerical model, which this study proposed to represent the entire landslide process from occurrence to post-failure . 
These proposed methods were applied at Baolai area, south Taiwan to track the evolution of the potential area. The failure scenario could be evaluated by the proposed numerical model. By this study, the investigation of underground can be evaluated and these results are very important information for the design of monitoring system.

How to cite: Tsai, Y. and Lee, W.: Application of the high-resolution digital elevation model in an integrated numerical method for the occurrence mechanism and post-failure behavior of the landslide, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-10764,, 2021.

In the pre-feasibility study stage, only a small amount of borehole data can be obtained. Since the available geological information is insufficient, the engineering geological conditions of the project can only be preliminarily and approximately estimated during this stage. In this study, we attempt to seek a method to make a preliminary analysis and evaluation of the stability of the surrounding rock masses of an underground rock carven project, which makes full and optimum use of the limited borehole data to accomplish the assessment of the investigated site. The basic information on rock fractures is extracted from the borehole Television logging data and the fracture extension directions are also determined. Providing that the cracks detected in the borehole would extend to the cavern area, the cracks with appropriate direction, larger width and larger hydraulic conductivity can be selected. These selected cracks are considered in the numerical model established using ANSYS, and the stability of surrounding rock of cavern is analyzed under this situation. In the absence of large amount of borehole data, this method, which set up an extreme case, can be used to analyze possible failure of rock mass under extreme adverse conduction in advance. In general, the proposed method for stability analysis could contribute to the design and construction practice of a tunnel project constructed in fractured rock masses.

How to cite: Tan, L. and Xiang, W.: Stability analysis of the surrounding rock of cavern under extreme situation based on limited borehole data., EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-230,, 2021.

EGU21-9667 | vPICO presentations | EMRP1.14

Cosmic-ray muography applications in underground tunnelling

Pasi Kuusiniemi, Marko Holma, and Zongxian Zhang

The novel geophysical remote imaging method of muography is based on cosmic-ray induced muon particles that are detected after passing through the media of interest. If the studied objects are solid, their sizes can vary from meters to up to kilometres. In terms of penetration capability, muography can be placed between methods based on X-rays and those using seismic waves. The most famous objects imaged with muography are pyramids (e.g., Khufu's Pyramid at Giza in Egypt) and volcanoes (e.g., Mt Etna in Italy). One clear advantage of muography compared with seismic methods is that muons, unlike seismic waves, do not reflect from geological interfaces. In addition, the scattering phenomenon is a minor issue and needs consideration only at low-energy muons. Raw data must be corrected according to topography. On the basis of extensive numeric simulations of Hivert et al. (2017), the lowest density variations observable for muography with a significant level of 3σ (a typical significance level in physics) are around 2% at 150 m, 4% at 300 m, and 10% at 700 m of depth, respectively. If these numbers are extrapolated to depths below 100 m, the mean density differences in the range of 1% are likely within the observation capability of muography. It is also worth to note that the 1% difference in a mean rock density results in an approximately 3% difference in the muon flux. This indicates that muon flux measurements are very sensitive to the density variations of rocks.

In underground tunnelling, muography has at least four applications: (1) muography can be used to detect a potential risk (such as a water reservoir, a weak zone with loose rocks, boulders, etc.) before or during tunnelling, (2) muography can be employed to monitor overburden rock behaviour during tunnelling operation to avoid risks like the roof cave-ins, (3) muography can be applied to monitor the overburdening rock masses in tunnels after they are excavated to predict and avoid the collapse of rock mass, and (4) muography can be used to estimate the size and volume of a rock mass collapse in a tunnel since the volume of the collapsed rocks must have markedly smaller density than original overburden rock mass. In an excavating tunnel project using a tunnel boring machine (TBM), a muon detector can be installed in the TBM during tunnelling. If there occurs a tunnel cave-in, muography can be employed in undamaged tunnels nearby (sideways or below) the collapse. If possible, the collapse can also be approached safely via an undamaged part of the collapsed tunnel. If none of these are available, borehole muography can be applied as a substitute solution. Whereas an undamaged underground tunnel is either filled by air or water, a collapsed tunnel segment is characterized by air and rock, or water and rock. In either case, the average density of the tunnel segment is increased. We are currently planning simulations and real-world tests to validate these assumptions.

How to cite: Kuusiniemi, P., Holma, M., and Zhang, Z.: Cosmic-ray muography applications in underground tunnelling, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9667,, 2021.

EMRP1.15 – Petrophysics and rock-physics across the scales: integrating models, laboratory experiments and field geophysical studies

Experimental rock physics allows the study of specific geological phenomena in a controlled manner. The experimental data are used to develop and calibrate predictive numerical models, which ultimately improve our understanding of natural processes and interpretation of field scale datasets. However, upscaling laboratory geophysical datasets to explain large scale geological complexes is challenging and by definition imprecise, as core-scale experiments are not fully representative of the events occurring in the field. This challenge gains in complexity with the increasing number of involved parameters and with changing environmental conditions. Nowadays, one of the most challenging rock physics areas is Carbon Capture Utilization and Storage (CCUS).

CCUS is a realistic global scale mitigation solution to tackle the excess of CO2 expelled from industrial production and sequestering it into deep reservoir formations, while attempting at generate energetic benefits from this process (e.g., enhanced oil recovery and geothermal energy). When CO2 is injected in rocks, the parental pore fluid is displaced by CO2 and therefore the effective bulk modulus of the fluid (and the rock) changes; if the rock and/or the fluid are reactive to CO2, then the rock frame and fluid composition are also changing parameters, affecting both the elastic and hydrodynamic properties of the original rock.

In this contribution, I examine the complexity of the geophysical interpretation of a specific CCUS-related process: CO2-induced salt precipitation in saline aquifers. In essence, injected CO2 dry out the brine and salt precipitates in the pores, leading to some degree of clogging that depends on, but not exclusively, the pore size, salt volume and pore connectivity. This phenomenon affects the injectivity of the CO2 storage site, which could lead to pressure built-up events and ultimately compromise the geomechanical integrity of the reservoir. Therefore, an early warning of CO2-induced salt precipitation is essential to apply specific mitigation strategies in a timely manner.

Salt precipitation is a rapid, self-enhanced process, with a footprint in the geophysical record that has to be isolated from that of the CO2 fluid replacement and other potential CO2-induced rock reactions. When the geophysical data include seismic and electromagnetic surveys, S-waves provide information about the changes in the rock frame, while P-waves and electrical resistivity would help to distinguish between CO2, brine and salt fractions. Here, I use a dataset generated in the laboratory using ultrasonic P- and S-waves attributes, resistivity and axial deformation, during a CO2-brine flow-through experiment with salt precipitation, to discuss the potential of the rock physics to detect and quantify dynamic changes in the bulk properties of reservoir rocks.  

How to cite: Falcon-Suarez, I. H.: Joint elastic-electrical monitoring for detecting and quantifying CO2-induced salt precipitation during geological carbon and storage, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9920,, 2021.

EGU21-9578 | vPICO presentations | EMRP1.15

Pore-scale hydrodynamic evolution within carbonate rock during CO2 injection and sequestration

Chi Zhang, Siyan Liu, and Reza Barati

The continuously rising threat of global warming caused by human activities related to CO2 emission is facilitating the development of greenhouse gas control technologies. Subsurface CO2 injection and sequestration is one of the promising techniques to store CO2 in the subsurface.  However, during CO2 injection, the mechanisms of processes like injectant immobilizations and trapping and pore-scale geochemical reactions such as mineral dissolution/precipitation are not well understood. Consequently, the multi-physics modeling approach is essential to elucidate the impact of all potential factors during CO2 injection, thus to facilitate the optimization of this engineered application. 

Here, we propose a coupled framework to fully utilize the capabilities of the geochemical reaction solver PHREEQC while preserving the Lattice-Boltzmann Method (LBM) high-resolution pore-scale fluid flow integrated with diffusion processes. The model can simulate the dynamic fluid-solid interactions with equilibrium, kinetics, and surface reactions under the reactive-transport scheme.  In a simplified 2D spherical pack, we focused on examining the impact of pore sizes, grain size distributions, porosity, and permeability on the calcite dissolution/precipitation rate. Our simulation results show that the higher permeability, injection rate, and more local pore connectivity would significantly increase the reaction rate, then accelerate the pore-scale geometrical evolutions. Meanwhile, model accuracy is not sacrificed by reducing the number of reactants/species within the system.

Our modeling framework provides high-resolution details of the pore-scale fluid-solid interaction dynamics. To gain more insights into the mineral-fluid interfacial properties during CO2 sequestration, our next step is to combine the electrodynamic forces into the model. Potentially, the proposed framework can be used for model upscaling and adaptive subsurface management in the future.  

How to cite: Zhang, C., Liu, S., and Barati, R.: Pore-scale hydrodynamic evolution within carbonate rock during CO2 injection and sequestration, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9578,, 2021.

EGU21-13987 | vPICO presentations | EMRP1.15

The Effects of Surface Roughness on Fluid Displacement Mechanisms and Residual Residual Trapping - A Pore Scale Investigation

Rumbidzai Nhunduru, Omid Shahrokhi, Krystian Wlodarczyk, Amir Jahanbakhsh, Susana Garcia, and Mercedes Maroto-Valer

Immiscible fluid displacement and the trapping of residual oil and gas phases in the pore spaces of reservoir rocks is critical to geological operations such as carbon geo-sequestration and enhanced oil recovery. In carbon geo-sequestration, residual trapping is advantageous because it ensures long-term storage security of carbon dioxide (CO2). In contrast, residual trapping can pose significant challenges during waterflooding in oil recovery operations where large volumes of oil may remain trapped in the interstitial spaces of the porous reservoir rock and cannot be extracted, thereby reducing the efficiency of the recovery process. In such operations, residual trapping is strongly influenced by the inherent surface roughness of the solid rock matrix amongst many factors. Surface roughness occurs in natural reservoir rocks as a result of geological processes that physically, chemically or biologically convert sediments into sedimentary rock (known as diagenesis) and weathering.

The effects of surface roughness on immiscible two-phase flow are currently not well understood. Previous investigations into residual trapping in porous media have mainly focused on the influence of factors such as pore geometry, wettability, fluids interfacial tension, mobility ratio and injection scenarios. Although some of these studies acknowledge the potential effect of surface roughness, there is still a lack of quantitative characterization and understanding of the influence of surface roughness on immiscible two-phase displacements in porous media.

In this study, the impacts of surface roughness on immiscible two-phase displacement are quantified. Immiscible two-phase displacement of air by water was conducted in a custom laser-manufactured glass microfluidic chip (micromodel). The glass chip comprised a 2.5D micro-structure analogous to the pore network pattern (micro-structure) of a natural reservoir rock, Oolitic limestone. The pore network pattern consisted of cylindrical pillars 400 µm in diameter arranged in a rhombohedra type of packing, generated on to a glass substrate using an ultrafast, pulsed picosecond laser. Surface roughness is an innate characteristic of laser machined surfaces and as a result, small variations in depth of the porous micro-structure were observed (50 ± 8 µm). The average surface roughness (Sa) of the laser-machined structure was measured to be 1.2 μm.

Experimental results for the rough micromodel exhibit high repeatability of fluid displacement patterns (preferential flow pathways) demonstrating that surface roughness has a strong influence on fluid invasion patterns and sweep efficiency and its effects must not be ignored. To ascertain the effects of surface roughness on the fluid displacement process, a direct numerical simulation (DNS) of the fluid displacement process was performed in OpenFoam using the Volume of Fluid (VOF) method assuming zero surface roughness. Comparing the experimental results with the numerical simulations, we show that surface roughness can significantly enhance residual trapping in porous media by up to 49.2%.


How to cite: Nhunduru, R., Shahrokhi, O., Wlodarczyk, K., Jahanbakhsh, A., Garcia, S., and Maroto-Valer, M.: The Effects of Surface Roughness on Fluid Displacement Mechanisms and Residual Residual Trapping - A Pore Scale Investigation, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-13987,, 2021.

EGU21-16072 | vPICO presentations | EMRP1.15

Seismic wave attenuation due to fluid pressure diffusion at the mesoscopic scale: an experimental and numerical study

Samuel Chapman, Jan V. M. Borgomano, Beatriz Quintal, Sally M. Benson, and Jerome Fortin

Monitoring of the subsurface with seismic methods can be improved by better understanding the attenuation of seismic waves due to fluid pressure diffusion (FPD). In porous rocks saturated with multiple fluid phases the attenuation of seismic waves by FPD is sensitive to the mesoscopic scale distribution of the respective fluids. The relationship between fluid distribution and seismic wave attenuation could be used, for example, to assess the effectiveness of residual trapping of carbon dioxide (CO2) in the subsurface. Determining such relationships requires validating models of FPD with accurate laboratory measurements of seismic wave attenuation and modulus dispersion over a broad frequency range, and, in addition, characterising the fluid distribution during experiments. To address this challenge, experiments were performed on a Berea sandstone sample in which the exsolution of CO2 from water in the pore space of the sample was induced by a reduction in pore pressure. The fluid distribution was determined with X-ray computed tomography (CT) in a first set of experiments. The CO2 exosolved predominantly near the outlet, resulting in a heterogeneous fluid distribution along the sample length. In a second set of experiments, at similar pressure and temperature conditions, the forced oscillation method was used to measure the attenuation and modulus dispersion in the partially saturated sample over a broad frequency range (0.1 - 1000 Hz). Significant P-wave attenuation and dispersion was observed, while S-wave attenuation and dispersion were negligible. These observations suggest that the dominant mechanism of attenuation and dispersion was FPD. The attenuation and dispersion by FPD was subsequently modelled by solving Biot’s quasi-static equations of poroelasticity with the finite element method. The fluid saturation distribution determined from the X-ray CT was used in combination with a Reuss average to define a single phase effective fluid bulk modulus. The numerical solutions agree well with the attenuation and modulus dispersion measured in the laboratory, supporting the interpretation that attenuation and dispersion was due to FPD occurring in the heterogenous distribution of the coexisting fluids. The numerical simulations have the advantage that the models can easily be improved by including sub-core scale porosity and permeability distributions, which can also be determined using X-ray CT. In the future this could allow for conducting experiments on heterogenous samples.

How to cite: Chapman, S., Borgomano, J. V. M., Quintal, B., Benson, S. M., and Fortin, J.: Seismic wave attenuation due to fluid pressure diffusion at the mesoscopic scale: an experimental and numerical study, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16072,, 2021.

EGU21-9230 | vPICO presentations | EMRP1.15

Laboratory experiments of water injection coupled with ultrasonic monitoring reveal wave-induced fluid flow in microporous carbonate rock

Davide Geremia, Christian David, Christophe Barnes, Beatriz Menéndez, Jérémie Dautriat, Lionel Esteban, Joel Sarout, Sara Vandycke, and Fanny Descamps

Monitoring of fluid movements in the crust is one of the most discussed topics in oil & gas industry as well as in geothermal systems and CO2 storage, but still remains a challenge. The seismic method is one of the most common ways to detect the fluid migration. However, the use of ultrasonic monitoring at the sample scale in laboratory experiments persists as the most effective way to highlight large scale observations in which the boundary conditions are not well constrained.

To unravel the fluid effect on P-wave and S-wave velocity, we performed mechanical experiments coupled with ultrasonic monitoring on Obourg chalk from Mons basin (Belgium). Water injection tests under critical loading, imbibition tests and evaporation tests provided a full spectrum of observations of fluid-induced wave alteration in term of propagation time and attenuation.

The analysis of these experimental results showed that significant velocity dispersion and attenuation developed through variations in water saturation, and that these processes are linked to the presence of patches of water and air in the pore space.

We used the White’s formulation to model the relaxation effects due to spherical pockets of air homogeneously distributed in a water-saturated medium. In this framework, the pressure induced by the passing wave, produces a fluid flow across the water-air boundary with consequent energy loss.

This model reproduces both qualitatively and quantitatively the experimental results observed on the water injection tests. Indeed, it is shown that the progressive water saturation or desaturation of this chalk, generates a shift of the critical frequency (from the undrained relaxed towards unrelaxed regimes) which at some point matches the resonance frequency of the piezoelectric transducers used in the experimental setup (0.5 MHz). This phenomenon allowed us to get a continuous recording of the relaxation processes induced by saturation variations.

The outcomes of this work can significantly improve the actual knowledge on coupled effects of waves and fluids which is a crucial aspect of fluid monitoring in the context of reservoir evaluation and production.

How to cite: Geremia, D., David, C., Barnes, C., Menéndez, B., Dautriat, J., Esteban, L., Sarout, J., Vandycke, S., and Descamps, F.: Laboratory experiments of water injection coupled with ultrasonic monitoring reveal wave-induced fluid flow in microporous carbonate rock, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9230,, 2021.

EGU21-12197 | vPICO presentations | EMRP1.15

Capillary suction effects on surface-wave dispersion in partially saturated soils

Santiago Solazzi, Ludovic Bodet, Klaus Holliger, and Damien Jougnot

Estimation of water content in the shallow subsurface using seismic data is a complex task of increasing importance in the overall field of hydrogeophysics. In this context, the velocities of compressional (Vp) and shear (Vs) waves can be used to infer strong water content variations in unconsolidated soils, such as, the presence of the water table, by means of Vp/Vs ratio estimations. This approach, which is based on first-arrival time data, generally does not permit a proper quantification of the water content distribution in the partially saturated zone. Conversely, field and laboratory measurements indicate that surface-waves are indeed remarkably sensitive to both the water table depth and the saturation characteristics in the overlaying capillary fringe. This apparent difference in sensitivity between body and surface waves cannot be explained using conventional models. Observations and experiments show that the effective stress of unconsolidated porous media is not only affected by the overburden stress and pore pressures, as classic models assume, but also by capillary forces, which tend to stiffen the soil at relatively low saturations. In this work, we extend seminal rock physics models to include capillary suction effects in the effective stress of the soil. This approach provides effective elastic moduli and, thus, Vp and Vs, which are depth- and saturation-dependent. Then, we solve the quasi-static fluid flow equations in a porous medium and obtain saturation profiles for a given water table depth. This information, combined with the proposed rock physics model, permits to simulate simple seismic data sets, that is, body-wave first-arrival times and surface-wave phase velocities, for different water table depths and soil textures. Our results clearly show that capillary effects allow to explain the apparent difference in sensitivity between body- and surface-wave signatures in response to small water content variations in the partially saturated zone. Capillary effects are primarily relevant in porous media composed by relatively small characteristic grain sizes. We conclude that the proposed framework has the potential to fundamentally improve our characterization of near-surface environments using both active and passive seismic methods.

How to cite: Solazzi, S., Bodet, L., Holliger, K., and Jougnot, D.: Capillary suction effects on surface-wave dispersion in partially saturated soils, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12197,, 2021.

EGU21-12483 | vPICO presentations | EMRP1.15

Seismic anisotropy in metamorphic rocks from the COSC-1 borehole, Sweden: A cross-scale investigation from thin section analysis to seismic scales

Felix Kästner, Simona Pierdominici, Alba Zappone, Luiz F. G. Morales, Anja M. Schleicher, Franziska D. H. Wilke, and Christian Berndt

Metamorphic and deformed rocks in thrust zones show particularly high seismic anisotropy causing challenges for seismic imaging and interpretation. A good example is the Seve Nappe Complex in Jämtland, Sweden, an exhumed orogenic thrust zone characterized by a strong but incoherent seismic reflectivity and considerable seismic anisotropy. However, only little is known about the origin of the anisotropy in relation to composition, structural influences, and implications for measurements at different seismic scales. We present an integrative study of the seismic anisotropy at different scales combining mineralogical composition, microstructural analyses and seismic laboratory experiments from samples of the 2.5 km-deep COSC-1 borehole. While there is a pronounced crystallographic preferred orientation in most of the core samples, variations in anisotropy correlate strongly with bulk mineral composition and dominant core lithology. Based on three major lithologic different facies (felsic gneiss, amphibole-rich rocks, and mica schists), we propose an anisotropy model for the full length of the borehole, which indicates two prevailing anisotropic units. Comparison of laboratory seismic measurements and electron-backscatter diffraction (EBSD) data reveals a strong scale-dependence, which is more pronounced in the highly deformed, heterogeneous samples. This highlights the need for comprehensive cross-validation of microscale anisotropy analyses with additional lithological data when integrating seismic anisotropy through seismic scales.

How to cite: Kästner, F., Pierdominici, S., Zappone, A., Morales, L. F. G., Schleicher, A. M., Wilke, F. D. H., and Berndt, C.: Seismic anisotropy in metamorphic rocks from the COSC-1 borehole, Sweden: A cross-scale investigation from thin section analysis to seismic scales, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12483,, 2021.

EGU21-10820 | vPICO presentations | EMRP1.15

Estimation of fracture compliance based on group delays inferred from full-waveform sonic log data

Zhenya Zhou, Eva Caspari, Nicolás D. Barbosa, Andrew Greenwood, and Klaus Holliger

Fractures, which are ubiquitous in the Earth’s upper crust, have significant impacts on a wide range of human activities, and, hence, their adequate characterization is of wide interest and importance. Seismic methods have significant potential for effectively addressing this objective. When a seismic wave propagates across a fluid-filled fracture, its amplitude is diminished and its travel time is increased. Based on the linear slip theory, the associated amplitude decays and phase delays can be used to estimate the mechanical compliance of fractures. Full-waveform sonic (FWS) log data are particularly well-suited for this purpose. While the amplitudes of FWS data acquired during standard continuous logging runs (tool being moved uphole at a constant logging speed) can be somewhat unstable, the associated first-arrival travel times are generally quite robust. In this work, we exploit the relation between the time delay that seismic waves experience across fractures and relate them to the associated compliances. Specifically, we estimate fracture compliance from the differences in group time delay of the refracted P-wave between fractured and non-fractured sections along a borehole. Numerical simulations indicate that the proposed method provides reliable compliance estimates not only for individual fractures, but also for sets of multiple discrete fractures. This finding is corroborated by applying our approach to FWS log data acquired in the course of standard logging runs in the Bedretto Underground Laboratory ( Our estimates are comparable to previously inferred compliance values in a closely comparable geological environment (Grimsel test site, The latter were inferred under rather ideal conditions, involving the quasi-static acquisition of the FWS data as well as the combination of amplitude and travel time information for their interpretation. An interesting and important open question, which we plan to address in the following, concerns the influence of the heterogeneity of the host rock embedding the fractures on compliance estimation in general and on the proposed method in particular.

How to cite: Zhou, Z., Caspari, E., D. Barbosa, N., Greenwood, A., and Holliger, K.: Estimation of fracture compliance based on group delays inferred from full-waveform sonic log data, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-10820,, 2021.

EGU21-3466 | vPICO presentations | EMRP1.15

Constrained Magnetic Vector Inversion of Scanning Magnetic Microscopy Data for Modelling Magnetization of Multidomain Mineral Grains

Peter Lelièvre, Zeudia Pastore, Nathan Church, Madeline Lee, Hirokuni Oda, and Suzanne McEnroe

We are using 3D magnetic vector inversion (MVI) of scanning magnetic microscopy (SMM) data to investigate the fine‐scale magnetization of rock samples, and particularly of their remanence carriers, which can record geologically meaningful information. Previous investigations of magnetite grains suggest variable remanence intensities and directions coherent with multidomain behaviour. This research seeks to improve our understanding of the contribution of different microstructures on remanence acquisition.

SMM offers a spatial resolution down to tens of micrometers, allowing detailed investigation of discrete magnetic mineral grains, or magnetic textures and structures. However, all magnetic measurements are, at some scale, bulk measurements. Further analysis of the data is required to extract information about the magnetization within the samples: for this, we employ state-of-the-art MVI methods. The MVI problem suffers from a high degree of nonuniqueness. Additional constraints are required to obtain accurate, reliable and interpretable results. Such constraints are readily available for this application.

SMM instruments use magnetic shields or Helmholtz coils to allow collection of data in controlled magnetic fields, enabling the removal of induced magnetization effects. Measurements can be taken both above and below the sample. Individual magnetized mineral grains are easily outlined through optical and electron microscopy. The internal geometry of the oxide mineral phases and compositions can also be constrained. Physical property information constrains the range of magnetization intensity. As such, there is a tremendous amount of constraining information invaluable for reducing the nonuniqueness of the inverse problem. We use a highly flexible and functional inversion software package, MAGNUM, developed jointly at Mount Allison University and Memorial University of Newfoundland, that allows incorporation of all available constraints.

We take a multitiered approach for investigating specific magnetized grains. First, coarse regional inversions are performed to assess and remove any effects of other magnetized grains in the vicinity. The entire grain is then modelled with a homogeneous magnetization to obtain an approximate but representative bulk magnetization. The grain is then modelled as a collection of independent subdomains, each with a different homogeneous magnetization direction. Subsequently, more heterogeneous scenarios are considered by relaxing inversion constraints until the data can be fit to the desired degree.

Obtaining reliable information about the magnetic mineralogy of rock samples is vital for an understanding of the origin of rock bulk behaviour in both the laboratory and larger scale magnetic surveys. This work is among the first to simultaneously invert SMM data collected above and below a thin sample, which is critical for improving depth resolution on thicker samples. It is also the first time we have been able to incorporate all available constraints into inverse modelling to improve results.

How to cite: Lelièvre, P., Pastore, Z., Church, N., Lee, M., Oda, H., and McEnroe, S.: Constrained Magnetic Vector Inversion of Scanning Magnetic Microscopy Data for Modelling Magnetization of Multidomain Mineral Grains, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-3466,, 2021.

EGU21-5887 | vPICO presentations | EMRP1.15

Seismic moduli dispersion and attenuation obtained using Digital Rock Physics

Simón Lissa, Matthias Ruf, Holger Steeb, and Beatriz Quintal

Seismic waves are affected by rock properties such as porosity, permeability, grain material and by their heterogeneities as well as by the fluid properties saturating the rocks. Consequently, seismic methods are a valuable tool for the indirect characterization of rocks. For example, at the microscale, the presence of compliant pores (cracks or grain contacts) in fluid-saturated rocks can cause strong seismic attenuation and velocity dispersion. In this case, the deformation caused by a passing wave induces a fluid pressure gradient between compressed compliant pores and much less compressed pores (stiff isometric pores or cracks having a different orientation than the most compressed ones) if they are hydraulically connected. The consequent fluid pressure diffusion (FPD) dissipates seismic energy due to viscous friction in the fluid.

Digital rock physics (DRP) aims to reproduce experimental measurements using numerical simulation in models derived from high resolution rock images. We developed a DRP workflow to calculate the frequency dependent seismic moduli dispersion and attenuation in fluid-saturated models derived from micro X-Ray Computed Tomography (µXRCT) images. Filtering, segmentation and meshing procedures are applied on sub-volumes of different rock images to create 3D numerical models. We apply our workflow to calculate seismic moduli attenuation due to FPD at the microscale (squirt flow). We consider a µXRCT image of a cracked (through thermal treatment) Carrara marble sample. A detailed visualization of the fluid pressure as well as of the energy dissipation rate in the 3D model helps to understand the squirt flow attenuation process at different frequencies.

How to cite: Lissa, S., Ruf, M., Steeb, H., and Quintal, B.: Seismic moduli dispersion and attenuation obtained using Digital Rock Physics, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5887,, 2021.

EGU21-6492 | vPICO presentations | EMRP1.15

Assessing elastic properties of cracks in rock samples subjected to thermo-hydro-chemo-mechanical damage

Anthony Clark, Tiziana Vanorio, Andrey Radostin, and Vladimir Zaitsev

An understanding of micro- and macrofracture behavior in low porosity rocks is pertinent to several areas of energy and environmental science such as petroleum production, carbon sequestration, and enhancement of technologies based on geothermal energy, etc. For example, the carbonate reservoirs in dolomitic or micritic formations with matrix porosities below 6% suggest the importance of fracture-augmented permeability in production. Similarly, hydrocarbons have been found on nearly every continent in tight basement rocks, all of which have little matrix porosity and their permeability therefore rely solely on hydraulic connectivity from fractures. For geothermal energy, various igneous and sedimentary rocks (granites, basalts, and limestones) are being exploited across the globe, with some of the lowest porosity and permeability. In all these cases, fractures are necessary to improve rock permeability and thermal exchange between the rock and working fluid, which can be enabled by hydraulic stimulation, as well as by secondary cracking due to extreme temperature gradients from the injection of cold water. The fracture geometry, density, and distribution all together control not only fluid and thermal transport in the rocks, but also their seismic attributes that can be used to extract information about the fractures. 
In order to accurately interpret the seismo-acoustic data (usually, the velocities of compression and shear waves) reliable rock physics models are required. Here, we report the results of interpretation of such experimental data for both as-cored rock samples and those subjected to thermo-hydro-chemo-mechanical damage (THCMD) in the laboratory. For interpretation, we use a convenient model of fractured rock in which fractures are represented as planar defects with decoupled shear and normal compliances. The application of such an approach makes it possible to assess and compare the elastic properties of fractures in the rocks before and after application of THCMD procedures. For the analyzed samples of granites, basalts, and limestones it has been found that for a significant portion of rocks, the ratio of normal-to-shear compliances of cracks significantly differ from the value typical of conventionally assumed penny-shape cracks. Furthermore, for some samples, this ratio appears to be noticeably different for fractures existing in the as-cored rock and arising in the same rock after THCMD procedures. These results indicate that damage to a rock typically changes its compliance ratio since the old and new cracks are likely to have different elastic properties. Our results are also consistent with the notion that a specific damage process occurring for a given microstructure will consistently create cracks with a particular set of elastic properties. The proposed methodology for assessment of elastic properties of cracks in rock samples subjected to thermo-hydro-chemo-mechanical damage has given previously inaccessible useful information about the elastic properties of fractures and can be extended to interpretation of seismic attributes of rocks for a broad range of other applications.

How to cite: Clark, A., Vanorio, T., Radostin, A., and Zaitsev, V.: Assessing elastic properties of cracks in rock samples subjected to thermo-hydro-chemo-mechanical damage, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-6492,, 2021.

EGU21-4147 | vPICO presentations | EMRP1.15

Digital rock physics: A geological driven segmentation workflow for Ruhr sandstone

Martin Balcewicz, Mirko Siegert, Marcel Gurris, David Krach, Matthias Ruf, Holger Steeb, and Erik H. Saenger

Over the last two decades, Digital Rock Physics (DRP) has become a complementary part of the characterization of reservoir rocks due to, among other things, the non-destructive testing character of this technique. The use of high-resolution X-ray Computed Tomography (XRCT) has become widely accepted to create a digital twin of the material under investigation. Compared to other imaging techniques, XRCT technology allows a location-dependent resolution of the individual material particles in volume. However, there are still challenges in assigning physical properties to a particular voxel within the digital twin, due to standard histogram analysis or sub-resolution features in the rock. For this reason, high-resolution image-based data from XRCT, transmitted-light microscope, Scanning Electron Microscope (SEM) as well as inherent material properties like porosity are combined to obtain an optimal spatial image of the studied Ruhr sandstone by a geologically driven segmentation workflow. On the basis of a homogeneity test, which corresponds to the evaluation of the grayscale image histogram, the preferred scan sample sizes in terms of transport, thermal, and effective elastic rock properties are determined. In addition, the advanced numerical simulation results are compared with laboratory tests to provide possible upper limits for sample size, segmentation accuracy, and a calibrated digital twin of the Ruhr sandstone. The comparison of representative grayscale image histograms as a function of sample sizes with the corresponding advanced numerical simulations, provides a unique workflow for reservoir characterization of the Ruhr sandstone.

How to cite: Balcewicz, M., Siegert, M., Gurris, M., Krach, D., Ruf, M., Steeb, H., and Saenger, E. H.: Digital rock physics: A geological driven segmentation workflow for Ruhr sandstone, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-4147,, 2021.

EGU21-4691 | vPICO presentations | EMRP1.15

Numerical determination of pressure-dependent effective thermal conductivity in Berea sandstone

Mirko Siegert, Marcel Gurris, and Erik Hans Saenger

Within the scope of the present work, the pressure-dependent effective thermal conductivity of rock samples is simulated. Our workflow can be assigned to the field of digital rock physics. In a first step, a 3D micro-CT scan of a rock sample is taken. Subsequently, the resulting greyscale images are analysed and segmented depending on the occurring phases. Based on this data set, a computational mesh is created and the corresponding thermal conductivities are assigned to each phase. Finally the numerical simulations can be carried out.
For the representation of the pressure dependency we use the approach proposed by Saenger [1]. By making use of the watershed algorithm, boundaries between the individual grains of the rock sample are detected and assigned to an artificial contact phase. In the course of several simulations, the thermal conductivity of the contact phase is continuously increased. Starting with the thermal conductivity of the pore phase and ending with the thermal conductivity of the grain phase. A linear correlation is used to match the thermal conductivity of the contact phase with the pressure of a given experimental data set. This enables a direct comparison between simulation and measurement.
In a further step, the numerical model is calibrated to optimise the agreement between experimental data and simulation results. In particular, starting from two calibration points of the experimental data set, an adjustment of the thermal conductivities in the numerical model is carried out. While the thermal conductivity of the pore phase is held constant during the whole calibration process, thermal conductivities of the grain and contact phase are adjusted.

[1] Saenger et al. 2016. Analysis of high-resolution X-ray computed tomography images of Bentheim sandstone under elevated confining pressures. Geophysical Prospecting, 64(4), 848–859.


How to cite: Siegert, M., Gurris, M., and Saenger, E. H.: Numerical determination of pressure-dependent effective thermal conductivity in Berea sandstone, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-4691,, 2021.

EGU21-14273 | vPICO presentations | EMRP1.15

Textural and mineralogical controls on temperature dependent SIP behavior during freezing and thawing

Jonas K. Limbrock, Maximilian Weigand, and Andreas Kemna

Geoelectrical methods are increasingly being used for non-invasive characterization and monitoring of permafrost sites, since the electrical properties are sensitive to the phase change of liquid to frozen water. Here, electrical resistivity tomography (ERT) is most commonly applied, using resistivity as a proxy for various quantities, such as temperature or ice content. However, it is still challenging to distinguish between air and ice in the pore space of the rock based on resistivity alone due to their similarly low electrical conductivity. Meanwhile, geoelectrical methods that utilize electrical polarization effects to characterize permafrost are also being explored. For example, the usage of the spectral induced polarization (SIP) method, in which the complex, frequency-dependent impedance is measured, can reduce ambiguities in the subsurface conduction properties, considering the SIP signature of ice. These measurements seem to be suitable for the quantification of ice content (and thus the differentiation of ice and air), and for the improved thermal characterization of alpine permafrost sites. However, to improve the interpretation of SIP measurements, it is necessary to understand in more detail the electrical conduction and polarization properties as a function of temperature, ice content, texture, and mineralogy under frozen and partially frozen conditions.

In the study presented here, electrical impedance was measured continuously using SIP in the frequency range of 10 mHz to 45 kHz on various water-saturated solid rock and loose sediment samples during controlled freeze-thaw cycles (+20°C to -40°C). These measurements were performed on rock samples from different alpine permafrost sites with different mineralogical compositions and textures. For all samples, the resistance (impedance magnitude) shows a similar temperature dependence, with increasing resistance for decreasing temperature. Also, hysteresis between freezing and thawing behavior is observed for all measurements. During freezing, a jump within the temperature-dependent resistance is observed, suggesting a lowering of the freezing point to a critical temperature where an abrupt transition from liquid water to ice occurs. During thawing, on the other hand, there is a continuous decrease in the measured resistance, suggesting a continuous thawing of the sample. The spectra of impedance phase, which is a measure for the polarization, exhibit the same qualitative, well-known temperature-dependent relaxation behaviour of ice at higher frequencies (1 kHz - 45 kHz), with variations in shape and strength for different rock texture and mineralogy. At lower frequencies (1 Hz - 1 kHz), a polarization with a weak frequency dependence is observed in the unfrozen state of the samples. We interpret this response as membrane polarization, which likewise depends on the texture as well as on the mineralogy of the respective sample. This polarization response partially vanishes during freezing. Overall, the investigated SIP spectra do not only show a dependence on texture and mineralogy, but mainly a dependence on the presence of ice in the sample as well as temperature. This indicates the possibility of a thermal characterization, as well as a determination of the ice content, of permafrost rocks using SIP.

How to cite: Limbrock, J. K., Weigand, M., and Kemna, A.: Textural and mineralogical controls on temperature dependent SIP behavior during freezing and thawing, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14273,, 2021.

EGU21-13527 | vPICO presentations | EMRP1.15

Correlations of magnetic pore fabrics with pore fabrics derived from high-resolution X-ray computed tomography and with permeability anisotropy in sedimentary rocks and synthetic samples

Yi Zhou, Michele Pugnetti, Anneleen Foubert, Pierre Lanari, Christoph Neururer, and Andrea Biedermann

Magnetic pore fabrics (MPF) are an indirect measure of the 3D pore structure. They are defined by measuring anisotropy of magnetic susceptibility after samples have been impregnated with ferrofluid. Previous studies proposed that MPFs target pores down to 10 nm. Therefore, the method complements X-ray computed tomography (XRCT) datasets, with resolution on the order of 1-10 µm. Empirical relationships exist between MPF and pore fabric, and between MPF and permeability anisotropy. This study investigates quantitative correlations between these three properties, and between measured quantities and digital-rock-model-simulations of permeability anisotropy and MPF. Samples used for this study include natural sedimentary rocks and synthetic samples. Sediments are Plio-Pleistocene calcarenite (Apulia, Italy) with ~50% porosity and complex pore structure, and Upper Marine Molasse sandstone (Belpberg, Switzerland) with 10-20% porosity and relatively homogeneous pore space properties. Synthetic samples were made from quartz sand and calcite powder in different proportions, to simulate sandstone and carbonate rocks. Samples were characterized by pycnometry, XRCT scans, MPF determination and directional permeability measurements to obtain porosity, digital rock models, MPFs and permeability anisotropy. Porosity, permeability anisotropy, and MPFs were also computed based on digital rock models derived from XRCT data, and compared to direct measurements. Permeability anisotropy and MPF are both second-order tensors, representing the average property of the entire sample. To directly relate the XRCT-derived individual pore properties to these second-order tensor quantities, a total shape ellipsoid was computed by adding the second-order tensors reflecting the best-fit ellipsoids of single pores. Once all properties were described by second-order tensors, they were correlated in terms of fabric orientation, degree and shape of anisotropy. The MPF and total shape ellipsoids are coaxial when the samples have sufficiently large pores to be resolved, and good impregnation efficiency, and as expected, total shape ellipsoids have larger anisotropy degree. Preliminary results further indicate that the permeability anisotropy is partly consistent with total shape ellipsoids and MPFs. The defined quantitative relationships facilitate the interpretation of MPF data, thus making the method more applicable to geological and fluid migration studies.

How to cite: Zhou, Y., Pugnetti, M., Foubert, A., Lanari, P., Neururer, C., and Biedermann, A.: Correlations of magnetic pore fabrics with pore fabrics derived from high-resolution X-ray computed tomography and with permeability anisotropy in sedimentary rocks and synthetic samples, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-13527,, 2021.

The magnetic pore fabrics method is a useful technique to investigate the pore fabric of rocks. The method is based on impregnating porous samples with ferrofluid, a colloidal suspension of magnetic nanoparticles (particle size of 10 nm) in water or oil carrier fluid, and measuring the anisotropy of magnetic susceptibility. It succesfully provides the average pore shape and orientation. This information is important to determine the preferred direction of fluid flow. A crucial step in magnetic pore fabric studies is ferrofluid impregnation; several studies pointed out the importance of forced impregnation methods to enhance impregnation efficiency. This study compares directional impregnation techniques (applying external forces along the ferrofluid flow along the axis of the core sample) with standard vacuum- and immersion-based methods. The newly developed or adapted techniques include: (1) pressure experiment: a cylindrical sample is placed in a metal tube under confining pressure of 12 bar, and an external pump-syringe system injects ferrofluid at a constant rate of 100 ml/min that generates a differential pressure of 5 bar; (2) resin flowthrough: vacuum is applied at the bottom of the sample and a mixture of resin and oil-based ferrofluid supplied at the top, so that the resin drags the fluid into the pores, where it hardens; (3) magnetically assisted flowthrough: fluid flow is enforced by the combined action of a hydraulic pressure gradient in the ferrofluid reservoir (~10 kPa) and the magnetic force exerted by the field gradient of about 2 A/m2 in the vicinity of an electric coil. These impregnation methods were tested on natural and synthetic samples, for which previous experiments employing standard impregnation methods exist. The natural samples include calcarenite from Apulia, Italy (50% porosity) and sandstone from Schupfheim, Swiss molasse (20% porosity). Synthetic samples consist of calcite and quartz sand in different proportions, consolidated with liquid glass (sodium silicate) in a cubic consolidation cell (specifically designed for the experiment), applying uniaxial pressure along the z axis, to create uniaxial anisotropy. The cube was dried in the oven for three days and three cylindrical cores were drilled along the x, y and z axes. For each impregnation method, the magnetic anisotropy of the samples was measured before and after impregnation. Impregnation efficiency was tested  using bulk susceptibility measurements, visual microscopic investigations and susceptibility profiles along the flow direction. Initial results show that (1) directional forced impregnation is more efficient than traditional methods in impregnating smaller pores,  avoids particle aggregation, and allows viscous fluid such as resin to acess the sample’s pores; (2) directional impregnation methods require less  fluid;  (3) the distribution of the ferrofluid after impregnation is more uniform, overcoming the difficulty of impregnating the  centre of the sample;  and (4) the fluid flow rate must be faster than the particle aggregation rate. For future studies, directional forced impregnation systems are recommended over standard vacuum- and immersion-based impregnation methods.


How to cite: Pugnetti, M., Zhou, Y., and Biedermann, A.: Experimental improvements for ferrofluid impregnation of rocks using directional forced impregnation methods: results on natural and synthetic samples, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14695,, 2021.

EMRP1.16 – Open Session in Rock Physics

EGU21-5053 | vPICO presentations | EMRP1.16

Seismic attenuation and rheology of crustal rocks: results from numerical tests

Maria Aurora Natale Castillo, Magdala Tesauro, and Mauro Cacace

Seismic attenuation of the rocks mainly depends on their intrinsic anelasticity, which is the dissipation of seismic energy as it propagates through the medium. Several studies have already demonstrated that seismic attenuation, described by the Q-factor, is intrinsically related to the rocks’ viscosity, considering their common dependency on composition, grain size, fluid content, and T-P conditions. However, viscous deformation of the rocks occurs through different mechanisms: diffusion creep, numerous mechanisms of the dislocation creep, pressure solution, which are expressed by several Arrhenius-type constitutive laws. This makes more complex the investigation of quantitative relationships between seismic attenuation and viscous rocks' rheology.

The main purpose of this study is to investigate the mutual dependence of the seismic attenuation and viscous deformation of the crustal rocks. To this aim, we performed several numerical tests to check the variability of some physical properties (e.g., elastic modulus, Poisson’s ratio) and improve the existing relationships between seismic attenuation and viscous deformation of several natural rock samples. The Burgers mechanical model and Arrhenius relation are included in these test series to achieve a closer approximation of the rocks’ viscous deformation.

In this way, it will be possible to predict the viscous rheology of rocks from the laws describing the seismic attenuation and viceversa. The obtained results will be used to (1) constrain the Q-factor and rheological (creep) parameters, which are still subjected to high uncertainties, (2) validate/modify the existing seismic attenuation and rheological laws, (3) increase the robustness of the geodynamic and rocks’ mechanics numerical codes and our understanding of the role that rocks’ rheology exerts on the tectonic processes.

How to cite: Natale Castillo, M. A., Tesauro, M., and Cacace, M.: Seismic attenuation and rheology of crustal rocks: results from numerical tests, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5053,, 2021.

EGU21-5585 | vPICO presentations | EMRP1.16

Time-lapse synchrotron X-ray imaging of deformation modes in organic-rich Green River Shale heated under confinement 

Maya Kobchenko, Anne Pluymakers, Benoit Cordonnier, Nazmul Mondol, and Francois Renard

Shales are layered sedimentary rocks, which can be almost impermeable for fluids and act as seals and cap-rocks, or if a shale layer hosts a fracture network, it can act as a fluid reservoir and/or conduit. Organic-rich shales contain organic matter - kerogen, which can transform from solid-state to oil and gas during burial and exposure to a suitable temperature. When hydrocarbons are expelled from the organic matter due to maturation, pore-pressure increases, which drives the propagation of hydraulic fractures, a mechanism identified to explain oil and gas primary migration. Density, geometry, extension, and connectivity of the final fracture network depend on the combination of the heating conditions and history of external loading experienced by the shale. Here, we have performed a series of rock physics experiments where organic-rich shale samples were heated, under in situ conditions, and the development of microfractures was imaged through time. We used the high-energy X-ray beam produced at the European Synchrotron Radiation Facility to acquire dynamic microtomography images and monitor different modes of shale deformation in-situ in 3D. We reproduced natural conditions of the shale deformation processes using a combination of axial load, confining pressure, and heating of the shale samples. Shales feature natural sedimentary laminations and hydraulic fractures propagate parallel to these laminae if no overburden stress is applied. However, if the principal external load becomes vertical, perpendicular to the shale lamination, the fracture propagation direction can deviate from the horizontal one. Together horizontal and vertical fractures form a three-dimensional connected fracture network, which provides escaping pathways for generated hydrocarbons. Our experiments demonstrate that tight shale rocks, which are often considered impermeable, could host transient episodes of micro-fracturing and high permeability during burial history.

How to cite: Kobchenko, M., Pluymakers, A., Cordonnier, B., Mondol, N., and Renard, F.: Time-lapse synchrotron X-ray imaging of deformation modes in organic-rich Green River Shale heated under confinement , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-5585,, 2021.

EGU21-11531 | vPICO presentations | EMRP1.16

Measurements of anisotropic poroelastic moduli in Westerly granite

Bobby Elsigood, Nicolas Brantut, Philip Meredith, Tom Mitchell, and David Healy

We measured both drained (dry and saturated) and undrained elastic and poro-elastic moduli in a sample of thermally cracked Westerly granite. Measurements were made at nominal effective confining pressures of 30 MPa and 100 MPa. Differential stress was also applied to the sample in order to induce increasing anisotropy (transverse isotropy). Radial and axial stress were independently cycled at increasing levels of differential stress to enable the measurement of anisotropic moduli in a sample with increasing anisotropy. Measurements were first conducted under dry conditions in order to obtain the drained dynamic moduli before the sample was saturated with water. We used newly developed, miniature differential pressure transducers were located directly around the sample surface, which allowed for direct measurement of Skempton's radial coefficient Bx and axial coefficient Bz. Wavespeeds measured on the dry sample were inverted to obtain the compliances Cijkl. The coefficients Bx and Bz were then calculated indirectly from the compliances using poro-elasticy theory for comparison with the direct measurements. Our results show that the values of Bz calculated from the compliances compare well with those measured directly, but that the calculated values of Bx significantly overestimate those measured directly, particularly at an effective pressure of 30 MPa. The discrepancy between the direct measurements and those obtained from drained moduli is likely due to nonlinear elastic effects, such as crack closure and opening, which occur during the small but finite pressure and stress steps.

How to cite: Elsigood, B., Brantut, N., Meredith, P., Mitchell, T., and Healy, D.: Measurements of anisotropic poroelastic moduli in Westerly granite, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-11531,, 2021.

EGU21-11824 | vPICO presentations | EMRP1.16

Prediction of Petrophysical Properties from Seismic Inversion and Neural Network: A case study

Siddharth Garia, Arnab Kumar Pal, Karangat Ravi, and Archana M Nair

Seismic inversion method is widely used to characterize reservoirs and detect zones of interest, i.e., hydrocarbon-bearing zone in the subsurface by transforming seismic reflection data into quantitative subsurface rock properties. The primary aim of seismic inversion is to transform the 3D seismic section/cube into an acoustic impedance (AI) cube. The integration of this elastic attribute, i.e., AI cube with well log data, can thereafter help to establish correlations between AI and different petrophysical properties. The seismic inversion algorithm interpolates and spatially populates data/parameters of wells to the entire seismic section/cube based on the well log information. The case study presented here uses machine learning-neural network based algorithm to extract the different petrophysical properties such as porosity and bulk density from the seismic data of the Upper Assam basin, India. We analyzed three different stratigraphic  units that are established to be producing zones in this basin.

 AI model is generated from the seismic reflection data with the help of colored inversion operator. Subsequently, low-frequency model is generated from the impedance data extracted from the well log information. To compensate for the band limited nature of the seismic data, this low-frequency model is added to the existing acoustic model. Thereafter, a feed-forward neural network (NN) is trained with AI as input and porosity/bulk density as target, validated with NN generated porosity/bulk density with actual porosity/bulk density from well log data. The trained network is thus tested over the entire region of interest to populate these petrophysical properties.

Three seismic zones were identified from the seismic section ranging from 681 to 1333 ms, 1528 to 1575 ms and 1771 to 1814 ms. The range of AI, porosity and bulk density were observed to be 1738 to 6000 (g/cc) * (m/s), 26 to 38% and 1.95 to 2.46 g/cc respectively. Studies conducted by researchers in the same basin yielded porosity results in the range of 10-36%. The changes in acoustic impedance, porosity and bulk density may be attributed to the changes in lithology. NN method was prioritized over other traditional statistical methods due to its ability to model any arbitrary dependency (non-linear relationships between input and target values) and also overfitting can be avoided. Hence, the workflow presented here provides an estimation of reservoir properties and is considered useful in predicting petrophysical properties for reservoir characterization, thus helping to estimate reservoir productivity.

How to cite: Garia, S., Pal, A. K., Ravi, K., and Nair, A. M.: Prediction of Petrophysical Properties from Seismic Inversion and Neural Network: A case study, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-11824,, 2021.

EGU21-12261 | vPICO presentations | EMRP1.16

Using AE based Machine Learning Approaches to Forecast Rupture during Rock Deformation Laboratory Experiments

Sergio Vinciguerra, Thomas King, and Philip Benson

The ability to detect precursors of dynamic failure in brittle rocks has key implications for hazard forecasting at the field scale. In recent years, laboratory scale rock deformation experiments are providing a wealth of information on the physics of the fracture process ranging from fracture nucleation, crack growth and damage accumulation, to crack coalescence and strain localization. Parametric analysis of laboratory Acoustic Emission (AE) data has revealed periodic trends and precursory behaviour of the rupture source mechanisms as a fault zone enucleates and develops, suggesting these processes are somehow repeatable and forecastable. However, due to the inherent anisotropy of rock media and the range of environmental conditions in which deformation occurs, finding full consistency between AE datasets and a prediction of rupture mechanisms from AE analysis is still an open goal. Here we apply a Time Delay Neural Network (TDNN) to Acoustic Emission (AE) sets recorded during conventional triaxial rock deformation tests. We forecast the Time-to-Failure using the discrete, non-continuous timeseries of AE rate, amplitude, focal mechanism and forward scattering properties. 4x10 cm samples of Alzo granite, a homogeneous medium-grained plutonic rock from NW Italy with an initial porosity as low as 0.72%, were triaxially deformed at strain rates of 3.6mm/hr under dry conditions until dynamic failure at confining pressures of 5, 10, 20 and 40 MPa respectively. Each sample was positioned inside an engineered rubber jacket fitted with ports where an array of twelve 1 MHz single-component Piezo-Electric Transducers were embedded, allowing to record AE during the experimentation. Several parameters were considered for the TDNN training: AE rate, deformation stages prior failure (elasticity, inelasticity and coalescence), AE amplitude, source mechanisms and scattering. All these parameters are key indicators of the evolving damage in the medium. Our training input consists of simplified timeseries of the previously discussed AE parameters from the experiments carried out at the lowest confining pressure (5 MPa). The inputs are classified as the stress-until failure and strain-until-failure for each AE. Once trained we then simulate the model on the untrained datasets to test it as a forecasting tool at higher confinements. At each step the model is simulated on AE data from the previous 0.2% of strain. At 10 MPa we observe a reliable forecast of failure that starts with the anelastic phase and becomes more accurate during strain-softening. At higher confining pressure, an increased limit of forecasting the solution is observed and interpreted with more complexity in the coalescence process. Despite these limitations, the model shows that when trained even on a limited input it is able to forecast dynamic failure in unseen data with surprising accuracy. Future studies should investigate AE spatial distribution for the TDNN training.

How to cite: Vinciguerra, S., King, T., and Benson, P.: Using AE based Machine Learning Approaches to Forecast Rupture during Rock Deformation Laboratory Experiments, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12261,, 2021.

EGU21-12490 | vPICO presentations | EMRP1.16

Electrical-elastic modelling of rocks using cross-property DEM

Phillip Cilli and Mark Chapman

The combination of electrical and elastic measurements can be useful in lowering subsurface characterisation uncertainty. The majority of rock physics relations which relate a rock’s electrical and elastic properties, however, rely on the estimation of porosity as an intermediate step. By combining differential effective medium schemes which relate a rock’s electrical and elastic properties to porosity and pore shape, we obtain cross-property expressions which are independent of porosity, depending only on pore aspect ratio. Analysing published joint electrical-elastic measurements shows the cross-property model works well for clean sandstones, and models Vp/Vs ratios as a function of resistivity without porosity. Although clay-bearing sandstones are more complex, our model can still identify the correct trends. On theoretical grounds, it seems our approach has the potential to produce additional cross-property relations; a topic upon which we speculate.

How to cite: Cilli, P. and Chapman, M.: Electrical-elastic modelling of rocks using cross-property DEM, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12490,, 2021.

EGU21-15915 | vPICO presentations | EMRP1.16

Is breaking through matter a hot matter? How to predict material failure by monitoring creep

Renaud Toussaint, Tom Vincent-Dospital, Alain Cochard, Eirik Grude Flekkøy, and Knut Jørgen Måløy

In any domain involving some stressed solids, that is, from seismology to rock physics or general engineering, the strength of matter is a paramount feature to understand.  The global failure of a mechanically loaded solid is usually dictated by the growth of its internal micro-cracks and dislocations. When this growth is rather smooth and distributed, the solid is considered to be in ductile condition. Alternatively, an abrupt propagation of localized defects leads to a brittle rupture of the full matrix.
It is then critical to understand what the physics and dynamics of isolated cracks are, when their tips are loaded at a given stress level. While the general elasticity theory predicts such stress to diverge, it is  acknowledged that some area around the crack fronts is rather plastic. In other words, some dissipation of mechanical energy, in a so-called process zone around a crack tip, prevents the - unphysical - stress divergence and shields the fronts from excessive load levels.

In this work, we focus on the local Joule heating, that significantly contributes to the energy dissipation. Analysing experimental data of the rupture of many materials, we indeed show that the scale for the thermal release around crack tips explains why the toughness of different media spans over orders of magnitude (we analysed materials spanning over 5 decades of energy release rate), whereas the covalent energy to separate two atoms does not.

We here discuss the ability of this simple thermally activated sub-critical model, which includes the auto-induced thermal evolution of crack stips [1], to predict the catastrophic failure of a vast range of materials [2]. It is in particular shown that the intrinsic surface energy barrier, for breaking the atomic bonds of many solids, can be easily deduced from the slow creeping dynamics of a crack. This intrinsic barrier is however higher than the macroscopic load threshold at which brittle matter brutally fails, possibly as a result of thermal activation and of a thermal weakening mechanism. We propose a novel method to compute the macroscopic critical energy release rate of rupture, Gc macroscopic, solely from monitoring slow creep, and show that this reproduces the experimental values within 50% accuracy over twenty different materials (such as glass, rocks, polymers, metals), and over more than four decades of fracture energy. We also infer the characteristic energy of rupturing bonds, and the size of an intense heat source zone around crack tips, and show that it scales as the classic process zone size, but is significantly (105 to 107 times) smaller.


[1] Vincent-Dospital, T., Toussaint, R., Santucci, S., Vanel, L., Bonamy, D., Hattali, L., Cochard,  A, Flekkøy, E.G. and Måløy, K.J. (2020). How heat controls fracture: the thermodynamics of creeping and avalanching cracks. Soft Matter, 2020, 16, 9590-9602. DOI: 10.1039/D0SM01062F

[2] Vincent-Dospital, T., Toussaint, R., Cochard, A., Flekkøy, E. G., & Måløy, K. J. (2020). Is breaking through matter a hot matter? A material failure prediction by monitoring creep. arXiv preprint arXiv:2007.04866.

How to cite: Toussaint, R., Vincent-Dospital, T., Cochard, A., Flekkøy, E. G., and Måløy, K. J.: Is breaking through matter a hot matter? How to predict material failure by monitoring creep, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15915,, 2021.

EGU21-15948 | vPICO presentations | EMRP1.16

High-speed tensile fractures in granite and gneiss: an experimental study

Beno J Jacob, Santanu Misra, Venkitanarayanan Parameswaran, and Nibir Mandal

Tensile fractures are ubiquitous in impact structures formed because of high strain rate deformations of the earth’s crust. At regions far from the point of meteorite impact, intense rupturing, fragmentation, and pulverisation are an implication of pressure waves limiting at the tensile strength of the host rock with little influence of shock deformation or shear failure. The branching and anastomosing of the fractures are controlled by the local stress state and anisotropy. Thus, a network of infilled fractures or impact breccia dikes is a common feature in the subsurface of impact sites.

We have investigated the failure processes under high strain rates responsible for the formation of Mode-I breccia dikes, at the laboratory scale. The control of planar fabric structures in the development of anastomosing tensile fracture networks was studied through high-strain-rate Brazilian disc tests on gneiss (foliated) and granite (isotropic) samples. A Split Hopkinson Pressure Bar, equipped with high-speed photography (~105 fps), was employed in the study. The gneissic foliation in the gneiss samples were oriented at θ = 0, 45 and 90° to the compression direction. The strength of granite lies between 24 and 26 MPa, and the gneisses failed in the range of 29-37MPa at about 70-90 μs. The fracture network formation was seen in the time series images. There is a stark disparity in the nature of failure of granite from gneiss and the geometry of clasts formed in each rock type. While granite samples fail with pulverised clasts localised along a single fracture spanning the diameter of the sample along the compression direction, the gneisses further developed a network of secondary fractures forming large elongate clasts. Preferential orientation of secondary crack growth in relation to the foliation is strongly influenced by θ in gneiss samples. The aspect ratio of the pulverised clasts (size < 10mm) formed in granite was about 1:2, whereas the gneisses produced larger clasts. The clasts in gneisses had an aspect ratio of 1:4 for θ = 45 and 90º, and 1:5 for θ = 0º.

The branching and anastomosing nature of fractures is similar in fracture networks observed from the field and in the experiments, thus providing an insight into the formation of high-speed impact breccia dikes in isotropic and foliated rocks. Our experiments demonstrate that monomict breccia dikes may by formed in situ inclusive of clasts, rather than by infilling in previously formed tensile fractures.

How to cite: Jacob, B. J., Misra, S., Parameswaran, V., and Mandal, N.: High-speed tensile fractures in granite and gneiss: an experimental study, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15948,, 2021.

EMRP1.18 – Coupled thermo-hydro-mechanical-chemical (THMC) processes in geological media

EGU21-14418 | vPICO presentations | EMRP1.18

Thermo-mechanical Effects on Fracture Growth and Apertures in Three-Dimensional Subsurface Fractured Rocks 

Adriana Paluszny, Robin N. Thomas, M. Cristina Saceanu, and Robert W. Zimmerman

A finite-element based, quasi-static growth algorithm models mixed mode concurrent fracture growth in three dimensions, leading to the formation of non-planar arrays and networks. To model the fully coupled THM model, equations describing mechanical deformation as well as heat transfer in the matrix and in the fractures are introduced in the formulation, simultaneously accounting for the effect of fluid flow and stress-strain response. This results in five separate, but two-way coupled model equations: a thermoporoelastic mechanical model; two fluid flow equations, one for the rock matrix and one for the fractures; two heat transfer equations, similarly for both the matrix and fractures. Fractures are represented explicitly as discrete surfaces embedded within a volumetric domain [1]. Growth is computed as a set of vectors that modify the geometry of a fracture by accruing new fracture surfaces in response to brittle deformation. Fracture tip stress intensity factors drive fracture growth. This growth methodology is validated against analytical solutions for fractures under compression and tension [2]. Thermal effects on the apertures and growth patterns will be presented. Isolated fracture geometries are compared with selected experimental results on brittle media. Accurate growth is demonstrated for domains discretised by refined and coarse volumetric meshes. Fracture and volume-based growth rates are shown to modify fracture interaction patterns. Two-dimensional cut-plane views of fracture networks show how fractures would appear on the surface of the studied volume.


[1] N. Thomas, A. Paluszny and R. W. Zimmerman. Growth of three-dimensional fractures, arrays, and networks in brittle rocks under tension and compression. Computers and Geotechnics, 2020. doi: 10.1016/j.compgeo.2020.103447

[2] Paluszny and R. W. Zimmerman. Numerical fracture growth modeling using smooth surface geometric deformation. Eng. Fract. Mech., 108, 19-36, 2013. doi: 10.1016/j.engfracmech.2013.04.012

How to cite: Paluszny, A., Thomas, R. N., Saceanu, M. C., and Zimmerman, R. W.: Thermo-mechanical Effects on Fracture Growth and Apertures in Three-Dimensional Subsurface Fractured Rocks , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14418,, 2021.

In this talk, I give an overview of our software ResFrac, which fully integrates a ‘true’ hydraulic fracturing simulator and a multiphase reservoir simulator (McClure et al., 2020a). Conventionally, these processes have been described with separate codes, using separate meshes, and with different physics. Integrating these two categories of software is advantageous because it enables seamless description of the entire lifecycle of a well. It is possible to seamlessly integrate wells with complex histories such as frac hits from offset wells, refracs, and huff and puff EOR injection.

ResFrac has been applied on 25+ studies for operators optimizing development of oil and gas resources in shale and has been commercially licensed by 15+ companies (;; The simulator has a modern user-interface with embedded help-documentation, wizards to help set up simulations, automated validators to identify issues with the setup prior to submitting, and plotting capabilities to preview 3D and tabular inputs. Simulations are run on the cloud and results are continuously downloaded to the user’s computer. This allows a user to easily run a large number of simultaneous simulations from their personal computer. The user-interface includes a custom-built and fully-featured visualization tool for 3D visualization and 2D plotting.

Hydraulic fracturing simulators must handle a diverse set of coupled physics: mechanics of crack propagation and stress shadowing, fluid flow in the fractures, leakoff, transport of fluid additives that impart non-Newtonian flow characteristics, and proppant transport. Proppant transport is particularly complex because proppant settles out into an immobile bed and may screen out at the tip. Many fracturing simulators approximate wellbore flow effects. However, because these effects are closely coupled to fracturing processes (especially in horizontal wells that have multiple simultaneously propagating fractures), we include a fully meshed, detailed wellbore model in the code, along with treatment of perforation pressure drop and near-wellbore tortuosity.

In the literature, separate constitutive relations are available to describe transport in open cracks, closed unpropped cracks, and closed propped cracks. However, there were not relations in the literature designed to describe transport under conditions transitional between these end-member states. A general numerical simulator must be able to describe all conditions (and avoid discontinuous changes between equations). To address this limitation, we developed a new set of constitutive equations that can smoothly transition between these end-member states – smoothly handling any general combination of aperture, effective normal stress, saturation, proppant volume fraction, and non-Newtonian fluid rheology (McClure et al., 2020).

The code solves all equations in a fully coupled way, using an adaptive implicit method. The fully coupled approach is chosen because of the tight coupling between many of the key physical processes. Iterative coupling converges very slowly and/or forces excessively small timesteps when tightly coupled processes are handled with iterative or explicit coupling.

McClure, Kang, Hewson, and Medam. 2020. ResFrac Technical Writeup (v5). arXiv.

How to cite: McClure, M.: Fully Integrating a Hydraulic Fracturing, Reservoir, and Wellbore Simulator into a Practical Engineering Tool, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9770,, 2021.

EGU21-7928 | vPICO presentations | EMRP1.18

A Geothermal Energy Concept based on Heat Storage in Geological Media

Rubén Vidal, Maarten W. Saaltink, and Sebastià Olivella

Aquifer Thermal Energy Storage (ATES) can help to balance energy demand and supply to make better use of infrastructures and resources. ATES consists of a pair or more wells that simultaneously inject or extract thermal energy into aquifers. The aim of ATES is to store the excess of energy during summer and to reuse it during winter, when there is an energy deficit. High-temperature Aquifer Thermal Energy Storage (HT-ATES) provides a good option to store water over 50°C, but it requires facing some problems, such as low efficiency recoveries and the uplift of the surface. Coupled thermo-hydro-mechanical (THM) modelling is a good tool to analyze the viability and cost effectiveness of the HT-ATES systems and understand the interaction of processes, such as heat flux, groundwater flow and ground deformation. We present the 3D THM modelling of a pilot HT-ATES system, inspired by one of the projects of HEATSTORE, which is a GEOTHERMICA ERA-NET co-funded project. The model aims to simulate the injection of hot water of 90°C in a central well and the extraction of water in four auxiliary wells during summer. In winter, the auxiliary wells inject water of 50°C and the central well extract water. The loading lasts longer than the unloading (8 months versus 4 months) and overall more heat is injected than extracted. We found that the system is more efficient in terms of energy recovery, the more years the system is operating. In the aquifer, both thermal loads and hydraulic loads have an important role in terms of displacements. At the surface, the vertical displacements are only a consequence of the hydraulic strains generated by the injection of water in the aquifer.

How to cite: Vidal, R., Saaltink, M. W., and Olivella, S.: A Geothermal Energy Concept based on Heat Storage in Geological Media, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-7928,, 2021.

EGU21-8561 | vPICO presentations | EMRP1.18

Flow-induced microfracturing of granite in superhot geothermal environments

Ryota Goto, Noriaki Watanabe, Kiyotoshi Sakaguchi, Youqing Chen, Takuya Ishibashi, Eko Pramudyo, Francesco Parisio, Keita Yoshioka, Kengo Nakamura, Takeshi Komai, and Noriyoshi Tsuchiya

Superhot geothermal environments with temperatures of approximately 400-500C at depth of approximately 2-4 km are expected as a new geothermal energy frontier. In order to efficiently exploit the superhot geothermal resources, fracture systems are necessary as flow path of working fluid. Hydraulic fracturing is a promising technique because it is able to create a new fracture system or enhance the permeability of preexisting fracture system. Laboratory-scale hydraulic fracturing experiments of granite have demonstrated the formation of densely distributed network of permeable fractures throughout the entire rock body at or near the supercritical temperature for water. Though the process has been presumed to involve continuous infiltration of low-viscosity water into preexisting microfractures followed by creation and merger of the subsequent fractures, plausible criterion for the fracturing is yet to be clarified. The possibility that the Griffith failure criterion is available to predict the occurrence of fracturing was shown by hydraulic fracturing experiments with acoustic emission measurements of granite at 400C under true triaxial stress. The present study provides a theoretical basis required to establish the procedure for hydraulic fracturing in superhot geothermal environment.

How to cite: Goto, R., Watanabe, N., Sakaguchi, K., Chen, Y., Ishibashi, T., Pramudyo, E., Parisio, F., Yoshioka, K., Nakamura, K., Komai, T., and Tsuchiya, N.: Flow-induced microfracturing of granite in superhot geothermal environments, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8561,, 2021.

Wellbore instability is one of the most serious drilling problems increasing well cost in well construction processes. It is widely known that many wellbore instability problems are reported in shale formations where water sensitive clay mineral exist. The problems become further complicated when the shale exhibits variation in strength properties along and across bedding planes. In this study, a coupled thermal-hydro-mechanical-chemical (THMC) model was developed for time-dependent anisotropic wellbore stability analysis considering chemical interactions between swelling shale and drilling fluids, thermal effects, and poro-elastoplastic stress-strain behaviors.

The THMC simulator developed in this work assumes that the shale formation behaves as an ion exchange membrane where swelling depends on chemical potential of drilling fluids invading from the wellbore to the pore spaces. The time-dependent chemical potential changes of water within the shale are evaluated using an analytical diffusion equation resulting in the evolution of swelling strain around the wellbore. On the other hand, the thermal and pressure diffusion equations are evaluated numerically by finite differences. The stress changes associated with thermal, hydro, and chemical effects are coupled to the 3D poroelastoplastic finite element model. The effects of bedding planes are also taken into account in the FEM model through the crack tensor method in which the normal and tangential stiffnesses of the bedding planes have stress dependency. The failure of the formation rock is judged based on the critical plastic strain limit.

The numerical analysis results indicate that the rock strength anisotropy induced by the existence of bedding planes is the most important factor influencing the stability of the wellbore among various THMC process parameters investigated in this work. The numerical results also reveal that an established theory to orient the wellbore in the direction of the minimum principal stress is not always a favorable option when the effect of the anisotropy of in-situ stresses and the distribution angle of bedding planes cancel out each other. Depending on both the distribution angle of bedding plane and ratio of the vertical to the horizontal stress, the trend of minimum mud pressure showed a great variation as predicted by the yield and failure criterion implemented in the model. Furthermore, the analysis results reveal that the distribution and evolution of plastic strains caused by the THMC processes have the time dependency, which can be controlled by the temperature and salinity of the drilling fluids.

The numerical wellbore stability analysis model considering shale swelling and bedding plane effects provides an effective tool for designing optimum well trajectories and determining safe mud weight windows for drilling complex shale formations. The time-dependent margins of safe mud weight window of drilling can be fine-tuned when the interaction among various parameters is fully considered as the THMC processes.

How to cite: Xue, J. and Furui, K.: Coupled Thermal-Hydro-Mechanical-Chemical Modeling for Time-Dependent Anisotropic Wellbore Stability Analysis, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9138,, 2021.

Numerous intersected rock fractures constitute the fracture network in enhanced geothermal systems. The complicated convective heat transfer behavior in intersected fractures is critical to the heat recovery in fractured geothermal reservoirs. A series of three-dimensional intersected fracture models are constructed to perform the flow-through heat transfer simulations. The geometries effects of dead-end fractures on the heat transfer are evaluated in terms of intersected angles, apertures, lengths, and the connectivity. The results indicate that annular streamlines appear in the rough dead-end fracture and cause an ellipsoidal distribution of the cold front. Compared to the steady flow in plate dead-end fractures, the fluid flow formed in the rough dead-end fracture enhances the heat transfer. Both the outlet water temperature Tout and heat production Q present the largest when the intersected angle is 90°. A larger intersected angle and longer length extension of the intersected dead-end fracture, raising Tout and Q, are beneficial to the heat production, while increasing the aperture is ineffective. Solely increasing numbers of dead-end fractures poses a little increase on Tout and Q. More significant heat extraction is obtained through connecting these dead-end fractures with the main flow fracture forming the flow network.

How to cite: chen, Y.: Convective heat transfer of water flowing in intersected rock fractures for enhanced geothermal extraction, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9347,, 2021.

EGU21-10538 | vPICO presentations | EMRP1.18

Characteristics and mechanism of CO2-based fracturing of granite in conventional and superhot geothermal environments

Eko Pramudyo, Noriaki Watanabe, Ryota Goto, Kiyotoshi Sakaguchi, Kengo Nakamura, and Takeshi Komai

Creation of fractured reservoir for enhanced geothermal system (EGS) in granitic rock at unconventional superhot geothermal environments (>400°C) has been found possible by injection of low viscosity supercritical water (Watanabe et al., 2017, Geoph. Res. Lett.; Watanabe et al., 2019, Sci. Rep.). Accordingly, the complex cloud-fracture network is formed through stimulation of pre-existing microfractures by the low viscosity (<50 μPa∙s) water. Nonetheless, water reactivity to rock forming minerals (Brown, 2000, Proc. 25th wksh. Resv. Engr.) and its high water footprint (Wilkins et al., 2016, Environ. Sci. Technol.) hinder water application. Therefore, CO2 is proposed to replace water, since it less reactive to rock forming minerals (Brown, 2000), and can reduce water footprint (Wilkins et al., 2016). CO2 has also low viscosity at conventional (c.a. 150 – 300°C) and superhot geothermal temperatures (based on Heidaryan et al., 2011, J. Supercrit. Fluid), which motivates to create similarly complex fracture pattern in those geothermal environments. A set of traditional-triaxial stress fracturing experiments was performed on cylindrical granite samples subjected to 200 – 450-°C temperature and varying differential stress to determine characteristics and mechanism of CO2-based fracturing in conventional and superhot geothermal environments. The experiments demonstrated that complex fracture pattern was formed at low pressure in all experimental-stress state and temperature. Breakdown pressure was also found to decrease with increasing differential stress. Hence, it was hypothesized that fracturing mechanism by injection of CO2 is governed by Griffith fracture theory. To unveil the fracturing process in detail, a CO2-based fracturing experiment was conducted on cubical granite sample with a realistic true-triaxial stress state and 300-°C temperature. The later experiment confirmed the above-mentioned hypothesis and showed that the fracturing occurs gradually.

How to cite: Pramudyo, E., Watanabe, N., Goto, R., Sakaguchi, K., Nakamura, K., and Komai, T.: Characteristics and mechanism of CO2-based fracturing of granite in conventional and superhot geothermal environments, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-10538,, 2021.

EGU21-9616 | vPICO presentations | EMRP1.18

On the consolidation behaviour of double-porosity geomaterials

Jinhyun Choo

Many natural and engineered geomaterials have double-porosity structure where two dominant pore systems coexist. Examples include structured soils where the two pore systems are inter-aggregate pores and intra-aggregate pores, and fissured rocks where the two pore systems are fissures and matrix pores. Although such double-porosity materials are frequently observed in geosciences and geoengineering applications, it remains mostly unclear how fluid flow and solid deformation interact differently in single- and double-porosity materials. The presentation explores this question through numerical simulation of consolidation – a paradigmatic problem in poromechanics – based on a recently developed modelling framework for fluid-infiltrated, inelastic materials with double porosity. Built on a combination of continuum principles of thermodynamics and standard plasticity theory, the framework can capture deformation, flow, and their coupling that occur individually in each pore system. Simulation results using this framework suggest that double-porosity structure gives rise to a two-staged consolidation behaviour, where the second stage appears similar to secondary compression in clays. It is also found that the simulated two-staged behaviour bears a striking semblance to experimentally observed consolidation processes in shales. These findings suggest that double-porosity structure may exert dominant control over the long-term hydro-mechanical behaviour of geomaterials.

How to cite: Choo, J.: On the consolidation behaviour of double-porosity geomaterials, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9616,, 2021.

This work presents the advancement of a self-developed Cellular Automata Software for engineering Rockmass fracturing processes (CASRock, in the applications of thermo-hydro-mechanical-chemical (THMC) processes in fractured geological media. It contains a series of previous developed numerical systems, namely EPCA for simulation of heterogeneous rock failure process, VEPCA for visco elastoplastic analysis, D-EPCA for rock dynamic response simulation, THMC-EPCA for coupled THMC processes in geological media and RDCA for simulation of rock cracking process from continuity to discontinuity. In CASRock, the non-isothermal, unsaturated fluid flow, mechanical process and chemical reaction are sequentially coupled by updating all the state variables using cellular automaton technique and finite difference method on spatial and temporal scale, respectively. The Lagrangian method is used to simulate the particle transport. The control equations, coupling scheme and numerical implementation are briefly introduced. Several applications, including  in the background of high level nuclear waste disposal are provided to show the abilities of CASRock in the simulation of coupling processes between physical fields. These applications include, (1) stability analysis of engineering rockmass under mechanical loading, (2) numerical study on coupled TM processes in hard rock pillar, (3) study on coupled THM processes in engineering barrier, (4) simulation of the THMC process in fractured rock.

How to cite: Pan, P.: CASRock—a code to simulate THMC processes in fractured geological media, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12203,, 2021.

EGU21-16373 | vPICO presentations | EMRP1.18

Thermodynamic-informed machine learning for pressure-sensitive plasticity model

WaiChing Sun and Nikolas Vlassis

This talk will present a machine learning framework that builds interpretable macroscopic surrogate elasto-plasticity models inferred from sub-scale direction numerical simulations (DNS) or experiments with limited data. To circumvent the lack of interpretability of the classical black-box neural network, we introduce a higher-order supervised machine learning technique that generates components of elasto-plastic models such as elasticity functional, yield function, hardening mechanisms, and plastic flow. The geometrical interpretation in the principal stress space allows us to use convexity and smoothness to ensure thermodynamic consistency. The speed function from the Hamilton-Jacobi equation is deduced from the DNS data to formulate hardening and non-associative plastic flow rules governed by the evolution of the low-dimensional descriptors. By incorporating a non-cooperative game that determines the necessary data to calibrate material models, the machine learning generated model is continuously tested, calibrated, and improved as new data guided by the adversarial agents are generated. A graph convolutional neural network is used to deduce low-dimensional descriptors that encodes the evolutional of particle topology under path-dependent deformation and are used to replace internal variables. The resultant constitutive laws can be used in a finite element solver or incorporated as a loss function for the physical-informed neural network run physical simulations.

How to cite: Sun, W. and Vlassis, N.: Thermodynamic-informed machine learning for pressure-sensitive plasticity model, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16373,, 2021.

EGU21-13883 | vPICO presentations | EMRP1.18

Wormhole Growth in Dissolving Limestones: Insights from 4D Tomography

Piotr Szymczak, Max P Cooper, Silvana Magni, Rishabh P Sharma, Tomasz P Blach, Andrzej P Radlinski, Marek Dohnalik, and Alessandro Tengattini

Dissolution of porous media introduces a positive feedback between fluid transport and chemical reactions at mineral surfaces leading to the formation of pronounced wormhole-like channels. While the impact of flow rate and reaction rate on the shapes of the wormholes is now well understood, much less is known about the dynamics of their propagation. In this study we capture the evolution of wormholes and their effects on flow patterns by in-situ X-ray microCT imaging of dissolving limestone cores. 4D tomography allows us in particular to correlate the permeability changes in a dissolving core with the advancement of the tip position of the wormhole. Surprisingly, we find that the relation between the two is highly nonlinear, with extensive periods of relatively fast growth of the wormhole which is nevertheless not reflected in any significant change in the overall permeability.  We hypothesize that this is caused by the presence of highly cemented regions in the core which act as permeability barriers for the flow. The presence of such regions is confirmed by a detailed analysis of the pore geometry based on the tomographic data. The results demonstrate that the analysis of the wormhole dynamics in 4d tomography can be used to probe the internal structure of the rock. 

How to cite: Szymczak, P., Cooper, M. P., Magni, S., Sharma, R. P., Blach, T. P., Radlinski, A. P., Dohnalik, M., and Tengattini, A.: Wormhole Growth in Dissolving Limestones: Insights from 4D Tomography, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-13883,, 2021.

EGU21-12226 | vPICO presentations | EMRP1.18

Hydro-geochemical processes in the emplacement cavern of a low and intermediate-level waste repository in an indurated clay-rock

Vanessa Montoya, Jaime Garibay-Rodriguez, and Olaf Kolditz

By 2080, Germany will have to store around 600 000 m3 of low and intermediate-level nuclear waste (L-ILW) with negligible heat generation. This kind of waste is largely made up of used parts of nuclear power stations such as pumps, pipelines, filters, etc. placed in various types of waste containers made from either steel, cast iron, or reinforced concrete in different designs and sizes (i.e. cylindrical or box shaped). It is already decided that a total of 303 000 of the 600 000 m3 L-ILW will be disposed in a final storage facility in the former iron ore mine Schacht Konrad which is under construction. However, it is still not clear where the L-ILW emplaced in in the old salt mine Asse (200 000 m3) will be stored in the future. The situation is particularly critical, as the waste have to be retrieved from the instable mine shafts partially flooded with groundwater, causing strong socio-political concerns as radioactive waste could contaminate the water nearby. For this reason, the new search for a nuclear waste repository for high-level waste (HLW), started in 2017, should also consider the possibility to accommodate the waste from Asse. Obviously, this is still subject to critics as this will make finding a final repository more difficult as storing HLW and L-ILW together requires different concepts and designs for each other and, above all, much more space.

In this context, in this contribution we have defined conceptual and numerical models to assess the hydro-chemical evolution of a L-ILW disposal cell in indurated clay rocks, involving the interaction of different components/materials and the expected hydraulic and/or chemical gradients over 100 000 years. The L-ILW disposal cell leverages a multi-barrier concept buried 400 m below the surface. The multi-barrier system is comprised of the waste matrix (i.e. backfilling the waste drums), the disposal container, the mortar backfill in the emplacement tunnel (where the disposal containers are located) and the clay host rock. The dimensions and design of the emplacement tunnel (e.g. 11 × 13 m) and disposal cells represent and consider some aspects taken into account in the designs of some European countries. In addition, tunnel walls reinforced with a shotcrete liner and the Excavation Damaged Zone is considered in the concept. The model is implemented in OpenGeoSys-6, an open-source version-controlled scientific software based on Finite Element Method which is capable of handling fully coupled hydro-chemical models by coupling OpenGeoSys to iPHREEQC. First calculation results, demonstrate that the most important processes affecting the near-field chemical evolution are i) the degradation of the concrete and cementitious grouts with porewater migrating inwards from the host rock and ii) the significant quantities of reactive and non-reactive gases (i.e. hydrogen, carbon dioxide and methane) that are generated as a result of: i) the anaerobic corrosion of metals present in the waste and containers and ii) the degradation of organic compounds by microbial and chemical processes. As a first approximation, some assumptions and simplifications have been considered, probably resulting in a wort case scenario.

How to cite: Montoya, V., Garibay-Rodriguez, J., and Kolditz, O.: Hydro-geochemical processes in the emplacement cavern of a low and intermediate-level waste repository in an indurated clay-rock, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12226,, 2021.

EGU21-13313 | vPICO presentations | EMRP1.18

Impact of dilatant strengthening and energy dissipation on fault slip stability

Antoine Jacquey, Manolis Veveakis, and Ruben Juanes

The temporal and spatial distribution of fluid pressure and temperature within a fault core are key determinants of the onset and nature (seismic or aseismic) of fault slip. Laboratory and field observations indicate that transient localization of fluid pressure and temperature often go hand in hand with strain localization upon seismic rupture: as slip occurs on a fault plane, temperature increases due to dissipated energy and fluid pressure decreases due to dilatant strengthening. An accurate description of this thermo-hydro-mechanical multiphysics coupling controlling slip mechanisms is therefore essential to characterize the stability of fault slip.

Here, we present results from analytical and numerical analyses of the stability of fault slip adopting a thermo-hydro-mechanical coupling scheme together with a rate-dependent plasticity formulation. In particular, we focus on the relevance of dilatant strengthening competing with energy dissipation as driving processes for stick-slip events and aseismic slip. We analyze the multiple steady states of the system and their respective stability by means of a numerical continuation technique, and we describe the dynamic evolution of deformation, fluid pressure and temperature fields by considering an associated transient problem.

The results presented here provide insights into the stability criterion for aseismic slip and the dynamic evolution of slip instability as a function of the physical (thermal and hydraulic) properties of the fault material and the boundary conditions (tectonic stresses and off-fault fluid pressure and temperature conditions). We identify two mechanisms for periodic slip, one driven by elastic loading and the other by multiphysics oscillations. We discuss the implications of these results for characterizing the transition from stable aseismic slip to unstable seismic slip in the context of natural and induced seismicity.

How to cite: Jacquey, A., Veveakis, M., and Juanes, R.: Impact of dilatant strengthening and energy dissipation on fault slip stability, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-13313,, 2021.

EGU21-12893 | vPICO presentations | EMRP1.18

CO2 fracturing using the phase field approach for the brittle fracture

Mostafa Mollaali, Vahid Ziaei-Rad, and Yongxing Shen

To simulate CO2 fracturing under an isothermal condition, we propose a phase field model. We take advantage of the ability of the phase field approach to predict fracture initiation and branching, as well as to avoid tracking the fracture path. We model the CO2 as a compressible fluid by modifying Darcy's law. In particular, we assume that the permeability is correlated to the value of the phase field by the exponential function. The dependence of the CO2 density as a function of the pressure is captured by the Span-Wagner state equation. The computed pressure breakdown values show good agreement with analytical solutions and experimental results.

How to cite: Mollaali, M., Ziaei-Rad, V., and Shen, Y.: CO2 fracturing using the phase field approach for the brittle fracture, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12893,, 2021.

EGU21-15088 | vPICO presentations | EMRP1.18

Hydro-Mechanical Modeling of the Year 2000 Hydraulic Stimulation of GPK2 Well, Soultz-sous-Forêts, France

Dariush Javani, Jean Schmittbuhl, and Francois Cornet

 Hydraulic stimulation of pre-existing fractures and faults plays a significant role in improving hydraulic conductivity of the fracture network around injection and production wells in geothermal reservoirs. It is therefore important to characterize the hydro-mechanical behavior of the faults against fluid injection. The Soultz-sous-Forêts site (France) has been an EGS pilot site where several major hydraulic stimulations have been performed and are well documented ( and

Here we use the 3DEC numerical modeling tool (Itasca) to analyze the year 2000 stimulation of GPK2 well where large scale seismic anomalies have been evidenced during the different stages of the stimulation using 4D-P-wave tomography (Calo et al, 2011). The specificity of the approach is to combine two modeling at different scales. First, a small-scale model (100x100x100 m3) is built to analyze the effective mechanical response of a stochastic discrete fracture network (DFN) following the statistical features of the observed fracture network (Massart et al, 2010). Second, a large-scale numerical model of the Soultz-sous-Forêts reservoir (5000x5000x5000 m3) containing the largest faults of the reservoir defined by Sausse et al., 2010, is developed including regional stresses. The objective is to constrain the large-scale mechanical properties of the surrounding matrix around the fault from the small-scale model, in particular, its hydro-mechanical behavior in terms of non-linear elastic response related to the stochastic DFN. As a first step only the largest fault (GPK3-FZ4770) is considered. The first stage of the stimulation is modelled as a constant flow rate of 30 ls-1 of water injected into the fault at the depth of approximately 4.7 km. We explored the effect of the normal and shear stiffness of the fault on the dynamical response of pore pressure along the fracture and the onset of slip. It is found that the increase of the aperture of the fault during the injection shows a slow migration (~2 cm/s) owing to poro-elastic effects. Also generated fluid pressure throughout the fault shows a long period oscillating behavior (~5 hr) sensitive to the magnitude of the fracture normal stiffness.

How to cite: Javani, D., Schmittbuhl, J., and Cornet, F.: Hydro-Mechanical Modeling of the Year 2000 Hydraulic Stimulation of GPK2 Well, Soultz-sous-Forêts, France, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15088,, 2021.

Understanding fluid flow in rough fractures is of high importance to large scale geologic processes and to most anthropogenic geo-energy activities. Here, we conducted fluid transport experiments on Carrara marble fractures with a novel customized surface topography. Transmissivity measurements were conducted under normal stresses from 20 to 50 MPa and shear stresses from 0 to 30 MPa. An open-source numerical procedure was developed to simulate normal contact and fluid flow through fractures with complex geometries. It was validated towards experiments. Using it, we isolated the effects of roughness parameters on fracture fluid flow. Under normal loading, we find that i) the transmissivity decreases with normal loading and is strongly dependent on fault surface geometry ii) the standard deviation of heights (hrms) and macroscopic wavelength of the surface asperities control fracture transmissivity. Transmissivity evolution is non-monotonic, with more than 4 orders of magnitude difference for small variations of macroscopic wavelength and roughness. Reversible elastic shear loading has little effect on transmissivity, it can increase or decrease depending on contact geometry and overall stress state on the fault. Irreversible shear displacement (up to 1 mm offset) slightly decreases transmissivity and its variation with irreversible shear displacements can be predicted numerically and geometrically at low normal stress only. Finally, irreversible changes in surface roughness (plasticity and wear) due to shear displacement result in a permanent decrease of transmissivity when decreasing differential stress. Generally, reduction of a carbonate fault’s effective stress increases its transmissivity while inducing small shear displacements doesn’t.

How to cite: Acosta, M., Maye, R., and Violay, M.: Effect of normal and shear loading on the hydraulic transport properties of calcite bearing faults with customized roughness., EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16294,, 2021.

EGU21-13536 | vPICO presentations | EMRP1.18

Selective mineral dissolution and permeability enhancement of fractured volcanic rocks by chelating agent flooding in geothermal environments 

Luis Salala, Noriaki Watanabe, Kaori Takahashi, Jose Erazo, and Noriyoshi Tsuchiya

Chemical stimulation using high-concentration hydrofluoric and hydrochloric acids has been a classic method to enhance the permeability of a geothermal reservoir. Our research group has recently proposed a new chemical stimulation using a weakly acidic (moderate-reactivity) aqueous solution containing an environmentally friendly chelating agent to create voids, which are sustained under crustal stress, by selective mineral dissolution with preventing precipitation by chelation of metal ions. In the present study, we have conducted chelating agent flooding experiments using an aqueous solution of pH 4 containing readily biodegradable chelating agent (GLDA) on various types of fractured volcanic rocks at 200 oC and effective confining stress of 15 MPa. The experiments have revealed fast permeability enhancement of up to approximately four times, from the initial value, in two hours. Further analyses have revealed phenocrysts of Fe-bearing minerals (ex. Hematite) dissolved faster than the groundmass of the rocks to create the voids. These results show the possibility of the new chemical stimulation.

Keywords: Chemical stimulation, Chelating agents, Geothermal energy, EGS

How to cite: Salala, L., Watanabe, N., Takahashi, K., Erazo, J., and Tsuchiya, N.: Selective mineral dissolution and permeability enhancement of fractured volcanic rocks by chelating agent flooding in geothermal environments , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-13536,, 2021.

EGU21-16079 | vPICO presentations | EMRP1.18

Dynamic permeability evolution during fluid injection and production in ultra-deep geothermal reservoirs

Samuel Scott, Alina Yapparova, Philipp Weis, and Matthew Houde

A geothermal well drilled into a reservoir at temperatures exceeding the critical point of pure water (>374 °C) could generate substantially greater quantities of energy than conventional geothermal wells. Although these temperatures can be found at shallow depths (<2-3 km) in high-grade geothermal resources located in volcanically active areas, similar temperatures are only found at depths >10 km beneath vast areas of continental crust with lower heat fluxes. Permeability decreases markedly with increasing depth below 2-3 km, so exploiting the tremendous heat resources of high temperature rock at such great depths will require permeability stimulation by the injection of high-pressure fluids. In this study, we use the CSMP++ platform to perform 3D simulations of transient permeability evolution around a geothermal doublet drilled to depths between 10-16 km. The simulations incorporate a well model initially devised by Peaceman (1978) to calculate well pressures and rates of fluid production/injection. The dynamic permeability model is based on Weis et al. (2012), initially developed to simulate the evolution of ore-forming magmatic-hydrothermal systems, and links a failure criterion for critically-stressed crust with depth-dependent permeability profiles characteristic for tectonically active crust as well as pressure- and temperature-dependent relationships describing hydraulic fracturing and the transition from brittle to ductile rock behavior. We investigate the permeability changes in response to high-pressure fluid injection in brittle and ductile rock, the timescales over which the zone of permeability stimulation migrates towards production wells, and dynamic permeability evolution in response to changes in injection and production parameters. These simulations aim to mitigate resource risks that could limit the ability to extract heat from geothermal resources in ductile upper crust and to help anticipate the conditions that would be required to make the exploitation of ultra-deep supercritical geothermal resources a reality. 


Peaceman, D. W. (1978) Interpretation of Well-Block Pressures in Numerical Reservoir Simulation. SPE 6893, 183–194.

Weis, P., Driesner, T., & Heinrich, C. A. (2012). Porphyry-copper ore shells form at stable pressure-temperature fronts within dynamic fluid plumes. Science, 338(6114), 1613–1616.

How to cite: Scott, S., Yapparova, A., Weis, P., and Houde, M.: Dynamic permeability evolution during fluid injection and production in ultra-deep geothermal reservoirs, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16079,, 2021.

EMRP2.10 – Space weather through ground and space magnetic data

EGU21-812 | vPICO presentations | EMRP2.10 | Highlight

Dynamics of ionospheric and telluric currents during large events of geomagnetically induced currents

Mirjam Kellinsalmi, Ari Viljanen, Liisa Juusola, and Sebastian Käki

Geomagnetic variations are mainly produced by external currents in the ionosphere and magnetosphere, and secondarily by induced (internal/telluric) currents in the conducting Earth. Large geomagnetically induced currents (GIC) are associated with large time derivatives of the horizontal magnetic field. Recent results show that the time derivative is typically dominated by the contribution from the telluric currents. Our study aims to find measures to quantify the behaviour of external and internal currents and their time derivatives during large GIC events. Results of this study show that strong external currents have quite narrow directional distributions. Angular variation is larger for internal currents, and especially for their time derivatives. For external currents angular variation is larger at higher latitudes. Similar behaviour is not seen with internal currents.

How to cite: Kellinsalmi, M., Viljanen, A., Juusola, L., and Käki, S.: Dynamics of ionospheric and telluric currents during large events of geomagnetically induced currents, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-812,, 2021.

EGU21-6934 | vPICO presentations | EMRP2.10 | Highlight

Sources of hazardous geoelectric fields in the United States during magnetic storms

Xueling Shi, Michael Hartinger, Joseph Baker, Paul Bedrosian, Benjamin Murphy, Anna Kelbert, Joshua Rigler, and Michael Ruohoniemi

Geomagnetic perturbations related to various phenomena in the near-Earth space environment can induce electric fields within the electrically conducting Earth. The geoelectric field is an important link between magnetospheric/ionospheric phenomena and geomagnetically induced currents in grounded electricity transmission networks. In evaluations of contiguous United Sates hazards, most previous studies have been focused on either 1-minute resolution geoelectric field measurements or geoelectric field time series derived from convoluting 1-minute geomagnetic field data with surface impedance tensors. To investigate sources of hazardous geoelectric fields during magnetic storms, including geoelectric fields induced by ultra-low frequency (ULF: 1 mHz to 1 Hz) waves, we use directly measured 1-second geoelectric field data from magnetotelluric survey stations that are distributed across the contiguous United States. Temporally-localized perturbations in measured geoelectric fields with a prominence of at least 0.5 V/km are detected during magnetic storms with a Dst minimum of at least -100 nT from 2008 to 2019. Most of these perturbations cannot have been resolved with 1-minute data since they correspond to phenomena that vary on smaller timescales and higher frequencies. The sources of geomagnetic perturbations inducing these extreme geoelectric fields can be categorized as interplanetary shocks, substorms, and ULF waves. We compare the geoelectric fields associated with the three sources and characterize their features. Extreme geoelectric fields related to these sources can have amplitudes of 1-2 V/km, comparable to the thresholds commonly used to identify hazardous events.

How to cite: Shi, X., Hartinger, M., Baker, J., Bedrosian, P., Murphy, B., Kelbert, A., Rigler, J., and Ruohoniemi, M.: Sources of hazardous geoelectric fields in the United States during magnetic storms, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-6934,, 2021.

EGU21-8410 | vPICO presentations | EMRP2.10

Investigation of the possibility of GIC development in Greece during the strongest magnetic storms of solar cycle 24 

Adamantia Zoe Boutsi, Georgios Balasis, Ioannis A. Daglis, Kanaris Tsinganos, and Omiros Giannakis

Geomagnetically Induced Currents (GIC) constitute an integral part of the space weather research and a subject of ever-growing attention for countries located in the low and middle latitudes. A series of recent studies highlights the importance of considering GIC risks for the Mediterranean region. Here, we exploit data from the HellENIc GeoMagnetic Array (ENIGMA), which is located in Greece, complemented by magnetic observatories in Italy, to calculate corresponding values of the GIC index, i.e., a proxy of the geoelectric field calculated entirely from geomagnetic field variations. We perform our analysis for the most intense magnetic storms (Dst<-150 nT) of solar cycle 24. Our results show a good correlation between the storm sudden commencement (SSC) and an increase of the GIC index value. These investigations indicate that despite the elevated amplitude of the GIC index the associated risk remains at low level for the power networks in Greece and Italy during the considered storm events.

How to cite: Boutsi, A. Z., Balasis, G., Daglis, I. A., Tsinganos, K., and Giannakis, O.: Investigation of the possibility of GIC development in Greece during the strongest magnetic storms of solar cycle 24 , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8410,, 2021.

EGU21-1083 | vPICO presentations | EMRP2.10

Spatio-temporal development of large scale auroral electrojet currents relative to substorm onsets

Sebastian Käki, Ari Viljanen, Liisa Juusola, and Kirsti Kauristie

The electric currents flowing in the ionosphere change rapidly and a large amount of energy is dissipated in the auroral ionosphere during auroral substorms. An important part of the auroral current system are the auroral electrojets whose profiles can be estimated from magnetic field measurements from low Earth orbit satellites. We have combined electrojet data derived from the Swarm satellite mission of ESA with the substorm database derived from the SuperMAG ground network data. We organize the electrojet data in relation to the location of the onset and obtain statistics for the development of the integrated current and latitudinal location for the auroral electrojets relative to the onset. Especially we show that just after the onset the latitudinal location of the maximum of the westward electrojet determined from Swarm satellite data is mostly located close to the onset latitude in the local time sector of the onset regardless of where the onset happens.

How to cite: Käki, S., Viljanen, A., Juusola, L., and Kauristie, K.: Spatio-temporal development of large scale auroral electrojet currents relative to substorm onsets, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-1083,, 2021.

EGU21-1263 | vPICO presentations | EMRP2.10

Introducing the “VALOR” analysis framework for auroral field-aligned current sheets using observations from Swarm Alpha and Charlie

Leonie Pick, Joachim Vogt, Adrian Blagau, and Nele Stachlys

The investigation of auroral field-aligned current (FAC) sheets is crucial in the context of space weather research since they serve as main transmitters of energy and momentum across geospace domains. Different magnetosphere-ionosphere coupling modes are reflected by the FACs’ multiscale nature with spatial scales, i.e., latitudinal extensions, ranging from below 1 km to hundreds of kilometers. The multiscale property can be addressed conveniently using ESA’s three-spacecraft mission Swarm. According to common practice a linear correlation analysis is performed on lagged and band-pass filtered scalar FAC density estimates from two nearby spacecraft.

We introduce the framework VALOR (Vectorial Association of Linearly Oriented Residua) which generalizes the common approach in two ways. First, VALOR utilizes the full magnetic field vector primarily observed at both spacecraft without filtering. Second, VALOR allows to test statistical association measures other than linear correlation in dependence of both time and along-track spacecraft lag. The method is further refined by considering the current sheet’s polarization, i.e., the directional preference of the associated magnetic field perturbation, which additionally constrains the sheet’s orientation.

Here, we apply VALOR to 1 Hz magnetic field observations from Swarm Alpha and Charlie and base the association measure on a vectorial version of the mean squared deviation. By means of a sample auroral oval crossing event we demonstrate that the incorporation of vectorial and polarization information helps to focus the association measure in the time-lag parameter plane leading to a smaller FAC spatial scale estimate. This result seems to hold in a statistical context including over 9000 quasi-perpendicular auroral oval crossings from 2014 to 2020. The fact that the VALOR derived FAC locations reflect the known ellipsoidal shapes of the auroral ovals speaks to the overall plausibility of the method as well as the independently supported finding that large-scale FACs (>300 km) dominate the dawn and dusk sectors while smaller scale FACs gain importance at noon and midnight. Among the various opportunities for future work are an application to 50 Hz high-resolution Swarm data as well as the investigation of the solar controlling parameters.

How to cite: Pick, L., Vogt, J., Blagau, A., and Stachlys, N.: Introducing the “VALOR” analysis framework for auroral field-aligned current sheets using observations from Swarm Alpha and Charlie, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-1263,, 2021.

EGU21-12667 | vPICO presentations | EMRP2.10

On Magnetic Helicity and Field-Aligned Currents in the Polar Ionosphere

Paola De Michelis, Giuseppe Consolini, Tommaso Alberti, Vincenzo Carbone, Roberta Tozzi, Igino Coco, and Fabio Giannattasio

Magnetic helicity, which is a measure of twist and linkage of magnetic field lines, is a useful quantity to investigate some processes occurring in space plasmas. In particular, there is a strong link between magnetic helicity, magnetic flux structures, turbulence and dissipation. We investigate the connection between the reduced magnetic helicity and the structure of field-aligned currents in the high-latitude ionosphere using high resolution (50 Hz) magnetic data collected on board the ESA Swarm constellation. We show the existence of a clear link between the multiscale coarse-grained structure of reduced magnetic helicity and the field-aligned currents. This finding strongly supports the idea that turbulence processes might be at the origin of the observed small-scale current structures. A discussion of the relevance of our results in the framework of the filamentary nature of the field-aligned current is also presented.

This work is supported by Italian PNRA under contract PNRA18_00289-A “Space weather in Polar Ionosphere: the Role of Turbulence ".

How to cite: De Michelis, P., Consolini, G., Alberti, T., Carbone, V., Tozzi, R., Coco, I., and Giannattasio, F.: On Magnetic Helicity and Field-Aligned Currents in the Polar Ionosphere, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12667,, 2021.

EGU21-2741 | vPICO presentations | EMRP2.10

High-latitude crochet (Solar flare effect) as a trigger of pseud-substorm onset

Masatoshi Yamauchi, Magnar Johnsen, Shin-Ichi Othani, and Dmitry Sormakov

Solar flares are known to enhance the ionospheric electron density and thus influence the electric currents in the D- and E-region.  The geomagnetic disturbance caused by this current system is called a "crochet" or "SFE (solar flare effect)".  Crochets are observed at dayside low-latitudes with a peak near the subsolar region ("subsolar crochet"), in the nightside high-latitude auroral region with a peak where the geomagnetic disturbance pre-exists during solar illumination ("auroral crochet"), and in the cusp ("cusp crochet").  In addition, we recently found a new type of crochet on the dayside ionospheric current at high latitudes (European sector 70-75 geographic latitude/67-72 geomagnetic latitude) independent from the other crochets.  The new crochet is much more intense and longer in duration than the subsolar crochet and is detected even in AU index for about half the >X2 flares despite the unfavorable latitudinal coverage of the AE stations (~65 geomagnetic latitude) to detect this new crochet (Yamauchi et al., 2020).  

The signature is sometime s seen in AL, causing the crochet signature convoluting with substorms.  From a theoretical viewpoint, X-flares that enhances the ionospheric conductivity may influence the substorm activity, like the auroral crochet.  To understand the substorm-crochet relation in the dayside, we examined SuperMAG data for cases when the onset of the substorm-like AL (SML) behavior coincides with the crochet.  We commonly found a large counter-clockwise ∆B vortex centered at 13-15 LT, causing an AU peak during late afternoon and an AL peak near noon at higher latitudes than the high-latitude crochet.  In addition, we could recognize a clockwise ∆B vortex in the prenoon sector, causing another poleward ∆B, but this signature is not as clear as the afternoon vortex.  With such strong vortex features, it becomes similar to substorms except for its local time.  In some cases, the vortex expends to the nightside sector, where and when nightside onset starts, suggesting triggering of onset.  Thus, the crochet may behave like pseudo-onset at different latitude than midnight substorms, and may even trigger substorm onset.

How to cite: Yamauchi, M., Johnsen, M., Othani, S.-I., and Sormakov, D.: High-latitude crochet (Solar flare effect) as a trigger of pseud-substorm onset, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-2741,, 2021.