A simple computer program for calculating stress and strain rate in 2D viscous inclusion-matrix systems

Computer-based numerical solutions of geomechanical problems are important to understand the processes forming rock structures as well as to quantify the associated pressure, stresses and strain rates. However, the development of such computer programs and the underlying mathematical methods are commonly not taught in a standard structural geology curriculum. Here, we present a simple computer program to calculate the stress, pressure, velocity and strain rate fields for two-dimensional (2D) viscous inclusion-matrix systems under pure shear and simple shear. The main aim of our contribution is to explain this computer program in a simple and transparent way, so that it can be used for introductory courses on geomechanical numerical modelling in structural geology. We present the governing equations of 2D viscous deformation and program the equations in the same order and style, so that the equations are still recognizable in the computer program. The computer program can treat stiff and weak inclusions of any shape and considers both linear and power-law viscous flow laws. We present numerical calculations for various inclusion-matrix scenarios. The program is written with the software MATLAB, is provided as supplementary material, and can also be run with the freely available software GNU Octave

Influence of cutting sequence and time effects on cutters and roof falls in underground coal mine -- numerical approach

Roof falls are among the most serious safety hazards faced by underground coal mines worldwide. Due to the stringent safety measures and development of the innovative support systems in the past few decades, their numbers were drastically reduced but have not been eliminated. Underground observations reveal that a number of larger roof falls are preceded by the development of shear failures near pillar ribs, termed cutters or guttering. In the past, many factors were identified as responsible for the development of cutters and ultimately roof falls. These factors can be broadly classified as stress related and non stress related. Although some useful work on the stress related aspects was conducted in the past, in this dissertation the cutter instability is investigated in more details while including some finer aspects of the mining process in particular the cutting sequence, which were not given due consideration before.;Three dimensional finite difference modeling has been carried out to accomplish the research objectives in this dissertation. The strain softening material behavior with cutting sequence has been used to realistically simulate the cutter formation as suggested by Gadde and Peng, 2005. A few cutting sequences employed by some U.S. coal mines have been considered in this study. This was done to understand if the cutting sequence has any significant influence over cutter formation. Apart from cutting sequence, factors such as the cut length, step cutting and the turning direction of crosscut into and away from major and minor horizontal stress are also examined for their effect on cutter development. Further, in contrast to past work, the effect of change in the immediate roof rock properties and horizontal stress directions are studied in several multiple entry models while simulating some realistic cutting sequences.;Field observations show that some cutters develop after a significant amount of time is elapsed since the area has been mined. While such time-dependent effects could be simulated with numerical modeling by using appropriate creep laws, due to the lack of knowledge on creep properties of coal measures rocks, realistic analysis is difficult at this stage. However, to consider the time effect on development of roof instability, a simple methodology has been suggested in this dissertation. The effect of different parameters like, entry and pillar width, intersection geometry has been correlated with the standup time for the observed roof falls at an IL Basin coal mine.;The combination of weak immediate roof and high horizontal stresses could have a devastating effect on roof stability. It\u27s very difficult to completely avoid roof instability in such conditions. In this research, based on the understanding of cutter development and roof falls, a few simple and practicable recommendations are made to minimize such instabilities. While these suggestions may not completely eliminate the roof failures, they may enhance the standup time to levels that will allow safer extraction of the reserve

Characterisation and modelling of natural fracture networks: geometry, geomechanics and fluid flow

Natural fractures are ubiquitous in crustal rocks and often dominate the bulk properties of geological formations. The development of numerical tools to model the geometry, geomechanics and fluid flow behaviour of natural fracture networks is a challenging issue which is relevant to many rock engineering applications. The thesis first presents a study of the statistics and tectonism of a multiscale fracture system in limestone, from which the complexity of natural fractures is illustrated with respect to hierarchical topologies and underlying mechanisms. To simulate the geomechanical behaviour of rock masses embedded with natural fractures, the finite-discrete element method (FEMDEM) is integrated with a joint constitutive model (JCM) to solve the solid mechanics problems of such intricate discontinuity systems explicitly represented by discrete fracture network (DFN) models. This computational formulation can calculate the stress/strain fields of the rock matrix, capture the mechanical interactions of discrete rock blocks, characterise the non-linear deformation of rough fractures and mimic the propagation of new cracks driven by stress concentrations. The developed simulation tool is used to derive the aperture distribution of various fracture networks under different geomechanical conditions, based on which the stress-dependent fluid flow is further analysed. A novel upscaling approach to fracture network models is developed to evaluate the scaling of the equivalent permeability of fractured rocks under in-situ stresses. The combined JCM-FEMDEM model is further applied to simulate the progressive rock mass failure around an underground excavation in a crystalline rock with pre-existing discontinuities. The scope of this thesis covers the scenarios of both two-dimensional (2D) and three-dimensional (3D) fracture networks with pre-existing natural fractures and stress-induced new cracks. The research findings demonstrate the importance of integrating explicit DFN representations and conducting geomechanical computations for more meaningful assessments of the hydromechanical behaviour of naturally fractured rocks.Open Acces

"Exhumation Mechanisms in Hot Orogens: insight from numerical modelling"

Hot orogens are characterised by high temperatures and large dimensions. The temperature is related to the crustal thickness achieved during the evolution of the orogenetic system. The high temperature achieved by this kind of orogens is high enough to induce a pervasive and massive partial melting of the middleupper crust. The partially melted rocks have a lower density and viscosity than the unmelted counterpart. The partially melted region that forms in the interior of the orogens affects the post collisional evolution of the orogen and affects the way in which the rocks exhume at the surface. The partially melted rocks self-organize in a hot channel that flows underneath the orogen, and affects large portions of the orogenetic system. Since this low viscosity channel has high mobility and it is highly sensitive to pressure changes, one of the most intriguing issues is to analyse if the changes of the topographic relief, due to the surface processes, affects the advection of the hot channel and therefore the exhumation patterns of the rocks. The best example of hot orogens is the Himalayan range in which there is abundant geological and geophysical evidence supporting widespread partial melting. The exhumation mechanism of the Great Himalayan Sequence (GHS) was one of the most important arguments of the Himalayan research in the last decades, and the models proposed have been based on the coupling between tectonics and climate. The most popular and predictive of these models was the channel flow model proposed by Beaumont [Beaumont et al., 2001] in which the GHS is interpreted as the coherent slice of Indian partially melted crust which has migrated from the Tibet, in the North, to the southern flank of the Himalayan range due to the established pressure gradient between the Tibetan plateu and the orogeny front. According to this model, the exhumation of the GHS is caused by the high and focused erosion in the souther flank. Recent researches have pointed out that the structural complexity of GHS is incompatible with the channel flow model. The aim of this thesis work is to reproduce all the principal features of hot orogens by using numerical modelling in order to i) test the influence of the surface processes on the partially melted crust and on the exhumation patterns of the rocks, and ii) compare the simulation results with the observations from the Himalayan range in order to verify if the exhumation patterns predicted by the channel flow model is a inherent feature of hot orogens. Results show that partial melting in the mid-upper crust always forms, owing to the accreted highly radiogenic sediments and the thick orogenic crust. Different patterns of exhumation and metamorphism are reproduced, most of which are consistent with different exhumation models proposed. Additionally, the results of the simulations give other important insights: i) the focused erosion is not the necessary condition to the exhumation of the partially melted rocks and of the channel flow; ii) the shape of the channel flow is affected by the timing of the erosion; iii) the sedimentation rate have first order implication on the development of the orogens

High-resolution imaging beneath the Santorini volcano

Volcanoes are surface expressions of much deeper magmatic systems, inaccessible to direct observation. Constraining the geometry and physical properties of these systems, in particular detecting high melt fraction (magma) reservoirs, is key to managing a volcanic hazard and understanding fundamental processes that lead to the formation of continents. Unfortunately, unambiguous evidence of magma reservoirs has not yet been provided due to the limited resolving power of the geophysical methods used so far. Here, a high-resolution imaging technique called full-waveform inversion was applied to study the magmatic system beneath the Santorini volcanic field, one of the most volcanically and seismically active regions of Europe. Quality-controlled inversion of 3d wide-angle, multi-azimuth ocean-bottom seismic data revealed a previously undetected high melt fraction reservoir 3 km beneath the Kolumbo volcano, a centre of microseismic and hydrothermal activity of the field. To enable the above method to handle land data, two major algorithmic improvements were added to the high-performance inversion code. First, to simulate instrument response of land seismometers, a pressure-velocity conversion has been implemented in a way that ensures reciprocity of the discretised 2nd-order acoustic wave equation. Second, the immersed-boundary method, originally developed for computational fluid dynamics, was implemented to simulate the wave-scattering off the irregular topography of the Santorini caldera. These advancements can be readily used to provide a higher-resolution image of the melt reservoir beneath the Santorini caldera already detected by means of travel-time tomography.Open Acces

Integrated stability mapping system for mines

The Integrated Stability Mapping System (ISMS) was developed as an engineering tool to quantify the geologic and geo-mechanical information of mines, and to integrate the critical stability influence factors into an overall stability index for use in mine planning and support design. It is generally understood that the inherent underground roof stability is determined by the interaction of both the given geologic characteristics and the local stress influences. Form this perspective, in this dissertation, the need for an integrated stability mapping system is established through investigating the traditional and current hazard mapping practices. In order to fulfill this need, computer aided hazard mapping techniques and popular numerical methods for geo-mechanical analysis are reviewed. Then, an integrated stability mapping system incorporating geology hazard mapping, geologic structural feature impacts, and advanced numerical stress analysis techniques into one solution has been developed.;The stability system is implemented inside the de-facto standard drawing environment, AutoCAD, and in compatible with widely used geology modeling software SurvCADD. This feature allows one to access numerous existing geologic data and mining information from present mine maps easily and directly. The LaModel stress calculation, a boundary element method, integrated within the mapping system can produce realistic and accurate stress and displacement analysis with its distinguished features such as the laminated overburden model, the true topography consideration and actual irregular pillar matching.;After the stability mapping system was developed, two case studies were performed to check for coding errors, calculation accuracy, and for demonstrating the functionalities and usefulness of the system. In the case studies, the composite stability index was compared with field observations. A good correlation has been found although only a few influence factors have been considered.;In the conclusion of this dissertation, it is suggested that the stability mapping system provides mining engineers with the ability to perform comprehensive, rapid and accurate multiple-factor stability mapping analysis. Then the resultant stability map can be a valuable guide to safer support designing and better mine planning, and ultimately increase the safety of mine design and reduce the injuries and fatalities associated with ground fall in underground mines

Wavefield Analysis of Rayleigh Waves for Near-Surface Shear-Wave Velocity

Shear (S)-wave velocity is a key property of near-surface materials and is the fundamental parameter for many environmental and engineering geophysical studies. Directly acquiring accurate S-wave velocities from a seismic shot gather is usually difficult due to the poor signal-to-noise ratio. The relationship between Rayleigh-wave phase velocity and frequency has been widely utilized to estimate the S-wave velocities in shallow layers using the multichannel analysis of surface waves (MASW) technique. Hence, Rayleigh wave is a main focus of most near-surface seismic studies. Conventional dispersion analysis of Rayleigh waves assumes that the earth is laterally homogeneous and the free surface is horizontally flat, which limits the application of surface-wave methods to only 1D earth models or very smooth 2D models. In this study I extend the analysis of Rayleigh waves to a 2D domain by employing the 2D full elastic wave equation so as to address the lateral heterogeneity problem. I first discuss the accurate simulation of Rayleigh waves through finite-difference method and the boundary absorbing problems in the numerical modeling with a high Poisson's ratio ( 0.4), which is a unique near-surface problem. Then I develop an improved vacuum formulation to generate accurate synthetic seismograms focusing on Rayleigh waves in presence of surface topography and internal discontinuities. With these solutions to forward modeling of Rayleigh waves, I evaluate the influence of surface topography to conventional dispersion analysis in 2D and 3D domains by numerical investigations. At last I examine the feasibility of inverting waveforms of Rayleigh waves for shallow S-wave velocities using a genetic algorithm. Results of the study show that Rayleigh waves can be accurately simulated in near surface using the improved vacuum formulation. Spurious reflections during the numerical modeling can be efficiently suppressed by the simplified multiaxial perfectly matched layers. The conventional MASW method can tolerate gentle topography changes with insignificant errors. Finally, many near-surface features with strong lateral heterogeneity such as dipping interfaces, faults, and tunnels can be imaged by the waveform inversion of Rayleigh waves for shallow S-wave velocities. This thesis consists of four papers that are either published (chapter 1) or in review (chapter 2, 3, and 4) for consideration of publication to peer-refereed journals. Each chapter represents a paper, and therefore inadvertently there will be a certain degree of overlap between chapters (particularly for the introduction parts, where references to many common papers occur)

Analogue and Numerical Models of Earthquake Rupture

Earthquakes represent one of the most important natural risks facing human populations in urban areas. Understanding the processes at the origin of these destructive events requires seismological observations, but also the use of laboratory analogues and numerical models for earthquake rupture. They allow for controlled conditions under which we can investigate the relative importance of different physical quantities involved in the system. The main points investigated in this thesis are the influence of loading rate on the nucleation of earthquakes, and the evolution of friction during dynamic ruptures. I conduct photoelastic experiments using polycarbonate plates, but also direct-shear experiments of precut granite blocks in a pressure vessel. I use finite-difference numerical models to reproduce and understand the dynamic laboratory ruptures, and I developed static finite element codes in order to reproduce the loading conditions induced by the experimental setup. The main results are that under certain conditions, increasing the loading rate makes the nucleation length shrink, and affects the nucleation position, which in this case is consistently situated on high coulomb stress areas. This is not necessarily the case for low loading rates. The shrinking of nucleation length may explain partly why some asperities in subduction zones can behave seismically or aseismically depending on the local tectonic loading velocity. Finally, I propose a method to estimate the dependence of friction on slip and slip velocity from strain gauge data during friction experiments. When conducted under realistic pressure conditions, this can provide useful constitutive laws to implement in numerical models simulating earthquakes. Eventually, the results presented in this thesis can be used in order to improve rupture scenarios, and short-term earthquake forecast