l'injection et de la circulation d'un fluide pouvant avoir un effet à long terme sont également étudiées séparément.

La sismicité induite est ensuite étudiée sous l'angle de l'utilisation de la loi de la conservation de l'énergie afin d'expliquer le mécanisme de génération d'ondes élastiques consécutif à la rupture. Le processus de conversion de l'énergie de déformation accumulée dans un système avant la rupture en énergie cinétique pendant la rupture, quelle que soit la vitesse avec laquelle la charge externe est appliquée au système, est clairement démontré. Cette approche permet de bien comprendre les facteurs qui ont un impact majeur sur l'intensité des ondes sismiques induites. Ensuite, une procédure de modélisation est proposée pour calculer les accélérations maximales induites à la surface du sol par le glissement de la faille. Ces accélérations sont utilisées pour évaluer l'intensité du séisme en relation avec la perception humaine des ondes sismiques et l'endommagement causé aux infrastructures, participant ainsi au processus d'acceptation par le public de tout nouveau projet d'exploitation géothermique.

1 Résumé La géothermie est une source d'énergie renouvelable pouvant servir de base pour la production d'électricité ou le chauffage, sans ou avec une émission limitée de gaz à effet de serre.

Cependant, le développement et l'exploitation des réservoirs géothermiques profonds s'accompagnent généralement d'une sismicité induite -un effet secondaire qui peut susciter de vives inquiétudes auprès du public et empêcher une utilisation à grande échelle de l'énergie géothermique. La pièce maîtresse d'une installation géothermique est un échangeur de chaleur, c'est-à-dire un réservoir, situé dans des roches chaudes en profondeur. Ce réservoir de roche est créé en utilisant la technique de la stimulation hydraulique, qui consiste essentiellement à injecter un fluide à des pressions élevées pour créer de nouvelles fractures hydrauliques qui se connecteront au réseau de fractures ou de failles préexistantes. Ce faisant de nouvelles voies pour la circulation du fluide seront créées et la surface d'échange de chaleur sera considérablement augmentée.

Plusieurs processus physiques couplés interviennent lors de la stimulation et de l'exploitation de réservoirs géologiques profonds. Il s'agit essentiellement de la déformation et de la rupture de la roche, de l'écoulement du fluide, et des variations de température. Ces processus couplés et dépendants du temps sont à l'origine de la propagation des fractures et de la réactivation de failles, responsable de la sismicité induite.

Cette thèse de doctorat est axée sur l'utilisation de simulations numériques dans le but de mieux comprendre le comportement thermo-hydro-mécanique des réservoirs géothermiques fracturés et de déterminer le potentiel de sismicité induite. La propagation des fractures hydrauliques et la réactivation de failles préexistantes sont modélisées à l'aide du concept de matériau cohésif en combinaison avec des éléments finis cohésifs et en tenant compte des processus physiques sous-jacents (par exemple, écoulement de fluide, déformation de la roche). Différents facteurs pouvant avoir un impact majeur sur la propagation de la fracture et la réactivation de la faille sont étudiés, tels que le frottement des failles préexistantes, le débit d'injection de fluide et l'orientation des failles. Les variations de température résultant de This doctoral research focuses on the use of numerical simulations to better understand the behavior of fractured geothermal reservoirs and the potential of induced seismicity. The propagation of hydraulic fractures and the reactivation of pre-existing faults are modeled using the cohesive material concept in combination with finite cohesive elements and taking into account the underlying physical processes (e.g., fluid flow, rock deformation). Various factors that may have a major impact on fracture propagation and fault reactivation are investigated, such as the friction condition of the pre-existing faults, the fluid injection rate, and the fault orientation. Temperature changes resulting from fluid injection and circulation that may have a long term effect are also separately investigated.

The induced seismicity is first studied from the standpoint of using the law of energy conservation in order to explain the mechanism of elastic waves generation due to failure. The process of converting the strain energy accumulated in a system prior to the failure into kinetic energy during failure, regardless of how slowly the external load is applied to the system, is clearly demonstrated. This approach provides a good understanding of the factors that have a major impact on the intensity of the induced seismic waves. A modeling procedure is then proposed to calculate the peak ground accelerations that are induced by the fault slip. The computed accelerations on the ground surface are used to assess the earthquake intensity, the human perception of the seismic waves, and the damage potential to infrastructures.

Chapter 1 Introduction

INTRODUCTION TO GEOTHERMAL ENERGY

Geothermal resources potential

Geothermal energy is the thermal energy contained in the Earth's crust. The sources of geothermal heat are associated with the upward conduction and convection of energy from the Earth's core and mantle, and with the decay of the naturally radioactive isotopes of uranium, thorium and potassium [START_REF] Barbier | Nature and technology of geothermal energy: A review[END_REF][START_REF] Boyle | Renewable energy[END_REF][START_REF] Brown | Mining the earth's heat: hot dry rock geothermal energy[END_REF]. These natural processes, which occur over geological time scales, maintain the renewal of the thermal energy content and make geothermal energy practically inexhaustible.

Geothermal energy exists in two major forms: natural hydrothermal systems and heat stored in hot rock formations. In natural hydrothermal systems, an interconnected network of open pores and fractures is present in the rock mass, which allows the flow of fluids. Natural groundwater is abundant in these systems; it is heated by the surrounding hot rock and is carried upward by natural convective circulation. These natural systems can be found in porous and/or naturally fractured rocks and are typically located at depths of 1 -4 km in a number of locations e.g., New Zealand, Indonesia, Northern Italy, Japan, Mexico, Iceland [START_REF] Dipippo | International Developments in Geothermal Power Production[END_REF][START_REF] Mock | Geothermal energy from the earth: its potential impact as an environmentally sustainable resource[END_REF][START_REF] Tester | Hot dry rock geothermal energy-a new energy agenda for the 21st century[END_REF]. However, natural hydrothermal systems are the exception rather than the rule as heat, natural groundwater and permeability/porosity must be present.

The heat stored in hot rock formations, on the other hand, is in principle available anywhere around the world if one drills deep enough and can appropriately stimulate the rock (i.e., create sufficient permeability) at these depths. The reason is that the temperature of the Earth increases with depth. On average, the geothermal gradient is about 30 °C/km [START_REF] Banks | An Introduction to Thermogeology: Ground Source Heating and Cooling[END_REF][START_REF] Barbier | Nature and technology of geothermal energy: A review[END_REF][START_REF] Brown | Mining the earth's heat: hot dry rock geothermal energy[END_REF]. However, there are many areas where thermal gradients are considerably greater. A thermal gradient of 50 °C/km is commonly available in the Western U.S [START_REF] Beckers | Levelized costs of electricity and direct-use heat from Enhanced Geothermal Systems[END_REF] where the minimum temperature at a depth of 5 km is estimated to be 250 ºC. The same thermal gradient was also found in Southern Australia (Geodynamics, 2015) where temperatures of 273 ºC and 283 ºC have been found at depths of around 5 km. A local geothermal gradient on the order of 65 °C/km was found at the Fenton Hill site in the Jemez Mountains of New Mexico, United States. For regions that have an average geothermal gradient of 30°C/km, a temperature of 150 °C is expected at depths of around 5 km, depths that are accessible with current drilling technology; this temperature is considered to be suitable for efficient heat exploitation [START_REF] Potter | Method of extracting heat from dry geothermal reservoirs[END_REF]. The rock formations at depths where suitable temperatures can be found have two typical characteristics: low permeability to almost impermeable and essentially drythe so-called hot dry rocks (HDR). HDRs hold an enormous amount of thermal energy that can be extracted for electricity production. A comparison of geothermal energy output with other energy sources is shown in Figure 1.1. [START_REF] Brown | Mining the earth's heat: hot dry rock geothermal energy[END_REF]. Note: 1 quad = 293 x 10 9 kWh; world energy consumption in 2010 is 524 quads It is obvious from Figure 1.1 that the geothermal energy stored in HDRs is second only to fusion energy, a source that is yet to be developed. It has been estimated that the total amount of thermal energy contained in HDR at accessible depths is about 300 times greater than the energy that oil, natural gas and coal combined could provide [START_REF] Duchane | Hot dry rock (HDR) geothermal energy research and development at Fenton Hill, New Mexico[END_REF][START_REF] Tester | Hot dry rock geothermal energy-a new energy agenda for the 21st century[END_REF]. This amount of energy is enough to supply all the energy needs of the world for centuries to come [START_REF] Edwards | Handbook of geothermal energy[END_REF].

Hot Dry Rock concept

The importance of HDR resources originates not just from its huge quantity but also from its abundant availability throughout the world. However, the low permeability and the absence of natural water in HDRs require that permeability-enhancing operations have to be carried out and a heat transport medium, e.g., water or brine, must be supplied artificially. The first attempts to exploit the heat contained in HDRs were made in the United States by a group of researchers from the Los Alamos National Laboratory [START_REF] Brown | Scientific progress on the Fenton Hill HDR project since[END_REF][START_REF] Smith | The Los Alamos scientific laboratory dry hot rock geothermal project (LASL Group Q-22)[END_REF], who conceived the HDR concept and began to develop the world's first HDR reservoir at Fenton Hill in New Mexico (USA) between 1974 and 1978 [START_REF] Murphy | Does hydraulic-fracturing theory work in jointed rock masses?[END_REF]. The original HDR concept, as shown in Figure 1.2, was simple [START_REF] Abé | Fundamentals of design concept and design methodology for artificial geothermal reservoir systems[END_REF][START_REF] Brown | Mining the earth's heat: hot dry rock geothermal energy[END_REF][START_REF] Potter | Method of extracting heat from dry geothermal reservoirs[END_REF]: A well is drilled into the HDR formation. Through this well, water is injected at pressures high enough to create a single plane vertical fracture, perpendicular to the direction of the minimum compressive stress, which is often horizontal at depths greater than 500 -1000 m. The vertical fracture would grow as a disc to a diameter of several hundred meters, creating a large reservoir of rock with sufficient area for heat exchange between the water and the hot rock. A second well is subsequently drilled to connect with the reservoir at some distance from the first well. Once the stimulation phase and the well completion are finished, water circulation and heat extraction can begin. Cold water is pumped into the reservoir through the first well. The water is heated by contact with the hot rock as it flows across the reservoir. The resulting superheated water is then pumped up to the surface via the second well and begins the power generation phase. After extracting the useful heat, the same water is re-injected back into the reservoir to begin a new circulation circle. Thus, an HDR system operates and recovers the heat from HDRs in a closed-loop circulation, making the production operations sustainable to some extent, although some water loss from the system may still occur. [START_REF] Brown | Mining the earth's heat: hot dry rock geothermal energy[END_REF].)

For large HDR systems, larger areas are required for heat exchange. This could be achieved by creating a set of multiple parallel fractures using hydraulic stimulation in open uncased sections of the borehole [START_REF] Nicol | Modelling the heat extraction from the Rosemanowes HDR reservoir[END_REF]. A schematic of such a system is illustrated in Figure 1.3.

The successful development of an HDR reservoir depends on the choice of a site that provides suitable rock at a depth that is reachable with limited drilling costs. The reservoir created should provide a sufficient area for rapid heat transfer between the rock and water and requires a reasonable amount of energy to pump water through it, i.e., resistance to flow must be relatively small. The reservoir must also operate with minimum water loss from the system and its temperature should not be drawn down too quickly.

S h S h

Figure 1.3 HDR concept for large geothermal systems that contains multiple parallel vertical fractures. Two deviated wells connect these fractures to create a closed-loop system (Adapted from [START_REF] Jung | EGS -Goodbye or Back to the Future[END_REF])

Enhanced Geothermal Systems

This original HDR concept was based on a consideration of the crystalline formation at great depths as an intact and almost impermeable rock mass. However, this consideration was refuted based on further understanding gained from the work at the Fenton Hill HDR project and the Camborne School of Mines (CSM) geothermal project at Rosemanowes, Cornwall, England [START_REF] Murphy | Does hydraulic-fracturing theory work in jointed rock masses?[END_REF]. The existence of natural joints and existing fractures has been confirmed and their influence on the development of the geothermal reservoir is considered to be significant. In fact, high fracture densities can be expected in most deep rock formations, although only a fraction of them are likely to be conductive [START_REF] Evans | Stress and rock mechanics issues of relevance to HDR/HWR engineered geothermal systems: Review of developments during the past 15 years[END_REF]. Microseismic evidence at the Fenton Hill project showed that instead of generating a disc-shaped fracture, hydraulic stimulation leads to the opening and/or shearing of the natural joints and existing fractures, and permanently improves the permeability of the formation [START_REF] Duchane | Hot dry rock (HDR) geothermal energy research and development at Fenton Hill, New Mexico[END_REF][START_REF] Evans | Stress and rock mechanics issues of relevance to HDR/HWR engineered geothermal systems: Review of developments during the past 15 years[END_REF][START_REF] Murphy | Does hydraulic-fracturing theory work in jointed rock masses?[END_REF][START_REF] Richards | The performance and characteristics of the experimental hot dry rock geothermal reservoir at Rosemanowes, Cornwall (1985-1988[END_REF]. The shearing and/or opening of existing fractures was also the main permeability enhancement mechanism observed in the CSM project [START_REF] Baria | Phase 2B Final Report of the Camborne School of Mines Project (Part 5 Sesmics. Section 1 Microseismic results)[END_REF][START_REF] Green | Phase 2B Final Report of the Camborne School of Mines Project (Part 5 Sesmics. Section 2 Crosshole seismic results)[END_REF]. In fact, the seismic investigations at the CSM project concluded that all detected microseismic activity indicated Injection well Production well Parallel fractures shear sources that are consistent with known rock joint orientations and in-situ stress magnitudes and orientations. These observations and further experimental studies prompted the concept of Enhanced Geothermal Systems (EGS) in which natural joints and existing fractures play important role in the development of the rock reservoir (Batchelor, 1989a;[START_REF] Cornet | Experimental investigations of forced fluid flow through a granite rock mass[END_REF]. Instead of being considered as a continuum medium, the rock is treated as a discrete medium that contains natural joints and fractures. The permeability of this medium could be enhanced by massive water injection in long uncased borehole sections to shear and open the natural fractures. The EGS concept is illustrated in Figure 1.4. Generally, an EGS reservoir must have a number of properties in order to be technically feasible and economically viable. These include the total volume, the total heat exchange surface, the flow impedance, and the thermal and stress-field properties. The reference limits for these properties are summarized in Table 1.1.

PIONEERING EGS PROJECTS

Since the creation of the world's first HDR reservoirs at Fenton Hill, many EGS projects have been developed in the United States, Europe and elsewhere. These projects have used a more or less similar approach to that used at Fenton Hill HDR reservoir, using the technique of hydraulic fracturing to create the large heat transfer areas necessary in a HDR reservoir. Major pioneering projects are summarized in the following paragraphs.

Fenton Hill HDR project

The first HDR geothermal project at Fenton Hill, New Mexico, United States, was initiated in 1971 by the Los Alamos National Laboratory (LANL). The primary objective of the project was to develop and demonstrate of an economical and commercially viable technology for extracting thermal energy stored in HDR. The site, located just outside of the ring fault structure of the Valles Caldera, was chosen primarily on the basis of favorable heat flow and the simplicity of the geologic structure in the anticipated reservoir rock [START_REF] Brown | Scientific progress on the Fenton Hill HDR project since[END_REF]. The estimated thermal gradient is on the order of 65 °C/km. The basement rock is at a depth of about 730 m, and contains Precambrian igneous and metamorphic rocks. The rock mass is extremely tight with hydraulic conductivity in the range of nano-darcies [START_REF] Brown | Mining the earth's heat: hot dry rock geothermal energy[END_REF]. A geologic cross section of the rock mass with the geothermal gradient is shown in ) and Phase II (EE-2 and EE-3) wellbore at the Fenton Hill HDR project [START_REF] Grigsby | Rock-water interactions in the Fenton Hill, new Mexico, hot dry rock geothermal systems I. fluid mixing and chemical geothermometry[END_REF] Two separate, confined HDR reservoirs were created at the Fenton Hill site in hot crystalline rock at different depths and then assessed and flow-tested for almost a year. Work on the shallower reservoir (Phase I reservoir) began in 1974 [START_REF] Brown | Mining the earth's heat: hot dry rock geothermal energy[END_REF]. The first borehole, GT-2 (Figure 1.5), was drilled into the granitic rock at a depth of 2900 m at a mean temperature of 195 °C. After a series of water injections to generate a hydraulic fracture from the bottom of the GT-2 borehole, the second borehole, EE-1, was directionally drilled beneath the bottom of the borehole GT-2 to intersect with the created hydraulic fracture. A good connection between the second borehole and the created reservoir was not immediately achieved, and sidetracking was necessary to establish contact with the initial well via a combination of induced and natural fracture pathways. In total, over 2000 m 3 of water was injected at average pressures of 10 MPa. Later, the reservoir was enlarged due to extensive testing. Using microseismic data, the final volume of the reservoir was estimated at 10 x 10 6 m 3 [START_REF] Brown | Hot dry rock geothermal energy: important lessons from Fenton[END_REF]. Circulation tests showed a thermal production power of 3 MW at a flow rate of about 6 L/s.

Development of the second, much larger, deeper, and hotter HDR reservoir began in 1980

(Phase II). The initial intent of this reservoir was to create a series of vertical, disc-shaped fractures connecting two inclined, vertically separated wellbores [START_REF] Brown | Scientific progress on the Fenton Hill HDR project since[END_REF].

This concept is presented in Figure 1.3. The two boreholes were drilled sequentially, one above the other. The first borehole EE-2 (Figure 1.5) was drilled to a vertical depth of about 4400 m, with the lower 1000 m directionally drilled at an angle of 35° to the vertical. The second borehole EE-3 (Figure 1.5) was then drilled in a similar manner to the first, but with the directionally drilled segment positioned 380 m above the directionally drilled portion of the first borehole. A number of hydraulic stimulations were carried out in an effort to connect these two boreholes by multiple vertical fractures. In total, over 21000 m 3 of water was injected at average pressures of 48 MPa. However, all these stimulations failed to connect the two boreholes. A satisfactory connection between the boreholes was then achieved after sidetracking the upper borehole into the stimulated region, which had been identified using induced microseismic data. The Phase II reservoir operated at a much higher temperature, above 300 °C, than the Phase I reservoir. A circulation test revealed a thermal power production of 10 MW at a circulation flow rate of 12 -14 L/s. Water loss during circulation was estimated to be extremely low due to the inherent tightness of the surrounding sealed rock mass, combined with the stress cage effect [START_REF] Brown | Mining the earth's heat: hot dry rock geothermal energy[END_REF]. However, the resistance to flow was 2.1 MPa/(L/s), which is unacceptably high compared to the reference limit value given in Table 1.1. It was believed that much of the impedance was concentrated in the zones immediately surrounding the boreholes. The Fenton Hill project was terminated in 1995 due to lack of financial support.

Camborne School of Mines Project

The Camborne School of Mines project at Rosemanowes Quarry in Cornwall, England started in 1977 and proceeded in 3 phases (Batchelor, 1989a;[START_REF] Parker | Overview, Hot Dry Rock geothermal energy, Phase 2B final report of the Camborne School of Mines Project[END_REF]. The Phase 1 of the project involved preliminary experiments at shallow depths. The objective of this phase was to understand the role of the natural fracture system and find the best way to establish and control flow paths within the fracture system. Four boreholes were drilled to a depth of 300 m;

explosives and hydraulic fracturing were used in combination to connect the four boreholes.

The circulation test later showed that the impedance was greatly reduced. The lowest impedance achieved was close to 0.1 MPa per L/s. Temperature gradients at the site were estimated at 30 -40 °C/km, the highest in England [START_REF] Parker | Overview, Hot Dry Rock geothermal energy, Phase 2B final report of the Camborne School of Mines Project[END_REF].

The next phase of the project, Phase 2A, began in 1980 and lasted until 1983. The goal of this phase was to learn how to create and manage a full-size reservoir that has all the parameters necessary for commercial operation, except the depth of the reservoir. During this phase, two boreholes (RH11 and RH12) deviated from the vertical by 30° and spaced by 300 m were drilled to a depth of 2000 m. The reservoir was created by explosive fracturing near the wellbore, followed by hydraulic stimulation of the lower borehole using water injected at flow rates up to 100 L/s to try to open up the existing joints to connect to the upper borehole [START_REF] Parker | Overview, Hot Dry Rock geothermal energy, Phase 2B final report of the Camborne School of Mines Project[END_REF]. However, the circulation test, in which water was injected through the lower borehole revealed, that a massive amount of water was lost from the system and the pumping pressure required was too high [START_REF] Baria | Phase 2B Final Report of the Camborne School of Mines Project (Part 5 Sesmics. Section 1 Microseismic results)[END_REF]. Seismic monitoring showed that during the stimulation and subsequent circulation, the reservoir grew vertically downwards, as shown in Figure 1.6, probably extending to a depth of two kilometers below the two boreholes.

During Phase 2B, a third borehole RH15 was drilled to a depth of 2600 m to access the stimulated zone below the existing two boreholes (Batchelor, 1989b). Its position was selected to penetrate through the microseismic cloud. Stimulation was then performed on RH15 using viscous gel. A substantial thermal drawdown was observed during the first year of a circulation test. Circulation tests later suggested that the optimum hydraulic performance at Rosemanowes would have an injection flow rate of 24 L/s, a flow impedance of 0.6 MPa/(L/s) and water loss of 21%. This is still far from the requirements, as given in Table 1.1, for this EGS reservoir to be considered to be economically viable. The Soultz granitic basement lies below 1.4 km of sediments, goes down to 5 km, and is strongly fractured [START_REF] Genter | Analysis of macroscopic fractures in granite in the HDR geothermal well EPS-1, Soultz-sous-Forêts, France[END_REF][START_REF] Ledésert | Fractal analysis of fractures applied to Soultz-sous-Forets hot dry rock geothermal program[END_REF][START_REF] Sausse | Hydromechanical properties and alteration of natural fracture surfaces in the Soultz granite (Bas-Rhin, France)[END_REF]. The joint systems are clustered, with a high density of joints in the fracture zones and a much lower density in the competent rock. The fracture zones and joints are mainly subvertical at a strike of 160°. The geostatic stress state includes a minimum horizontal stress S h , which is about 54% of the vertical stress S v (S h = 0.54 S v ), and a maximum horizontal stress S H , which is almost equal to the vertical stress (S H = S v ) [START_REF] Cornet | How close to failure is a granite rock mass at a 5 km depth?[END_REF][START_REF] Ledésert | Fractures, hydrothermal alterations and permeability in the Soultz Enhanced Geothermal System[END_REF][START_REF] Sausse | Types of permeable fractures in granite[END_REF]. Temperature anomalies have been documented at the site and in the region around Soultz [START_REF] Dèzes | Evolution of the European Cenozoic Rift System: interaction of the Alpine and Pyrenean orogens with their foreland lithosphere[END_REF][START_REF] Ziegler | European Cenozoic rift system[END_REF], which indicate that some of the faults are permeable, transporting water from great depth up to the cap rock. The temperature at 5 km in the granitic basement is about 200 °C [START_REF] Freymark | The deep thermal field of the Upper Rhine Graben[END_REF][START_REF] Held | Economic evaluation of geothermal reservoir performance through modeling the complexity of the operating EGS in Soultz-sous-Forêts[END_REF]Hettkamp et al., 2004) and the lower part of the temperature profile indicates a conductive heat transfer regime.

In total, five deep boreholes have been drilled at the Soultz site [START_REF] Hébert | The Enhanced Geothermal System of Soultz-sous-Forêts: A study of the relationships between fracture zones and calcite content[END_REF]. L/s. The maximum wellhead pressure produced was 17.9 MPa. The borehole GPK4 was stimulated by injecting 21500 m 3 of water at a flow rate from 30 to 45 L/s and the peak wellhead pressure was 17.5 MPa [START_REF] Evans | A survey of the induced seismic responses to fluid injection in geothermal and CO2 reservoirs in Europe[END_REF]. The permeability of the rock mass was significantly improved after stimulation. Transmissibility of single fractures after the stimulation could exceed 1 d.m and was only slightly pressure dependent [START_REF] Schill | Hydraulic performance history at the Soultz EGS reservoirs from stimulation and long-term circulation tests[END_REF].

Several circulation tests were later performed. The production flow rate for the two production boreholes GPK2 and GPK4 reached 30 L/s, which is still far from the requirement of 50 -100 L/s (Table 1.1). Long term circulation of the triplet system with power production commenced in 2010 [START_REF] Genter | The Soultz geothermal adventure: 20 years of research and exploration of deep crystalline fractured rocks for EGS development[END_REF].

(a) (b)

INDUCED SEISMICITY IN EGS DEVELOPMENT

The centerpiece of an EGS project is the heat exchanger created in the hot rock at depth.

Hydraulic stimulation has become a common technique to achieve this result. Fluid, such as water, is injected at sufficiently high pressure into the rock mass to create new fractures or activate the pre-existing ones. Given the high pressures and the massive amount of water needed in the stimulation, the seismic risk has increasingly become a major concern [START_REF] Evans | A survey of the induced seismic responses to fluid injection in geothermal and CO2 reservoirs in Europe[END_REF][START_REF] Zang | Analysis of induced seismicity in geothermal reservoirs -An overview[END_REF]. In fact, one of the main challenges for the development of an EGS is associated with assessing the induced seismicity that may be triggered by hydraulic pressurization during the stimulation and circulation phases [START_REF] Ellsworth | Injection-Induced Earthquakes[END_REF][START_REF] Pellet | Rock mechanics is meeting the challenges of geo-energies[END_REF][START_REF] Wohlenberg | Monitoring and interpretation of seismic observations in hot dry rock geothermal energy systems[END_REF]. Thousands of seismic events are generated during the stimulation and exploitation of EGS reservoirs. In some sites, events of high magnitude were induced and could be felt on the ground surface.

At Soultz-sous-Forêts, massive water injections into the boreholes GPK2, GPK3 and GPK4

(located at depths from 4400 to 5000 m) resulted in approximately 90,000 micro-earthquakes.

The stimulation of the borehole GPK2 alone triggered over 700 events of magnitude between M L 1.0 and 2.5 during injection. A large event of magnitude M L 2.6 occurred after the shut-in [START_REF] Dorbath | Seismic response of the fractured and faulted granite of Soultz-sous-Forêts (France) to 5 km deep massive water injections[END_REF][START_REF] Evans | A survey of the induced seismic responses to fluid injection in geothermal and CO2 reservoirs in Europe[END_REF]. More than 200 event of magnitude between M L 1.0 and 2.5 occurred during the stimulation of borehole GPK3; the largest event that had a magnitude M L 2.9 also occurred after the shut-in. Stimulation of the borehole GPK4 also resulted in multiple seismic events, some reaching a magnitude M L 2.7; in total, more than 30 events of magnitude above M L 2.0 were triggered [START_REF] Rybach | The future of geothermal energy" and its challenges[END_REF]. The circulation tests performed in 2005 and 2008 also induced seismicity with the largest event having a magnitude of M L 2.3 [START_REF] Cuenot | Analysis of the Microseismicity Induced by Fluid Injections at the EGS Site of Soultz-sous-Forêts (Alsace, France): Implications for the Characterization of the Geothermal Reservoir Properties[END_REF][START_REF] Dorbath | Seismic response of the fractured and faulted granite of Soultz-sous-Forêts (France) to 5 km deep massive water injections[END_REF].

At the Basel EGS site in Switzerland, the stimulation of a 5000 m deep borehole in 2006 triggered multiple seismic events that were clearly perceived by the public. The first event of magnitude M L 2.6 occurred during fluid injection, followed by two larger events of magnitude M L 2.7 and 3.4 that occurred during shut-in. After the bleed-off, three aftershocks of M L > 3 were recorded [START_REF] Deichmann | Earthquakes induced by the stimulation of an enhanced geothermal system below Basel (Switzerland)[END_REF][START_REF] Häring | Characterisation of the Basel 1 enhanced geothermal system[END_REF]. These induced seismic events raised concerns among the public, which led to the cancelation of the project by the Basel government at the end of 2009 and the borehole was closed in 2011.

Induced seismicity has also occurred during the development phase of many other EGS projects. For instance, multiple induced seismic events of magnitude approaching M L 3.7 were recorded at the Cooper Basin EGS project in Australia after the injection of 20,000 m 3 of water [START_REF] Asanuma | Microseismic monitoring of hydraulic stimulation at the Australian HDR project in Cooper Basin[END_REF][START_REF] Baisch | Induced seismicity during the stimulation of a geothermal HFR reservoir in the Cooper Basin, Australia[END_REF]. As recently as 2017, an earthquake of magnitude M W 5.5 that struck South Korea in November 2017 was believed to be linked to the stimulation activities at the nearby Pohang EGS project [START_REF] Grigoli | The November 2017 Mw 5.5 Pohang earthquake: A possible case of induced seismicity in South Korea[END_REF].

Many induced microseismic events have also been reported during the production phases of various EGS projects. At the Hellisheidi geothermal power plant, which is located within the Hengill volcanic system in Iceland, microseismic activity increased immediately after reinjection of geothermal waste water began [START_REF] Flóvenz | Reinjection and induced seismicity in geothermal fields in Iceland[END_REF][START_REF] Juncu | Injection-induced surface deformation and seismicity at the Hellisheidi geothermal field, Iceland[END_REF], with two largest induced seismic events of magnitude 4. Later investigations showed that the fluid reinjection caused an increase in pore pressure which led to increased seismicity and fault slip [START_REF] Juncu | Injection-induced surface deformation and seismicity at the Hellisheidi geothermal field, Iceland[END_REF]. Both the pore pressure increase and the fault slip resulted in the observed surface deformation. An induced event of magnitude M L = 2.7 was also recorded in August 2009 during the circulation phase at the geothermal project in Landau, Germany [START_REF] Gaucher | Induced seismicity in geothermal reservoirs: A review of forecasting approaches[END_REF].

The risk of induced seismicity has to be thoroughly assessed and quantified before an EGS project can begin, especially for projects located near urban areas. Otherwise, this issue can impede the development or even lead to the cancelation of the project as was seen in Basel, Switzerland. Seismicity induced during the injection and withdrawal of fluids is believed to be associated with the slip of pre-existing faults and the propagation of new tensile fractures [START_REF] Pearson | The relationship between microseismicity and high pore pressures during hydraulic stimulation experiments in low permeability granitic rocks[END_REF]. Both the fault slip and the fracture propagation are related to the in-situ stress variation, which is induced by thermo-poromechanical processes involved during the injection and withdrawal of fluids [START_REF] Ghassemi | Thermo-poroelastic effects on reservoir seismicity and permeability change[END_REF][START_REF] Ghassemi | Analysis of Fracture Propagation under Thermal Stress in Geothermal Reservoirs[END_REF][START_REF] Ghassemi | Poro-thermoelastic mechanisms in wellbore stability and reservoir stimulation[END_REF][START_REF] Segall | A note on induced stress changes in hydrocarbon and geothermal reservoirs[END_REF]. The thermo-mechanical coupling (i.e., changes in in-situ stress due to temperature variations) is important on long time scales from months to years [START_REF] Ghassemi | A Review of Some Rock Mechanics Issues in Geothermal Reservoir Development[END_REF][START_REF] Ghassemi | Analysis of Fracture Propagation under Thermal Stress in Geothermal Reservoirs[END_REF]. While the poro-mechanical effect due to changes in pore pressure affects the deformation of the rock mass more quickly.

Consequently, seismic events induced during the fluid injection are usually considered to be caused by the shear slip on natural faults in response to a reduction of the shear strength the faults due to an increase in pore pressure. However, thermoelastic perturbations can cause an increase of the deviatoric stress in the cooled portion of the reservoir, potentially provoking an unstable condition [START_REF] De Simone | Thermal coupling may control mechanical stability of geothermal reservoirs during cold water injection[END_REF]. Thus, a comprehensive approach to assess the risk of induced seismicity should take into consideration the effects of both thermo-and poromechanical processes. A schematic of the fault slip due to fluid injection, the propagation of the induced seismic waves and the potential impact on surface structures is illustrated in Figure 1.8. The intensity of induced seismicity can be estimated based on the peak accelerations or peak velocities recorded on the ground surface and using the modified Mercalli intensity scale [START_REF] Wald | Relationships between peak ground acceleration, peak ground velocity, and modified Mercalli intensity in California[END_REF][START_REF] Wood | Modified Mercalli intensity scale of 1931[END_REF]. Figure 1.2 provides a relationship between the peak ground accelerations or peak ground velocities and the seismic intensity, human perception and potential structural damage. The consequences of the induced seismicity can range from a nuisance to the population to structural damages to the infrastructure.

SUMMARY

This chapter first describes the central problem of this doctoral research which is the propagation of hydraulic fractures and their interactions with pre-existing faults during the development of geothermal reservoirs. The chapter then provides a summary of the theoretical basis for the physical processes that are involved. The equations of the theory of fluidsaturated thermo-poroelastic media are presented first. These include equations for deformation of porous media, the fluid flow in the pores and in the fractures and faults, and the heat transfer. Next, the criteria for fracture propagation, based on fracture mechanics and the concept of cohesive material, are presented. The dynamic aspect of fracture propagation is also included. The phenomenon of fault reactivation and the criterion for fault reactivation are then reviewed with an emphasis on using the Coulomb friction criterion for assessing fault reactivation. The equations for single hydraulic fracture models are also presented. Those models serve as reference solutions for the validation of the numerical models. Finally, the numerical methods for modeling fracture propagation are presented, which include the extended finite element method and the cohesive element technique.

INTRODUCTION

This doctoral research is centered on the numerical modeling of the propagation of hydraulic fractures and the potential reactivation of pre-existing faults during hydraulic stimulation of deep geothermal reservoirs. The process of hydraulic stimulation is schematically presented in Figure 2.1. A borehole is drilled into a hot dry rock formation and fluid is injected at pressures high enough to create new hydraulic fractures (Figure 2.1a). As the hydraulic fractures propagate they may intersect pre-existing faults (Figure 2.1b). The fluid then flows into these faults, leading to an increase in the fluid pressure in the faults, which in turn causes the faults to slip. The fault slip is both beneficial and hazardous. On the one hand, the fault slip will improve the fault permeability due to dilatancy effects and fault roughness, thus enhancing the overall permeability of the reservoir which is beneficial. On the other hand, the fault slip may also be unstable, which means that the fault slip may induce seismic waves; if the induced seismic waves are strong enough, they can cause harm to people and damage structures on the ground, which is hazardous. The propagation of hydraulic fractures or the reactivation of pre-existing faults in a permeable porous rock mass that is fully saturated due to the injection of fluid is a coupled problem that involves multiple physical processes:

(i) deformation of the porous rock mass;

(ii) flow of the fluid within the pores;

(iii) flow of the fluid within the hydraulic fractures and the pre-existing faults;

(iv) propagation of the hydraulic fractures or the reactivation of the pre-existing faults;

(v) heat transfer.

In-situ conditions such as the presence of natural fractures, incomplete data regarding the initial conditions (e.g., stresses and pore pressure) or physical properties of the rock mass, may add additional complexities to the problem. Furthermore, the propagation of brittle fractures is inherently a dynamic process that releases energy in the form of elastic waves. The effects of these waves may be significant over a short time period (on the order of seconds) following the sudden propagation of the fracture. Meanwhile, the diffusion of fluid within the pores and the fluid injection operation itself occur over much longer time period (on the order of hours or even days). Thus, coupling the inertial effects with the fluid diffusion is an unrealistic task.

A modeling procedure to quantify the induced dynamic effects is presented in Chapter 5.

Thermal effects also need to be accounted for since hydraulic stimulation of EGS involves temperature changes of significant orders as cold water is pumped into hot rock formations.

The above phenomena and the equations for describing these phenomena are presented in the next sections.

GENERAL GOVERNING EQUATIONS

Deformation of porous rock mass

The rock mass is macroscopically considered as a homogeneous, isotropic, and poroelastic medium whose constitutive relation is expressed as [START_REF] Coussy | Poromechanics[END_REF]:

    ,0 0 0 2 23 3 ij ij kk ij ij ij ij K G G K T T b p p                    (2.1)
where ij  is the current total stress (compression is negative), ,0 ij  is the initial stress; T and 0 T are the current and initial temperatures, respectively; K and G are the bulk modulus and the shear modulus of the skeleton, respectively; α is the linear thermal expansion coefficient of the skeleton; and b is the Biot's coefficient, which is related to the bulk modulus K of the skeleton and the bulk modulus K S of the solid phase as:

1 S K b K  (2.2)
The bulk modulus K and the shear modulus G are related to the Young's modulus E and the Poisson's ratio ν of the skeleton as:

    3 1 2 2 1 EE KG    (2.3)
Equation (2.1) can be interpreted as follows: on the left hand side is the total stress, while on the right hand side is the stress supported by the skeleton, i.e., the first two terms, and the stress supported by the fluid, which is the last term. Thus, the effective stress, which is the stress supported by the skeleton, is expressed as:

  "" ,0 0 2 23 3 ij ij kk ij ij ij K G G K T T                 (2.4)
Comparing equation (2.4) with equation (2.1), the effective stress can be expressed as function of the total stress and the pore pressure as:

" ij ij ij bp     (2.5)
Equation (2.5) is the expression for Biot's effective stress, which applies when the solid phase is either compressible or incompressible. For an incompressible solid phase (i.e., K S →∞), according to equation (2.2) the Biot's coefficient is b = 1, and thus equation (2.5) becomes the familiar expression for Terzaghi's effective stress:

" ij ij ij p     (2.6)
All the field quantities (stresses, strains, pressure, and temperature) in equation (2.1) are functions of coordinates 1 2 3 ( , , ) x x x x and time t. Under conditions of small perturbations, the strain ij  can be related to the displacement components i u as:

  ,, 1 2 ij i j j i uu   (2.7)
where the comma denotes the partial derivative with respect to the spatial variable.

The skeleton of the rock mass must always be in a state of equilibrium.

k ki i kk i i G K u Gu K T bp         
(2.9)

The unknown variables in equation (2.9) include the displacements i u , the pore pressure p and the temperature T.

Pore fluid flow

The continuity equation for the pore fluid under small perturbations conditions is given by, [START_REF] Coussy | Poromechanics[END_REF][START_REF] Selvadurai | Thermo-poroelasticity and Geomechanics[END_REF]:

  1 3 . 0 kk m pT bv t M t t             (2.10)
where M is the Biot modulus; α m is a coefficient; v  is the fluid velocity vector; and .()

 is the divergence operator.
The Biot modulus M is related to the bulk modulus of the fluid f K , the bulk modulus of the solid phase S K , the initial porosity 0  , and the Biot's coefficient b as:

00 1 fS b M K K    (2.11)
The coefficient α m is a function of the linear thermal expansion coefficients α s and α f of the fluid and the solid phases [START_REF] Selvadurai | Thermo-poroelasticity and Geomechanics[END_REF]:

  m f S b       
(2.12)

The fluid flow within the pores is assumed to follow Darcy's law. The velocity vector v  is expressed as:

k vp      (2.13)
where k is the hydraulic conductivity, in m 2 ; η is the dynamic viscosity of the fluid, in Pa.s Substituting equation (2.13) into equation (2.10) leads to the field equation that governs the fluid flow in the pores:

2 1 3 kk m p T k bp t M t t             (2.14)
It should be noted from equation (2.14) that the partial derivative with respect to time of the volumetric strain kk  shows the coupling between the mechanicals response of the skeleton and the fluid diffusion.

Fluid flow in hydraulic fractures and existing fractures

In permeable reservoirs, fluid loss from fractures is considered to be a pressure-independent process. The Carter's leakoff model gives an expression for the infiltration rate of fluid from the fracture as [START_REF] Carter | Optimum fluid characteristics for fracture extension[END_REF]:

0 2 ( , ) L C gt tt   x (2.15)
where ( , ) gt x is the infiltration rate from both faces of the hydraulic fracture;

L C is the Carter's leakoff coefficient; 0 t is the time at which the hydraulic fracture arrives at the point x.

The Carter's model is obtained from the pressure diffusion equation, assuming that the leakoff is one-dimensional. The model is used to describe the physical phenomenon that occurs as the fracturing fluid deposits a thin layer of material on the fracture faces. This phenomenon of material deposition is usually observed in hydraulic stimulation of oil and gas reservoirs where proppants and chemical additives (e.g., polymers) are used. In hydraulic stimulation of deep geothermal reservoirs, fresh water or brine is usually used. Thus, we assume that there is no material deposition on the surfaces of the hydraulic fracture. Therefore, the leakoff flow rate will depend only on the permeability of the rock mass and the difference in pressure between the inside of the fracture and inside of the reservoir, and can be expressed as:

0 () fr g C p p  (2.16)
where C 0 is the constant leakoff coefficient; p f and p r are the pressure inside the fracture and inside the reservoir (Figure 2.2), respectively.

With leakoff included, the continuity equation for the incompressible fluid flow in the fracture is given by: 0 wq g ts

     (2.17)
where w is the fracture aperture; q is the longitudinal fluid flow rate, as shown in Figure 2.2.

When two surfaces of the fracture are considered as parallel and smooth plates and the fluid is Newtonian and incompressible, the longitudinal flow rate q in steady state is related to the pressure gradient as (i.e., the Poiseuille equation):

3 12 f p w q s     (2.18)
where η is the dynamic viscosity of the fluid. From equation (2.18), the transmissivity of the fracture is defined as:

Heat transfer

As stated previously, hydraulic stimulation of EGS reservoirs involves temperature changes of significant magnitudes. Rock formations, e.g., crystalline rock masses, where EGS reservoirs are developed, usually have a very low permeability. In these formations, the predominant mode of heat transfer in the rock mass is conduction. The convective heat transfer associated with the slow diffusion fluid flows is negligible [START_REF] Selvadurai | Thermo-poroelasticity and Geomechanics[END_REF]. Therefore, Fourier's law of heat conduction can be used to describe the heat transfer in the rock matrix as:

qT     (2.21)
where q is the heat flux vector; λ is the thermal conductivity of the rock mass.

The thermal energy conservation equation for the rock mass is given as:

.0 At the interface between a flowing fluid and a solid, such as fluidrock interface, heat transfer takes place by convection. The basic relationship for heat transfer by convection is:

P T Cq t      (2.
() a q n h T T    (2.24)
where q is the heat flux vector, n is the unit normal vector, T a is the temperature of the fluid, T is the temperature on the surface of the solid, and h is the convective heat transfer coefficient

FRACTURE PROPAGATION

For brittle rocks such as granite under in-situ conditions of stress, fracture is the main mechanism of failure. Other processes that may lead to failure include corrosion and wear.

However, these processes happen slowly over long time periods and are less of a concern, especially for rock masses where EGS rock reservoirs are to be created. In general, the fracture process is divided into three distinct phases [START_REF] Broberg | Cracks and Fracture[END_REF]: (1) loading without crack growth, (2) stable crack growth, and (3) unstable crack growth. Stable crack growth can be controlled by controlling the applied load and a prescribed slow crack growth may be obtained. This is not possible for unstable crack growth, which occurs spontaneously and is usually accompanied induced dynamic effects such as elastic waves.

In the vicinity of the tip of a crack, two distinctly different zones exist, as shown in Figure 2.3.

These include a small zone immediately next to the crack tip known as the process zone and a plastic zone outside the process zone. The fracture process always starts from the process zone, where stress and strain concentrations occur. When subjected to high stress, the nucleation and growth of micro-separations takes place in the process zone, i.e., material separation, such as the breakage of atomic bonds. Micro-separations are nucleated at inhomogeneities in the material, e.g., flaws or boundaries between grains. As these microseparations extend and propagate, voids are created, which in turn propagate and coalesce with other voids and with the main crack, allowing the main crack to grow. The nucleation and growth of microseparations take place in the process zone, leading to the creation of voids. The propagation and coalescence of voids allow the crack to grow. After [START_REF] Broberg | Cracks and Fracture[END_REF] For polycrystalline materials such as rocks, the micro-separations are in the form of microcracks. The nucleation and growth of these micro-cracks in the process zone is the kernel for the propagation of macro-crack in these materials. [START_REF] Atkinson | Fracture Mechanics of Rock[END_REF] used the loading of a blunt notch, as shown in

Fracture mechanics

The theory of fracture mechanics revolves around two major approaches, namely the Griffith energy balance approach, which was developed by [START_REF] Griffith | The phenomenon of rupture and flow in solids[END_REF], and the Irwin stress analysis approach, which was introduced by [START_REF] Irwin | Analysis of stresses and strains near the end of a crack traversing a plate[END_REF]. These approaches and related issues are presented in detail in the following paragraphs.

Griffith energy balance approach

Griffith used thermodynamic arguments to find the criterion for a crack to grow due to an external load. By using the results of stress analyses of a plate containing a crack-like elliptical hole [START_REF] Inglis | Stresses in a plate due to the presence of cracks and sharp corners[END_REF], Griffith showed a quantitative connection between the fracture stressthe applied stress at which the crack starts to propagate, and the crack size [START_REF] Griffith | The phenomenon of rupture and flow in solids[END_REF].

According to the first law of thermodynamics, the total energy of a system remains constant as the system moves from one equilibrium state to another. Griffith argued that the quasi-static propagation of a crack can be seen as occurring under equilibrium conditions. Thus, during the propagation of the crack the decrease in strain energy must be sufficient to overcome the energy necessary to create new crack surfaces.

Consider a Mode I fracture as shown in Figure 2.6, the total energy of the system is composed of the potential energy  and the total surface energy s W . Thus:

s EW    (2.25)
The potential energy  , which is supplied by the internal strain energy and the work of the external loads, is given by:

1 2 ij ij ij j S d n dS          (2.26)
The energy release rate G I , which is a measure of the energy available for an increment of crack extension in Mode I, is defined as:

I d G dA   (2.27)
where A is the crack area

The Mode I fracture energy G IC , which is the energy necessary to create new crack surfaces for the same increment of crack extension, is defined as:

s IC dW G dA 
(2.28)

The Griffith energy balance for an incremental increase in the crack area dA can be written as:

Irwin stress analysis approach

This stress analysis-based approach for studying fracture was introduced by [START_REF] Irwin | Analysis of stresses and strains near the end of a crack traversing a plate[END_REF], who recognized that for certain cracked configurations subjected to external loads, it is possible to obtain closed-form expressions for the stresses in the region near the crack tip, assuming that the material is isotropic and linear elastic. The stresses at a point M near the crack tip and located by the polar coordinates (r, θ), as shown Figure 2.7, are given by:

2a σ σ () 2 ij ij K f r    (2.36)
where K is the stress intensity factor that depends on the loading mode; f ij (θ) is a well defined dimensionless function of θ and depends on the loading mode. Functions f ij (θ) for Modes I and II are given in Table 2.1

f ij (θ) Mode I Mode II f xx 3 cos 1 sin sin 2 2 2                          3 sin 2 cos cos 2 2 2                          f yy 3 cos 1 sin sin 2 2 2                          3 sin cos cos 2 2 2                      f xy 3 cos sin cos 2 2 2                      3 cos 1 sin sin 2 2 2                         
There are three stress intensity factors K I , K II , and K III that correspond to the three loading modes I, II, and III. In a mixed-mode problem (i.e., when more than one loading mode is present), the principle of linear superposition applies and the stresses are given as:

( ) ( ) ( ) I II III ij ij ij ij       
(2.37) For a given loading condition, the stress intensity factor is the only parameter that characterizes the stress field, and thus the strain field, of the region near a crack tip. The propagation of the crack depends entirely on the stress intensity factor. For instance, under pure Mode I loading, the crack propagation criterion is expressed as:

I IC KK  (2.38)
where K IC is the critical stress intensity factor, which is unique to each material and can be obtained through experiments. Expressions similar to equation (2.38) also apply to pure fracture Modes II and III.

Relationship between the energy release rate and stress intensity factor

The fracture criterion can be expressed as the energy release rate as in equation ( 2 the more convenient way is to use the fracture criterion defined in terms of the energy release rate in a generalized form because there exists a relationship between the stress intensity factor and the energy release rate for mixed-mode loading conditions, which is given by, [START_REF] Anderson | Fracture Mechanics: Fundamentals and Applications[END_REF]:

2 2 2
''

I II III K K K G EE     (2.41)
where μ is the shear modulus. Equation (2.41) can be interpreted as: under mixed-mode conditions, the total energy release rate G is equal to the sum of work G I , G II , and G III necessary to close the crack along the normal, in-plane, and out-of-plane directions, respectively, i.e.,:

I II III G G G G    (2.42)
where 2 2 2 ; ; ''

I II III I II III K K K G G G EE     (2.43)
The fracture criterion under mixed-mode loading conditions now becomes:

C GG  (2.44)
where G C is the equivalent fracture energy.

Equivalent fracture energy

Equation (2.44) is similar to equation (2.40), except that the equivalent fracture energy G C is no longer a constant, depending not only on the material as is the case of G IC (or G IIC , G IIIC ), but also on the proportions between G I , G II , and G III . Several formulae have been proposed to compute G C , including the following:

 The BK law [START_REF] Benzeggagh | Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus[END_REF])

  m II III C IC IIC IC I II III GG G G G G G G G         (2.45)
where m is a material parameter  The Power law [START_REF] Wu | Crack extension in fiberglass reinforced plastics[END_REF])

a b c I II III C IC IIC IIIC G G G G G G G G                      (2.46)
where a, b, c are material parameters  The Reeder law [START_REF] Reeder | Postbuckling and growth of delaminations in composite plates subjected to axial compression[END_REF])

    m III II III C IC IIC IC IIIC IIC II III G G G G G G G G G G G G                 (2.47)
where m is a material parameter

Crack orientation

Under mixed-mode loading conditions, the direction of crack propagation also needs to be determined; several proposals exist in the literature [START_REF] Erdogan | On the crack extension in plates under plane loading and transverse shear[END_REF][START_REF] Qian | Mixed mode fatigue crack growth: a literature survey[END_REF][START_REF] Sih | Mixed mode fatigue crack growth predictions[END_REF]. The propagation direction determined based on the maximum tangential stress or minimum strain energy is presented below.

 Maximum tangential stress criterion:

This criterion states that the crack will propagate along the radial direction determined by c   on which the tangential stress σ θ becomes maximum. Mathematically, θ c is determined from: (2.49)

Comments on the Griffith energy and Irwin stress approach

Both the Griffith energy balance and the Irwin stress analysis approaches form the foundation of linear elastic fracture mechanics. Both are based on the assumption that the material is linear elastic and thus singular stresses are allowed to exist at the crack tip. This assumption cannot be true for real materials, which all have a certain stress limit. Therefore, for real materials plastic deformation will always occur in the region around the crack tip where there are high stress concentrations. However, as long as the plastic zone is small in relation to the crack size, its effect is negligible and the above approaches are considered acceptable. This is true for rocks, especially brittle rocks that show very little plastic deformation before fracture occurs, as demonstrated in many experiments [START_REF] Keshavarz | Contribution to experimental study of mechanical and thermal damage in crystalline hard rocks[END_REF][START_REF] Keshavarz | Damage and changes in mechanical properties of a Gabbro thermally loaded up to 1,000°C[END_REF].

When plastic deformation is no longer negligible, the fracture criteria proposed by Griffith and

Irwin are no longer valid. The energy consumed by plastic deformation has to be taken into consideration. In these cases, the fracture criterion can be formulated based on the concept of the J-integral, a more general concept introduced by Rice (1968) that can be applied to both linear and non-linear materials in monotonic loading conditions. With a crack shown in Figure 2.9, the J-integral is defined as [START_REF] Rice | A path independent integral and the approximate analysis of strain concentration by notches and cracks[END_REF]:

i i u J wdy T ds x         (2.50)
where w is the strain energy density; T i are components of the traction vector; u i are components of the displacement vector; ds is the length increment along the contour Γ. [START_REF] Rice | A path independent integral and the approximate analysis of strain concentration by notches and cracks[END_REF] showed that the value of the J-integral is independent of the path of integration around the crack and is equal to the energy release rate G defined in the previous section. In fracture problems where poroelastic effects are involved (e.g., fracture propagation due to hydraulic stimulation), the poroelasticity needs to be considered. It was shown by [START_REF] Atkinson | Plane strain fracture in poroelastic media[END_REF] that when taken into account the poroelastic effects retard the propagation of a fracture as pore pressures reduce the stress intensity factor. Using computational modeling, [START_REF] Selvadurai | Computational modeling of the indentation of a cracked poroelastic half-space[END_REF] and [START_REF] Mahyari | Enhanced consolidation in brittle geomaterials susceptible to damage[END_REF] also showed that the time-dependent variation of stress intensity factors at a propagating fracture can be computed using the transient analysis of the associated hydromechanical problem. For instance, for the problem of fracture propagation due to hydraulic stimulation the transient coupled stress-diffusion analysis is performed first and the effective stresses are then used to compute the stress intensity factors.

Cohesive zone model

An alternative for modeling the fracture is the cohesive zone model (CZM). This model was introduced by [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF], who tried to find a way to describe the damage zone at the crack tip. The model is based on the observation that there is always a partially damaged zone at the tip of a crack (i.e., the cohesive zone) through which stress can continue to be transferred until the separation between the two faces of this cohesive zone reaches a certain value. The cohesive zone (Figure 2.9) at the crack tip is the transition zone between the traction-free crack and the intact material and can be seen as the process zone. For the crack presented in Figure 2.9, the separation reaches a critical displacement δ f at the tip of the traction-free crack. The tractions are zero at this point, but are equal to the tensile strength T R at the tip of the cohesive zone. The cohesive zone model is defined by a tractionseparation law. The damage initiation occurs when the traction reaches the tensile strength of the material. As the external load is increased the separation increases while the traction across this cohesive surface decreases and eventually vanishes. The traction-free crack is thus created and allows the crack to grow. A damage variable is used to track the damage of the cohesive zone during the separation process. The damage variable is zero until damage initiation, monotonically increases during the damage evolution, and reaches one as separation reaches the critical separation.

An example of the cohesive model for Mode I fracture is shown in Figure 2.10. This linear cohesive model can be described by three parameters: the tensile strength R T , the critical separation δ f and the initial stiffness K. The energy dissipated during the separation process for a unit increment of crack extension is by definition the fracture energy. It is obvious from Figure 2.10 that:

1 2

IC T f GR   (2.51)
The normal traction that can still be transferred across the cohesive zone during the separation process is calculated as:

(1 ) p D  

(2.52)

where D is the damage variable and σ p is a traction predictor. These two parameters are defined as (with 0 /

T RK   ) 0 0 () () f p f DK           (2.

DYNAMIC FRACTURE PROPAGATION AND WAVE PROPAGATION

The propagation of brittle fractures is inherently a dynamic process that involves the release of strain energy in the form of elastic waves [START_REF] Hardy | Acoustic Emission/Microseismic Activity: volume 1: Principles, Techniques and Geotechnical Applications[END_REF][START_REF] Pollock | Acoustic emission -2: Acoustic emission amplitudes[END_REF][START_REF] Scholz | Experimental study of the fracturing process in brittle rock[END_REF]. The criterion for Mode I fracture propagation under dynamic conditions is written as follows [START_REF] Anderson | Fracture Mechanics: Fundamentals and Applications[END_REF]:

() I ID K t K  (2.54)
where K I (t) is the instantaneous stress intensity factor, K ID is the dynamic fracture toughness of the material. In the dynamic context, K I (t) depends not only on the current external loads but also the reflecting stress waves that interact with the local stress field at the crack tip while K ID is a function of the fracture propagation speed V F . The dependence of K ID on the fracture propagation speed can be expressed by an empirical formula as follows [START_REF] Anderson | Fracture Mechanics: Fundamentals and Applications[END_REF]:

1 1 m F ID IA L V KK V          (2.55) Separation Traction R T δ f D = 0 D = 1 K δ 0 1 2 IC T f GR   δ σ 54
where V F is the fracture propagation speed, V L is the limiting fracture propagation speed in the material which is usually taken to be equal to the speed of Rayleigh waves [START_REF] Freund | Dynamic Facture Mechanics[END_REF], m is a constant, and K IA is the arrest toughness of the material corresponding to V F = 0.

The equation of wave propagation is the equilibrium equation taking into account the inertia forces and is written as:

2 2 . t       u F (2.56)
where u is the displacement vector; F is the body forces; and  is the Cauchy stress tensor.

Under the condition of small perturbations and a linear elastic material, one has the following relations:

 Kinematic relation between strains and displacements:

  1 2 T      uu (2.57)
 Linear relation between stresses and strains:

( ) 2 tr I      (2.58)
where ()  is the gradient operator; λ and μ are the Lamé coefficients.

Combining equations (2.56), (2.57), and (2.58) leads to (without body force, F = 0):

2 2 2 ( ) ( . ) t             u uu (2.59)
Using the Helmholtz decomposition the displacement vector can be decomposed as follows:

with

P PS S        u u u u u rot ψ (2.60)
where φ is a scalar potential; and ψ is vector potential; rot () is the rotation operator. It's worth noting that the displacement vector u P is irrotational (i.e., rot u P = 0) and the divergence of u S is zero. Using these properties of u P and u S and equation (2.59), the uncoupled wave equations can be written in terms of u P and u S as follows:

2 22 2 2 22 2 0 0 P PP S SS C t C t                 u u u u (2.61)
where C P and C S are the speed of dilatational waves (or longitudinal waves) and shear waves (transverse waves), respectively, and are determined by: 2 and

PS CC      
(2.62)

The wave equations can also be written in terms of the two potentials φ and ψ as follows:

2 22 2 2 22 2 0 0 P S C t C t                   ψ ψ
(2.63)

FAULT REACTIVATION

Faults are shear fractures of rock at geological scales. Depending on the role of the vertical stress relative to the fault shear displacement, faults are classified into three major categories and are illustrated in Figure 2.11. These include [START_REF] Cornet | Elements of Crustal Geomechanics[END_REF][START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF][START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF]: The block is subjected to two external forces: a normal force N normal to the surface and a shear force T parallel to the surface. Two reaction forces (Figure 2.12) appear to counterbalance the applied forces: the normal reaction force Q and shear resistance R. As long as the shear force T does not exceed the shear resistance R the solid stay still, or in other word no relative slip between the solid and the surface occurs. Relative slip is initiated when the balance between the shear force T and the shear resistance R is broken. Experiments showed that the shear resistance R is proportional to the normal force N via a coefficient called friction coefficient μ [START_REF] Engelder | The role of asperity indentation and ploughing in rock friction-II: Influence of relative hardness and normal load[END_REF][START_REF] Scholz | Experimental study of the fracturing process in brittle rock[END_REF]:


RN   (2.64)
The slip criterion is expressed as:

TN   (2.65)
By dividing both sides of equation (2.65) by the effective area of contact, the slip criterion expressed in shear and normal stress τ and σ is obtained as follows:

   (2.66) σ 1 = σ v σ 1 σ 3 σ 3 σ 3 = σ v σ 3 σ 1 σ 1 σ 3 σ 3 σ 1 σ 1 57
The right hand side of equation (2.66) can be regarded as the local shear strength. The shear strength determined by equation (2.66) is usually referred to as the Coulomb friction criterion.

It's worth noting that when pore pressure is present, the normal stress σ is taken as the effective stress, which wasdefined in equations (2.5) and (2.6), and is repeated here for an incompressible solid phase as:

' p   (2.67)
Thus the Coulomb friction criterion becomes:

() p     (2.68) Figure 2.
12 Solid block on a horizontal surface, subjected to a force N normal to the surface and a force T parallel to the surface. R is the shear resistance Equation (2.66) is usually used as a standard criterion for assessing fault reactivation.

Equation (2.66) is obtained by considering the global equilibrium of the whole block.

However, equation (2.66) itself is expressed in terms of quantities that are defined locallythe shear stress τ and the normal stress σ n . It is obvious that even though equation (2.66) is satisfied at several locations of the interface, it does not mean that the whole interface will slip. All frictional interfaces, including faults, establish contact through asperities e.g., on several asperities of the interface. Slow slip first initiates locally where resistance to rupture (i.e., shear strength) is the smallest along the interface or where shear stress is locally high.

The extension of the local slow slip to a certain state will ignite the global and potentially dynamic slip of the interface [START_REF] Selvadurai | Laboratory-developed contact models controlling instability on frictional faults[END_REF]. This phenomenon of slow slip occurring prior to dynamic slip was experimentally observed and investigated in a series of meticulously programmed experiments by [START_REF] Glaser | Asperity generation and its relationship to seismicity on a planar fault: a laboratory simulation[END_REF]. It should be noted

T N R Q
that the influence of scale effects on the shear strength of frictional interfaces also needs to be considered when assessing the stability of the interface [START_REF] Vallier | A shear model accounting scale effect in rock joints behavior[END_REF].

In geothermal reservoirs, the slip of pre-existing faults may be induced by the injection of fluid. The fault slip can be seismic, i.e., generating seismic waves or aseismic, i.e., slip without generating seismic waves. Thermo-poro-mechanical couplings play an important role in the stability of faults and the seismic slip [START_REF] Rattez | Effect of Thermo-Poro-Mechanical Properties on the Stability of Faults[END_REF][START_REF] Sulem | Multiphysics Couplings and Stability of Fault Zones[END_REF]. Many experiments on rock friction [START_REF] Byerlee | Friction of rocks[END_REF][START_REF] Byerlee | Stick slip, stable sliding, and earthquakes-effect of rock type, pressure, strain rate, and stiffness[END_REF][START_REF] Hoskins | A medium-scale direct friction experiment[END_REF] have shown that for certain types of rock and surface roughness, the relative displacement of a shear fracture can occur in a stick-slip manner, as illustrated in Figure 2.13. Under a constantheld normal load N, as the shear force T applied along a shear fracture increases the shear elastic relative displacement increases, until, at some point, unstable displacement occurs abruptly with a concurrent drop in the shear force. The fracture surfaces subsequently lock together and the shear force can increase and the stick-slip cycle occurs again. The friction of rock is therefore characterized by two different friction coefficients: a coefficient for the stick phase, called coefficient of static friction μ, and a coefficient for the slip phase, known as the coefficient of dynamic friction μ d . [START_REF] Brace | Stick-slip as a mechanism for earthquakes[END_REF] suggested that the stick-slip phenomenon may provide a mechanism for generation of seismic waves from fault slip. [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF] Rock friction is also dependent on time and slip rate, which was experimentally demonstrated by [START_REF] Dieterich | Time-dependent friction in rocks[END_REF] and [START_REF] Jafari | Evaluation of shear strength of rock joints subjected to cyclic loading[END_REF]. The time-and rate-dependency of the rock friction is expressed through the rate-and-state model as follows:

d 59 0 0 0 ln CC V ab V V dV e dt D D            (2.69)
where V is the slip rate; μ 0 is the steady state friction coefficient at a slip rate V = V 0 ; a and b are friction parameters; D C is the characteristic slip distance; and Θ is a state variable. The relation between a and b dictates whether the fault is velocity-weakening or velocity strengthening and thus whether the slip is stable or unstable.

The time-dependency of the friction coefficient in the rate-and-state model expressed through the state variable is attributed to the creep effect that causes an increase in the actual area of asperity contact between the two faces of the fault [START_REF] Engelder | The role of asperity indentation and ploughing in rock friction-II: Influence of relative hardness and normal load[END_REF]. The slip ratedependency is believed to be caused by the loss and recovery of contact between asperities of the fault as its two faces slip relative to each other.

The rate-and-state friction model has been used to model the slip behavior of a single fault due to direct fluid injection by Cueto-Felgueroso et al. ( 2017), [START_REF] Dublanchet | Interactions and triggering in a 3-D rate-and-state asperity model[END_REF] and [START_REF] Mcclure | Investigation of injection-induced seismicity using a coupled fluid flow and rate/state friction model[END_REF]. The model is able to explain various observational features of earthquake phenomena [START_REF] Segall | Earthquake and Volcano Deformation[END_REF]. However, recent laboratory experiments on fault slip due to fluid injection [START_REF] Scuderi | Frictional stability and earthquake triggering during fluid pressure stimulation of an experimental fault[END_REF] concluded that slip instabilities can be induced whenever the critical stress state for reactivation is met by an increase in fluid pressure, even if the fault is characterized by velocity strengthening frictional behavior. Moreover, these authors stated that the decrease of shear strength due to the reduction in the effective normal stress largely prevails over the impact of the rate-and-state effects on fault frictional strength.

In the current research, the local shear strength of natural fractures or faults is modeled using the Coulomb friction law. This law is implemented in the finite element code ABAQUS using the subroutine USDFLD and the field variable (ABAQUS, 2016).

HYDRAULIC FRACTURE MODELS AND INTERACTION OF FRACTURES

Hydraulic fracture models

As stated in the previous chapter, the development of a rock reservoir is crucial to the success of a EGS project. Excessive permeability will lead to an unacceptable loss of water from the system, while high impedance (i.e., resistance to flow) will require a large amount of energy to pump water through the reservoir. Thus, in order to estimate the permeability of the reservoir it is necessary to predict the development of fracture geometry based on the stimulation conditions. Three single hydraulic fracture models have been derived to approximately define the development of the fracture geometry. These include two plane strain modelsthe Khristianovic-Geertsma-de Klerk (KGD) model [START_REF] Geertsma | A rapid method of predicting width and extent of hydraulically induced fractures[END_REF] and the Perkins-Kern-Nordgren (PKN) model [START_REF] Nordgren | Propagation of a vertical hydraulic fracture[END_REF], and a three-dimensional axisymmetric modelthe radial penny-shaped model [START_REF] Sneddon | The distribution of stress in the neighbourhood of a crack in an elastic solid[END_REF]. Figure 2.14 illustrates the geometry of these three fracture models. Early solutions for the fracture geometry (e.g, fracture length, aperture) and wellbore pressure have been derived for those models by assuming that (i) the material is homogeneous, isotropic, linear elastic, and fully impermeable, (ii) the fluid is Newtonian and its flow within the fracture is laminar everywhere, (iii) the fluid injection is considered as a line source for the KGD and KPN models and as a point source for the penny-shaped model, and (iv) linear elastic fracture mechanics applies [START_REF] Geertsma | A rapid method of predicting width and extent of hydraulically induced fractures[END_REF][START_REF] Geertsma | A Comparison of the Theories for Predicting Width and Extent of Vertical Hydraulically Induced Fractures[END_REF]. These early solutions have been used by many researchers as a benchmark for comparison with their numerical models. However, it should be noted that several inconsistencies were ignored in order to obtain those solutions. For instance, to obtain solutions of fracture geometry for the KGD and penny-shaped models, a smooth and finitecurvature fracture tip is assumed. This is obviously inconsistent with the linear elastic fracture mechanics, which considers the fracture tip as sharp and having infinite curvature.

More recently, these fracture models (KGD and penny-shaped models) were studied with the incorporation of more consistent considerations for the fracture tip and boundary conditions incorporated [START_REF] Bunger | Toughness-dominated hydraulic fracture with leak-off[END_REF][START_REF] Detournay | Propagation regimes of fluid-driven fractures in impermeable rocks[END_REF][START_REF] Savitski | Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions[END_REF]. For instance, the fracture tip is treated as having a zero-aperture with no fluid flow. In these studies, the following assumptions were made: (i) the medium is infinite, (ii) the material is fully impermeable and linear elastic, (iii) linear elastic fracture mechanics applies, and (iv)

Carter's leakoff model is used to describe the fluid loss across the walls of the hydraulic fracture. With these assumptions, the governing equations of the problem then reduce to (i) the equilibrium equation for the linear elastic material, which can be given by an integral that relates the hydraulic aperture and the net fluid pressure [START_REF] Spence | Self-similar solutions for elastohydrodynamic cavity flow[END_REF], (ii) the local and global mass balance equations for the injected fluid, and (iii) the fracture propagation criterion, which is based on linear elastic fracture mechanics that relates the fluid pressure and the fracture energy. A scaling technique was then used to obtain approximate analytical solutions in the form of regular asymptotic expansions [START_REF] Bunger | Toughness-dominated hydraulic fracture with leak-off[END_REF][START_REF] Detournay | Propagation regimes of fluid-driven fractures in impermeable rocks[END_REF][START_REF] Garagash | Plane-strain propagation of a fluid-driven fracture during injection and shut-in: Asymptotics of large toughness[END_REF][START_REF] Savitski | Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions[END_REF]. It was shown that different regimes of fracture propagation exist, depending on the energy dissipation mechanisms and the fluid balances. Energy can be dissipated either due to the flow of the viscous fluid (viscositydominated regime) or the propagation of the fracture (toughness-dominated regime).

Meanwhile, the fluid can either be stored in the fracture (storage-dominated regime) or leak into the surrounding medium (leakoff-dominated regime). These four regimes of fracture propagation can be represented in a two dimensional space where each limiting regime has In the following paragraphs, the main results for fracture geometry and injection pressure for the KGD and penny-shaped fracture models near the M and K vertices (viscosity/storagedominated regime and toughness/storage-dominated regime) are briefly summarized. These regimes were chosen because they correspond to real conditions that could be encountered in deep geothermal reservoirs where the granitic basement has a very low permeability, which corresponds to the storage-dominated regime. Under these conditions, the approximate solutions for the fracture geometry and wellbore pressure for both the KGD and the pennyshaped fracture models are given in a generic form as [START_REF] Detournay | Propagation regimes of fluid-driven fractures in impermeable rocks[END_REF][START_REF] Savitski | Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions[END_REF]:

 Fracture length: ( ) ( ) L t L t   (2.70)  Fracture aperture: ( , ) ( ) ( ) w x t t L t   (2.71)  Wellbore pressure: M K M K Storage-dominated Leakoff-dominated Viscosity dominated toughness dominated 0   Energy Mass ( ) ' p t E   (2.72)
In the above equations, L(t) is the length scale of the same order as the fracture length; expressions for () t  , L(t), and the dimensionless numbers γ, Ω, and  are defined below, dependent on the fracture model and the propagation regime; E' is defined as:

2 ' 1 E E    (2.73)
where E and ν are the Young's modulus and Poisson's ratio of the medium, respectively.

The following definitions are also introduced for convenience:

1/2 2 ' 12 , ' 4
IC KK       (2.74)
where μ is the viscosity of the fluid and K IC is the critical Mode I stress intensity factor

KGD model

Viscosity-dominated regime (M vertex):

Expressions for () t  and L(t) in the viscosity-dominated regime are given by Detournay ( 2004):

1/6 1/3 34 0 ' ' , '' mm E Q t L Et           (2.75)
where Q 0 is the injection rate per unit height of the fracture (m 2 /s). By combining equation (2.75) with equations (2.70), (2.71), and (2.72), the expression for the fracture geometry and the wellbore pressure can be obtained as:

  1/6 23 1/3 0 0 1/6 3 2/3 0 0 1/3 2 1/3 0 ' ' ' ' '' m m m Q wt E EQ Lt p E t               (2.76)
where the dimensionless numbers γ m0 , Ω m0 , and 0 m  can be expressed in the form of series expansion and can be found in [START_REF] Adachi | Self-similar solution of a plane-strain fracture driven by a power-law fluid[END_REF].

Toughness-dominated regime (K vertex):

Expressions for () t  and L(t) in the toughness-dominated regime are given by Detournay ( 2004):

1/3 2/3 4 0 4 0 ' ' , '' kk E Q t K L E Q t K        
(2.77) By combining equation (2.77) with equations (2.70), (2.71), and (2.72), the expressions for the fracture geometry (aperture and length) and the wellbore pressure are:

1/3 2 1/3 0 0 2 2/3 2/3 0 0 1/3 4 1/3 0 0 ' ' ' ' ' ' k k k KQ wt E EQ Lt K K pt EQ               (2.78)
where the dimensionless numbers γ, Ω, and  are given by:

  1/3 1/2 2 0 0 0 2/3 1/3 21 , , 1 8 k k k          (2.79)
with ρ = x/L, x and L shown in Figure 2.14.

Transition between propagation regimes:

A dimensionless number to distinguish the propagation regimes is defined as [START_REF] Adachi | Self-similar solution of a plane-strain fracture driven by a power-law fluid[END_REF][START_REF] Detournay | Propagation regimes of fluid-driven fractures in impermeable rocks[END_REF]:

  1/4 3 0 ' '' K K EQ   (2.80)
Strictly speaking, the viscosity-and toughness-dominated regimes of propagation correspond to K ≪ 1 and K ≫ 1, respectively. However, according to [START_REF] Detournay | Propagation regimes of fluid-driven fractures in impermeable rocks[END_REF], the fracture can be considered as propagating in the viscosity-dominated regime when K < 1 and in the toughness-dominated regime when K > 4.

Penny-shaped fracture model

Two regimes of propagation are also considered for the penny-shaped fracture model: the viscosity-and toughness-dominated regimes.

Viscosity-dominated regime (M vertex):

Expressions for () t  are L(t) in the viscosity-dominated regime are given by [START_REF] Savitski | Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions[END_REF]:

1/9 1/3 34 0 ' ' , '' mk E Q t L Et           (2.81)
where Q 0 is the injection rate in m 3 /s By combining equation (2.81) with equations (2.70), (2.71), and (2.72), the expressions for the fracture geometry and the wellbore pressure can be obtained as:

  1/9 23 1/9 0 0 2 1/9 3 4/9 0 0 1/3 2 1/3 0 ' ' ' ' '' m m m Q wt E EQ Lt p E t               (2.82)
where the dimensionless numbers γ m0 , Ω m0 , and 0 m  can be expressed in the form of a series expansion and can be found in O' Keeffe et al. (2018) and [START_REF] Savitski | Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions[END_REF].

Toughness-dominated regime (K vertex):

Expressions for () t  and L(t) in the toughness-dominated regime are given by [START_REF] Savitski | Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions[END_REF]:

1/5 2/5 6 0 6 0 ' ' , '' kk E Q t K L E Q t K         (2.83)
By combining equation (2.83) with equations (2.70), (2.71), and (2.72), the expression for the fracture geometry (aperture and length) and the wellbore pressure are:

1/5 4 1/5 0 0 4 2/5 2/5 0 0 1/5 6 1/5 0 0 ' ' ' ' ' ' k k k KQ wt E EQ Lt K K pt EQ               (2.84)
where the dimensionless numbers γ k0 , Ω k0 , and 0 k  are given by:

  2/5 1/5 1/5 1/2 2 0 0 0 33 , , 1 8 12 8 2 k k k                          (2.85) with ρ = r/R, r and R shown in Figure 2.14

Transition between propagation regimes:

A dimensionless number K used to distinguish different propagation regimes is also defined as [START_REF] Savitski | Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions[END_REF]:

1/18 2 5 3 13 0 ' '' t KK QE      (2.86)
The viscosity-and toughness-dominated regimes correspond to K < 1 and K > 3.5, respectively. Unlike the KGD model, the dimensionless number K in this case depends not only on the injection rate and material properties but also on the injection time. From equation (2.86), it is obvious that the fracture propagates in the viscosity-dominated regime at an earlier time (when t is small) and gradually transitions to the toughness-dominated regime at an later time.

Interaction of hydraulic fractures and natural fractures

It has been recognized that natural fractures may have an important impact on the development of EGS reservoirs, especially after the development of the Camborne School of Mines project, where instead of generating a vertical fracture propagating upwards as predicted the reservoir migrated downwards into the region where there were multiple natural fractures. The interaction between natural fractures and the advancing hydraulic fractures is a key factor that produces complex fracture patterns. Even in rock formations where the natural fractures are sealed or cemented, the natural fractures may still serve as weak pathways to potentially divert advancing hydraulic fractures (Dahi [START_REF] Dahi Taleghani | How natural fractures could affect hydraulic-fracture geometry[END_REF].

When a hydraulic fracture approaches a natural fracture, two potential scenarios can occur, as Figure 2.16 illustrates: (i) the hydraulic fracture is diverted into the natural fracture, or (ii) the hydraulic fracture crosses the natural fracture. The hydraulic fracture may also be "arrested" by the natural fracture under specific conditions (e.g., frictional strength of the natural fracture). However, as long as the fluid continues to be injected, the hydraulic fracture will eventually continue to propagate, either crossing the natural fracture or being diverted into it.

The interaction between hydraulic fracture and natural fracture has been studied using both experimental, theoretical, and numerical simulation approaches. The main objective is to predict under which conditions the hydraulic fracture will cross or divert into the natural fracture. Theoretical studies, primarily based on the theory of linear elastic fracture mechanics, have focused on quantifying the alteration of stress intensity factor at the tip of the hydraulic fracture caused by the stress perturbations due to the presence of the natural fracture in the medium [START_REF] Gu | Criterion For Fractures Crossing Frictional Interfaces At Nonorthogonal Angles[END_REF][START_REF] Taleghani | Overview of numerical models for interactions between hydraulic fractures and natural fractures: Challenges and limitations[END_REF][START_REF] Thomas | Quantification of fracture interaction using stress intensity factor variation maps[END_REF][START_REF] Weng | Impact of Preexisting Natural Fractures on Hydraulic Fracture Simulation[END_REF].

One such study by [START_REF] Renshaw | An experimentally verified criterion for propagation across unbounded frictional interfaces in brittle, linear elastic materials[END_REF] (2.87) Equation (2.87) suggests that the crossing would always occur when the friction coefficient μ is larger than 0.5, which is a typical value for friction coefficient of rocks [START_REF] Cornet | Elements of Crustal Geomechanics[END_REF].

However, the interaction of fractures is known to be very sensitive to the intersection angle.

(Gu and Weng, 2010) used a similar approach to obtain a crossing criterion for non-orthogonal intersection conditions. They found that as the intersection interaction decreases from 90° (the orthogonal intersection case) the hydraulic fracture is more likely to be diverted into the natural fracture. A plot of their crossing criterion is shown in Figure 2.17. It is obvious from

Figure 2.17 that when the fracture intersection angle is smaller than 30°, the diversion of the hydraulic fracture into the natural fracture always occurs in most practical situations.

Experiments on the interaction between hydraulic fractures and natural fractures under different conditions of stress and friction were carried out by [START_REF] Blanton | An Experimental Study of Interaction Between Hydraulically Induced and Pre-Existing Fractures[END_REF], [START_REF] Renshaw | An experimentally verified criterion for propagation across unbounded frictional interfaces in brittle, linear elastic materials[END_REF], [START_REF] Warpinski | Influence of Geologic Discontinuities on Hydraulic Fracture Propagation (includes associated papers 17011 and 17074 )[END_REF], and [START_REF] Zhou | Experimental Investigation of Fracture Interaction between Natural Fractures and Hydraulic Fracture in Naturally Fractured Reservoirs[END_REF]. They found that the contrast in in-situ stresses and the angle of intersection are the most important factors that control the intersection behavior. Under normal conditions of stresses found in-situ, the hydraulic fracture tends to divert into the natural fracture; the crossing only occurs in cases of high differential stresses and high intersection angles.

NUMERICAL METHODS FOR MODELING FRACTURE PROPAGATION

Extended finite element method

The extended finite element method (XFEM) was introduced by [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF] and [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. The method allows discontinuities such as cracks to be represented independently from the finite element (FE) mesh. Thus, the crack growth is also independent of the FE mesh and remeshing is no longer needed. The central idea in the formulation of this method is that discontinuities are incorporated in a finite element by enriching the traditional finite element approximation with special enriched functions in conjunction with additional degrees of freedom. The enrichment functions typically consist of the near-tip asymptotic functions that capture the singularities around the crack tip and a discontinuous function that represents the jump in displacement across the crack surfaces. The approximation for a displacement vector at any location x is given by: 4 1 ( ) ( ) ( ) ( )

I I I I IN N H F            u x x u x a x b (2.
F r r r r r             (2.89)
where (r, θ) is a polar coordinate system with its origin at the crack tip and θ = 0 is tangent to the crack at the tip, as Figure 2.18 illustrates.

The Heaviside function H(x) defined with reference to Figure 2.18 is:

1 if ( *). () 1 otherwise H       x x n x (2.90)
where x is the location where displacement is currently approximated, and x* is the closest projection of x on the crack surface, n is the unit normal vector at x*. XFEM is very flexible in the sense that it allows cracks to be modeled along arbitrary paths without the need for remeshing. However, it is very difficult to obtain a converged solution, especially when pore fluid pressure is involved and there is intersection of cracks.

H(

Unfortunately, these two phenomena (fluid injection and the interaction of fractures) are the main focus of this work. Thus, XFEM is only used for the coupled thermos-mechanical simulations.

Cohesive element method

Cohesive elements are an alternative for modeling fracture propagation and the intersection of fractures. For the cohesive element method, the trajectory of the fracture to be modeled is typically defined a priori and cohesive elements are placed only along the fracture trajectory and between finite continuum elements. Cohesive elements can also be inserted between all finite elements of the mesh to model fractures propagating along unknown paths. Another method for unknown fracture paths is to combine cohesive element with remeshing. For each calculation step, the propagation direction is determined based on the stress and strain fields from the previous step. The mesh is then re-generated such that one edge of two continuum elements lies along the propagation path and a cohesive element is inserted between these two continuum elements.

A cohesive law, such as that presented in Section 2.3.2, is necessary to describe the mechanical response of the cohesive element. The cohesive law consists of a damage initiation criterion and a damage evolution criterion. The damage initiation can be based on either stress or strain. For instance, a criterion based on the maximum stress can be expressed as: For the damage evolution, to account for all fracture modes including mixed-mode fracture, an effective separation δ m is defined as:

2 22 12 m n s s        (2.92)
where δ n , δ s , δ t are the separations along the normal, first shear and second shear direction, respectively. Then the damage variable for a linear cohesive law is expressed as:

0 0 () () f m m m f m m m D          (2.93)
where f m  is the mixed-mode fracture separation, defined as:

2 f C m eqv G R   (2.94)
in which G C is the equivalent fracture energy and can be estimated using, for example, the BK law as presented in Section 2.3.1; R eqv is the equivalent strength, defined as:

2 2 2 12 eqv T s s R R R R    (2.95)
The actual stress components in the cohesive elements are then calculated via the damage variable D as:

(1

) p t D t  (2.96)
where t can be the normal stress σ n or the shear stresses σ s1 or σ s2 ; t p is the corresponding predictor which is estimated with the strain as if the material was undamaged (similar to equation Erreur ! Source du renvoi introuvable.). dispositif expérimental déployé in-situ est modélisée à l'aide de la méthode des éléments finis étendus [START_REF] Ngo | Numerical modeling of thermally-induced fractures in a large rock salt mass[END_REF]. Cette méthode permet une représentation de la géométrie de la fracture indépendante du maillage d'éléments finis, tout en permettant de prédire l'orientation de la propagation de la fracture en fonction de l'état de contrainte autour de l'extrémité de la fracture. Les résultats numériques montrent un bon accord avec les observations et les résultats expérimentaux.

Vient ensuite l'étude de propagation des fractures due à l'injection de fluide à partir d'un puits de forage. Le problème est simplifié et seul l'effet thermique est pris en compte. Les simulations montrent que les changements de température peuvent provoquer la fracturation de la roche, même dans des conditions de confinement élevé. Il a été démontré que la contrainte initiale avait un impact significatif sur la croissance des fractures d'origine thermiques.

INTRODUCTION

Thermal fracturing can play an important role in the development of geothermal reservoirs.

The injection of cold water into the geothermal reservoir during the stimulation and circulation phase will lead to significant temperature decreases in the rock. In response to these temperature decreases, the rock will contract but is refrained by the surrounding rock. Thus, tensile stresses are induced in the rock close to the wellbore. The thermally-induced tensile stresses in the rock, whose amplitude is proportional to the temperature gradient, may exceed the in-situ stresses of the reservoir and result in the formation of fractures. Thermal fractures can be initiated in the vicinity of a wellbore where the rock is in direct contact with the injected cold water and high temperature gradients exist. Secondary thermal fractures can also be initiated from the walls of a major fracture through which cold water is circulated during the heat extraction phase [START_REF] Bruel | Heat extraction modelling from forced fluid flow through stimulated fractured rock masses: application to the Rosemanowes hot dry rock reservoir[END_REF]. The thermal fractures initiated near the wellbore can contribute to a reduction of the the peak hydraulic pressure required to initiate fracture in the rock during the hydraulic stimulation [START_REF] Cha | Laboratory system for studying cryogenic thermal rock fracturing for well stimulation[END_REF][START_REF] Zhao | Experimental investigation on thermal cracking, permeability under HTHP and application for geothermal mining of HDR[END_REF][START_REF] Zhou | Experimental study on hydraulic fracturing of granite under thermal shock[END_REF].

While the secondary fractures initiated from the walls of the major hydraulic fractures may enhance the permeability of the reservoir, which in turn can contribute to facilitating the fluid circulation or increasing water loss from the reservoir [START_REF] Bruel | Impact of induced thermal stresses during circulation tests in an engineered fractured geothermal reservoir: example of the Soultz-sous-Forets European hot fractured rock geothermal project, Rhine Graben, France[END_REF][START_REF] Ghassemi | A Review of Some Rock Mechanics Issues in Geothermal Reservoir Development[END_REF][START_REF] Huang | Modeling 3D thermal fracture propagation by transient cooling using virtual multidimensional internal bonds[END_REF]. Thermal effects as a result thermal fracturing have been suggested as a means for the so-called "soft stimulation" of geothermal reservoirs.

The effect of temperature has been shown to "soften" the fracturing of rock (Zhang et al., 2018a;Zhang et al., 2018b;[START_REF] Zhou | Experimental study on hydraulic fracturing of granite under thermal shock[END_REF]. [START_REF] Zhou | Experimental study on hydraulic fracturing of granite under thermal shock[END_REF] conducted laboratory experiments in which rock specimens were heated before performing hydraulic fracturing, showed that the initiation pressure of hydraulic fracturing decreases gradually with an increase in the temperature of the specimens. [START_REF] Cha | Laboratory system for studying cryogenic thermal rock fracturing for well stimulation[END_REF] utilized liquid nitrogen at temperature of nearly -200 °C to create a thermal shock to stimulate rock specimens. The thermallystimulated specimens were then used in hydraulic stimulation tests. The results showed a significantly lower initiation pressure for the specimens that had been thermally stimulated compared to normal specimens. A series of experiments by [START_REF] Keshavarz | Damage and changes in mechanical properties of a Gabbro thermally loaded up to 1,000°C[END_REF] in which rock specimens were heated to high temperatures also showed a drastic decrease in rock strength. In all of these experiments, microcracking that originated from temperature changes was the main cause of the reduced initiation pressure of hydraulic fracturing. Crystalline rocks are composed of multiple minerals that have different thermal expansion characteristics. When subjected to temperature changes, either increasing or decreasing, these minerals deform differently, leading to microcracking at the mineral interface. This microcracking significantly weakens the rock and shifts the fracturing mode of the rock from brittle fracturing to continuous fracturing during fluid injection.

The effect of the temperature on the fracturing of rock is numerically investigated in this chapter and is presented in the next sections. Different numerical techniques are used, including the extended finite element method in combination with fracture mechanics and the cohesive element technique based on the cohesive material concept.

SIMULATION OF ROCK FRACTURE IN A THERMAL COOLING EXPERIMENT

In this section, the extended finite element method is used to model an experiment of cooling effect on the fracturing of rock salt performed in an underground gallery. The results of this analysis were published in a journal article [START_REF] Ngo | Numerical modeling of thermally-induced fractures in a large rock salt mass[END_REF].

The main outcomes of the experiment are summarized, followed by the numerical modeling.

The modeling of thermal fracturing of the rock salt consists of two sequential simulations. A heat transfer simulation is first carried out and then, the temperature field obtained from the heat transfer simulation is exported and used as the loading for the fracture propagation simulation.

Summary of the cooling experiment

The primary objective of this cooling experiment was to verify whether a temperature drop is capable of generating fractures in the walls of salt caverns. If so, the secondary objective was to quantify the characteristics of the fractures (e.g. depth, aperture) in order to assess the actual consequences for a gas storage facility during exploitation with a fast cycle of storage and withdrawal, pressure drawdown and associated temperature drops [START_REF] Hévin | Perform a Thermo-mechanical Test in a Salt Mine, as part of SMRI's Research Program on High Frequency Cycling of Salt Storage Caverns[END_REF].

Test setup:

The experiment was conducted in a 120-meter deep salt mine located in a thick layer of rock salt from the Keuper lithostratigraphic unit. In the main gallery (Figure 3 In order to relax the initial horizontal stresses in the salt block, two parallel slots were dug close to the niche wall (Figure 3.1a). Before starting the cooling test, these two slots (1.5 m in depth) were backfilled with salt powder that has a porosity of about 30%.

Testing procedure:

During the experiment, three cooling-warming cycles were performed using the refrigeration system and the four fans installed in the cooling chamber. Each cycle consisted of a 28 day cooling stage, during which air temperature inside the cooling chamber was decreased from 14.5 °C (initial temperature in the gallery) to around -9 °C, followed by a 28 day warming phase. A fourth cooling phase, which also lasted 28 days, followed the first 3 cycles, and the

Heat transfer simulation

The modeling was performed for the first cooling stage since the fracturing process due to 3.1.

For the simulation, the initial temperature in the whole domain was considered to be 14.5°C.

The temperature was held constant at 14.5°C at all boundary planes: x = -5.7 m, x = 6.3 m, y = 6 m, y = -11.5 m, and z = -6.5 m (see Figure 3.4 for the coordinate system).

Fracture propagation simulation

The fracture propagation simulation consists of two steps. The first step is to initialize the initial stress conditions in the rock mass surrounding the salt block where the cooling experiment was performed. The second step is to model the fracture propagation; the extended finite element method (XFEM) is used for this purpose. XFEM does not require pre-defining the fracture trajectory. Instead, the fracture trajectory is determined based on the stress condition at the fracture tip. Here the fracture propagation direction is supposed to be perpendicular to the maximum tangential stress at the fracture tip, while the fracture propagation criterion is based on a comparison between the strain energy release rate G and the fracture energy G C as presented in equation (2.44), which is repeated here for convenience:

C GG  (3.1)
The strain energy release rate G can be estimated using the Virtual Crack Closure Technique (VCCT) (ABAQUS, 2016), which is based on the assumption that the strain energy released when a crack is extended by a certain amount is the same as the energy required to close the crack by the same amount. Figure 3.9 illustrates the principle of VCCT. In Figure 3.9a, the fracture extends from node m to node n, accompanied by the release of a certain amount of strain energy. This amount of energy is supposed to be equal to the work necessary to close the fracture from node n to node m in Figure 3.9b.

For a pure Mode I fracture with the fracture tip located at node 2 (or node 5) as shown in Figure 3.10, the energy release rate when the fracture propagates from node 2 is calculated as follows:

,25 16

1 2 v I Fv G bd   (3.2)
where F v,25 is the nodal force at nodes 2 and 5; v 16 is the displacement between nodes 1 and 6;

d is the length of the element at the crack front; and b is the thickness of the element (not presented in Figure 3.10). The fracture will begin to propagate when G I = G IC where G IC is the Mode I fracture energy.

91

The mixed-mode fracture energy G C in equation (3.1) can be estimated using different laws presented in subsection 2.3.1.4. In this simulation the BK law is used and is repeated in the equation below:

  m II III C IC IIC IC I II III GG G G G G G G G         (3.3)
The

Initial stress initialization:

As explained above, the mechanical simulation of the fracture propagation was performed using the same sub-domain as in the heat transfer simulation. However, the initial stress state applied. The geostatic stress is the lithostatic stress and is supposed to be isotropic. The second step simulates the excavation of the main gallery and the niche. In the third step, the two slots are excavated. The density of the salt and marl were presented in Table 3.1. The density of the upper layer of sandstone is 2500 kg.m -3 . The thermo-mechanical properties of all 3 rocks are listed in Table 3.2. The stresses computed in the domain of interest will be extracted and used as initial stresses in the sub-model for fracture propagation using the submodeling technique in ABAQUS (ABAQUS, 2016). The distributions of stresses σ x and σ y on the floor of the cooling chamber at different times are shown in Figure 3.12. Note that the excavation of the niche has the same effect on both σ x and σ y as they both decrease from -2.85 MPa to about -2.0 MPa (compression is negative).

Conversely, the excavation of the 2 slots modifies σ x and σ y differently. As it can be seen in The cracks then propagated rapidly and coalesced into one crack across the chamber after 5h, as shown in Figure 3.13c. The location and the time at which the first cracks appeared are consistent with the observations made on site. Indeed, [START_REF] Hévin | Perform a Thermo-mechanical Test in a Salt Mine, as part of SMRI's Research Program on High Frequency Cycling of Salt Storage Caverns[END_REF] reported that a macro-crack was visible after about 6h of cooling at x = 1.5 m.

Stress relaxation was also noticed around the cracks at 5h (σ x on the left of Figure 3.13). After 5h, the macro-crack on the floor did not propagate further. The region of maximum tensile stress was shifted farther from the fan box, as can be seen in the contour plot of σ x in Figure 3.13c,d. At the end of the first cooling stage, the temperature distribution became more homogeneous and steady. This caused a decrease in the amplitude of both σ x and σ y at the end of the first cooling stage. For instance, the maximum σ x , which was about 2.4 MPa after 5h, reduced to close to 1.0 MPa at 28 days.

Although the crack stopped propagating on the floor of the cooling chamber after 5h, the temperature decreases caused the crack to propagate deeper into the salt block. Figure 3.14

shows the crack area on the plane y-z at different times. The crack's depth at 5h was approximately 0.3 m. As the cooling continued, the crack gradually penetrated deeper into the rock mass; after 1 day it was about 0.8 m deep. The downward crack penetration slowed and stabilized after 5 days. This evolutionary trend was also observed in the experiment [START_REF] Hévin | Perform a Thermo-mechanical Test in a Salt Mine, as part of SMRI's Research Program on High Frequency Cycling of Salt Storage Caverns[END_REF]. The maximum crack depth at 28 days was approximately 1.15 m, which is consistent with the depth observed on site, about 1.25 m.

The evolution over time of the crack depth and the crack aperture on the floor of the cooling chamber at y = 1.8 m are illustrated in Figure 3.15. Note that both the depth and aperture of the crack increased rapidly during the first few hours, after which the increase rate slowed down and stabilized after 5 days. The crack aperture reached a maximum of 0.65 mm after 2 days.

The computed crack aperture is of the same order as that observed experimentally, which was about 1 mm [START_REF] Hévin | Perform a Thermo-mechanical Test in a Salt Mine, as part of SMRI's Research Program on High Frequency Cycling of Salt Storage Caverns[END_REF]. The crack aperture then slightly decreased after 5 days; this crack evolution is consistent with that observed experimentally. The reason for this decrease is that the temperature distribution became more homogeneous after 5 days, causing the stresses (including σ x ) to relax and decrease the aperture.

THERMAL FRACTURE FROM A WELLBORE

As mentioned in the introduction to this chapter, the injection of cold water through a wellbore into a hot rock formation will cause a sudden temperature decreases in the vicinity of the wellbore, which in turn will lead to the contraction of the rock. Tensile stresses will be induced as the rock around the wellbore is refrained from contracting by the surrounding rock. As a result, fractures can be initiated if the induced tensile stresses exceed the in-situ stress and the tensile strength of the rock. This phenomenon will be studied in this section using 2D simulations.

Model setup

The fracture propagation from a wellbore due to temperature changes was studied using a 2D model presented in For the boundary conditions, the temperature of all outer edges of the 2D domain and of the wellbore is maintained at an initial temperature of 20 °C.

The physical and thermal properties necessary for the heat transfer simulation are summarized in Table 3.3. Fracture propagation simulation:

The temperature field obtained from the heat transfer simulation is used as the loading condition for the fracture propagation simulation. The fracture propagation is modeled using the cohesive material concepts in combination with cohesive elements. The plane deformation condition is assumed and the rock mass is supposed to be homogeneous and linear elastic. The temperature field obtained from the heat transfer is symmetric about the center of the wellbore because of the symmetry of the boundary conditions. Therefore, in this fracture propagation simulation only Mode I fractures that propagate radially from the wellbore are expected. For this reason, 24 layers of cohesive material are evenly placed around the wellbore as illustrated in Figure 3.17. The mechanical properties of the rock formation and the initial stresses are given in Table 3.4.

For the boundary conditions, the displacements of all rock surfaces are confined in the direction perpendicular to them. The initial temperature of the rock is 220 °C, which is similar to the initial temperature used in the heat transfer simulation. The temperature of the rock at time t > 0 is imported from the heat transfer simulation. The fracture propagation simulation is also run for 10 days, similar to the time duration of the heat transfer simulation.

Simulation results

Results of heat transfer simulation

The contour plots of temperature at different times are presented in Figure 3.18. It can be seen that the temperature field is symmetric about the center of the wellbore as predicted. The cold front spreads away from the wellbore as time increases.

Results of fracture propagation simulation

The fracture propagation due to temperature changes is presented in The fracture pattern at the end of 10 days for all cases is shown in Figure 3.22. In all cases except the case of S h = 36 MPa, multiple fractures are initiated at first, but only two fractures that are perpendicular to the minimum principal stress S h continue to propagate further from the wellbore. For the case of S h = 36 MPa, all the initiated fractures continue to propagate in a symmetrical manner without any preferred direction; this can be observed more clearly in Figure 3.23, which shows the damage variable at the end of 10 days. All the fractures propagate radially and the fracture length is more or less identical in all directions.

The time evolution of the length of the major fracture is shown in Figure 3.24 for all values of S h . In all cases the fracture length increases rapidly at first and then the rate of increase slows 106 down. This is because the temperature gradients are largest in the early hours and decrease as time goes on. Another remark is that as the minimum stress S h increases from 12 MPa to 36

MPa, the length of the major fracture decreases from 2.5 m to 0.3 m.

The fractures induced near the wellbore due to temperature changes resulted from injection of cold fluid may reduce the hydraulic impedance and wellhead pressure [START_REF] Bruel | Impact of induced thermal stresses during circulation tests in an engineered fractured geothermal reservoir: example of the Soultz-sous-Forets European hot fractured rock geothermal project, Rhine Graben, France[END_REF][START_REF] Murphy | Thermal stress cracking and the enhancement of heat extraction from fractured geothermal reservoirs[END_REF][START_REF] Swenson | The effects of thermal deformation on flow in a jointed geothermal reservoir[END_REF] but may also lead to more water loss due creation of new unwanted flow paths. Thermal fractures near the wellbore may also contribute to the instability of wellbore, rendering the wellbore inaccessible as observed at Soultz EGS reservoir [START_REF] Schill | Hydraulic performance history at the Soultz EGS reservoirs from stimulation and long-term circulation tests[END_REF]. increases. Thus, in order to reduce the peak hydraulic pressure, it may be necessary to increase the length of fractures near the wellbore. One way to achieve this is by decreasing the temperature of the injected fluid, which will result in higher temperature gradients and thus induce higher tensile stresses. Many researchers are experimentally investigating an even more "radical" approach by using liquefied nitrogen at temperatures of nearly -200 °C to create thermal fractures [START_REF] Cha | Laboratory system for studying cryogenic thermal rock fracturing for well stimulation[END_REF]Zhang et al., 2018a;Zhang et al., 2018b). It is expected that as the temperature of the injected fluid is lowered, the thermal fracture will propagate further. The question is: How effective is the lowering of the temperature of the injected fluid in increasing the length of the thermal fracture? To answer this question, the following simulation scenario was studied. Water at three different temperatures T 0 = 5 °C, 20 °C, and 35 °C is injected into a rock formation that has an initial temperature of 220 °C. The simulation procedure is similar to that described in subsections 3.3.1. The main points are summarized below:

 The modeling of thermal fracturing consists of two simulations: a heat transfer simulation and a fracture propagation simulation  For the heat transfer simulation, the temperature of the wellbore is assumed to be constant and equal to the temperature of the injected water T 0 . The material properties needed for the heat transfer simulation are taken from Table 3.3.

 The fracture propagation simulation uses the temperature field obtained from the heat transfer simulation as the thermal load. The mechanical properties and initial stresses needed for the fracture propagation simulation are taken from Table 3.4.

Simulation results:

The results of the heat transfer simulation for the case T 0 = 20 °C are presented in subsection advance. This is the advantage of using XFEM compared to the classic modeling approach using finite element method, which is not able to model the fracture opening and is therefore limited to identifying the locations of first fractures.

The injection of cold water into a rock mass through a wellbore was shown to lead to the propagation of fractures, even under high compressive stresses. Under anisotropic stress conditions, the cooling of the wellbore leads to the initiation of multiple fractures but only have both displacements and pore pressure as degree of freedom (DOF)), while the fracture is modeled by a thin layer of cohesive material that is discretized by coupled fluid pressure/deformation cohesive elements COH2D4P (i.e., elements with nodes at four corners that have both displacements and pore pressure as DOF and two mid-edge nodes that have only fluid pressure as DOF). The continuum elements and the cohesive elements are connected through shared nodes. Figure 4.2 illustrates the finite element mesh used.

The mesh is refined near the fracture plane where high gradients of stress and pore pressure are expected. In order to obtain accurate results when using cohesive elements, the traction in the cohesive zone ahead of the fracture tip must be properly represented. This requires the size of the cohesive elements to be small relative to the size of the cohesive zone. In other words, the cohesive zone must to be discretized into multiple cohesive elements. The length of the cohesive zone and the size of the cohesive element are defined as L cz and L e , respectively; the number of cohesive elements in the cohesive zone is:

cz e e L N L  (4.1)
The length of the cohesive zone L cz can be calculated from the properties of the cohesive material as [START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF]:

2 IC cz T EG LM R  (4.2)
where E is the Young's modulus of the rock mass; G IC and R T are the Mode I fracture energy and the tensile strength of the cohesive material, respectively; M is a constant. Several authors have proposed different values for M; in general, M can vary from 0.21 to 1 [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF][START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF][START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF].

According to [START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF], three cohesive elements in the cohesive zone are sufficient to obtain accurate results. In all the simulations presented here, this suggestion is followed and a conservative value of M = 0.21 is used. In most cases, the cohesive zone is discretized with up to ten cohesive elements.

Following the guidance above and with the selected material properties presented in the next section, the finite element mesh was generated. Elements of size 5 cm are used close to the fracture, while elements of up to 5 m are used near the boundaries. In total, the mesh has 7076 nodes and 6603 elements.

Material properties and injection rate

The material properties typical for granitic rock masses were chosen and are summarized in Table 4.2. The fluid is injected at a constant rate Q = 0.001 m 3 /s over 20 seconds. The leakoff coefficient is set to zero, which ensures a storage-dominated regime of propagation.

Initial and boundary conditions

The poroelastic medium is assumed to be initially fully saturated with zero initial stress in all directions. The initial pore pressure can have an arbitrary value, but is set to zero in this simulation.

The boundary conditions are illustrated in Figure 4.2b. Constant pore pressure is maintained at the boundaries, except for the left boundary where a no fluid flow condition is enforced, which implies that the left edge of the model is a symmetric boundary. A zero-displacement condition is imposed for the normal direction of all boundaries.

Convergence issues and viscous regularization

Modeling fracture propagation using cohesive elements involves the progressive stiffness degradation and a softening response of the cohesive material, which are both known to lead to convergence difficulties. A common technique to avoid these problems and aid convergence is the use of viscous regularization of the constitutive equations, which allows the tangent stiffness matrix of the softening material to be positive for sufficiently small time increments.

The regularization process involves the use of a viscous stiffness degradation variable, D v , which is defined as:

  1 vv D D D   (4.3)
where μ is the numerical viscosity (s -1 ) that represents the relaxation time of the viscous system and D is the actual damage variable evaluated in the backbone model (see equations

(2.93)). The damaged cohesive response is then given by:

(1 )

vp t D t  (4.4) 123
where t can be the normal or shear stresses of the damaged cohesive elements; t p is the corresponding predictor that is estimated with the strain as if the cohesive material is undamaged.

Using this viscous regularization technique, the stresses on the damaged cohesive elements are allowed to be outside the limits set by the cohesive traction-separation law. However, these stresses will converge to those in the case without viscous regularization when / t  ,

where t represents time.

Some of the energy of the system is dissipated due to the use of the viscous regularization. To ensure the accuracy of the solution results, the dissipated energy associated with the viscous regularization must be small compared to the damage energy of the cohesive elements. The final solution results may be altered if this criterion is not met. A numerical viscosity μ = 10 -4 s -1 is used in this simulation.

Results and discussions

The time evolution of the hydraulic aperture at the injection point, the injection pressure, and the fracture length is presented in To assess the effect of the viscous regularization on the simulation results, the above simulation was re-run with a numerical viscosity μ = 10 -2 s -1 . The results of the re-run simulation for damage energythe dissipated energy associated with viscous regularization, and the hydraulic aperture at the injection point are presented in Figure 4.9. It can be seen that the dissipated energy associated with viscous regularization is approximately 17% of the damage energy (Figure 4.9a). The simulated hydraulic aperture is much higher than the analytical value. This shows that an over-large value for the numerical viscosity can alter the simulation results. Thus, care must be taken when viscous regularization is used. Simulations should start with a small numerical viscosity, e.g., 10 -4 s -1 , and the ratio between the dissipated energy associated with viscous regularization and the damage energy should be relatively small. The material properties, initial conditions, and boundary conditions are identical to those previously used for the KGD model. The fluid is injected at constant rate of 0.001 m 3 /s over 40 s. With the selected parameters, the dimensionless number K, calculated according to equation (2.86), is equal to 3.5 at t = 0.01 s. This ensures that the fracture will propagate in the toughness-dominated regime for most of the injection duration as intended. In fact, the viscosity-dominated and transition regimes only last 0.01 s from the beginning of the injection, which is negligible compared to the total injection time of 40 s. To aid convergence, a numerical viscosity μ = 10 -4 s -1 is also used.

Results and discussions

Simulation results and analytical predictions for the hydraulic aperture at the injection point, the injection pressure, and the fracture radius are presented in The contour plots of pore pressure and the effective stress in direction 2 are presented in

Concluding remarks

In section 4.1 simulations of the two most used hydraulic fracture models are presented -the KGD and the penny-shaped models -using the cohesive element technique. Good agreement between the numerical simulations and analytical predictions was obtained for both models.

Several numerical aspects associated with using cohesive elements for modeling fracture propagation were examined. These include viscous regularization and finite element size.

When the viscous regularization technique is used to aid convergence, the dissipated energy associated with viscous regularization should be significantly smaller than the damage energy, while the size of the cohesive element should be small enough to properly characterize the stress in the cohesive zone.

PROPAGATION OF A HYDRAULIC FRACTURE AND SLIP OF A PRE-EXISTING

FAULT

The propagation of a hydraulic fracture in a low-permeability granitic rock mass and the slip of a pre-existing fault are studied in this section. The slip behavior of a fault that has been intersected by a hydraulic fracture is carefully investigated. This scenario was chosen to mimic the geological conditions that may be encountered in igneous faulted basements targeted for the development of deep geothermal reservoirs; faulted areas are known to be a favorable control for geothermal upwelling [START_REF] Faulds | Characterizing structural controls of geothermal reservoirs in the Great Basin, USA, and Western Turkey: Developing successful exploration strategies in extended terranes[END_REF][START_REF] Meixner | Structural control of geothermal reservoirs in extensional tectonic settings: An example from the Upper Rhine Graben[END_REF]. Some of the most important aspects that may influence fault slip are examined. These include the fault friction coefficient, the rate at which fluid is injected, and the orientation of the fault with respect to the initial stresses [START_REF] Aochi | Self-induced seismicity due to fluid circulation along faults[END_REF][START_REF] Meller | The significance of hydrothermal alteration zones for the mechanical behavior of a geothermal reservoir[END_REF][START_REF] Meller | Is there a link between mineralogy, petrophysics, and the hydraulic and seismic behaviors of the Soultz-sous-Forêts granite during stimulation? A review and reinterpretation of petro-hydromechanical data toward a Better Understanding of Induced Seismicity and Fluid Flow[END_REF].

Model setup

A 2D plane strain model is used, consisting of two potential hydraulic fractures HF1 and HF2

and one pre-existing fault F1 (Figure 4.16). This configuration is a simplification of a common pattern known as splay fracture, which results from a fluid driven nucleation and propagation process from a permeable fault [START_REF] Zhang | Fluid-driven nucleation and propagation of splay fractures from a permeable fault[END_REF]. The model is 120 m wide and 200 m high. The length of fracture HF1 and fault F1 are 10 m and 30 m, respectively. The initial stress state consists of major and minor effective far-field stresses S h and S H , which act in directions 1 and 2, respectively. It is worth noting that two fractures HF1 and HF2 are perpendicular to the minor initial stress S h and the fault F1 is oriented at an angle θ from the major initial stress S H (hereafter called the orientation angle). The orientation angle is subjected to parametric study later (i.e., θ is changed from 10° to 45°). with S H as it propagates away from the fault, much like tensile or wing cracks do in a strikeslip stress environment. This is the rationale for positioning fracture HF2 in alignement with the stress S H .

The material properties used in the modeling (i.e., rock mass, cohesive materials, and fracturing fluid), initial stresses, and initial pore pressure are summarized in Table 4.4.

The boundary conditions for the single fault model are illustrated in Figure 4.17. The displacements normal to all the outer boundaries are constrained. The pore pressure is kept constant at its initial value of 23.7 MPa for all outer boundaries, except for the left boundary where no fluid flow condition is imposed. The poroelastic rock mass is assumed to be initially fully saturated (saturation degree equals 1) and remains saturated during the injection time.

To generate the finite element mesh, the rock mass is discretized by linear plane strain coupled pore fluid pressure/deformation continuum elements CPE4RP, while the cohesive layers are discretized by coupled fluid pressure/deformation cohesive elements COH2D4P. At the intersection of two cohesive layers (e.g., HF1 with F1 or F1 with HF2), a mid-edge node with pore pressure as a degree of freedom is shared between those cohesive elements to support which are from 0.6 to 0.65 [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF]. However, in faults that are filled with clay minerals or quartz, which are the main product of the hydrothermal alteration of rock, the friction coefficient can be reduced to as low as 0.3 [START_REF] Morrow | Chapter 3 Frictional Strength and the Effective Pressure Law of Montmorillonite and lllite Clays[END_REF][START_REF] Zoback | The importance of slow slip on faults during hydraulic fracturing stimulation of shale gas reservoirs[END_REF].

Thus, μ f = 0.35 represents the lower range for the friction coefficient of faults. The results of the simulation are then presented and discussed in detail.

Parametric studies for the friction coefficient, the injection rate, and the orientation angle are presented in subsections 4.2.3, 4.2.4, and 4.2.5, respectively. The role of these parameters on the behavior of fault slip will then be discussed. MPa (see Table 4.4). The initial minor total stress S hT is S hT = S h + p 0 = 29 + 23.7= 52.7 MPa.

The tensile strength of rock mass is R T = 2 MPa. Thus, the summation of S hT and R T is 52.7 + 2 = 54.7 MPa, which is very close the stable fluid pressure 55 MPa.

When the fracture HF1 intersects the fault F1 at time 35.1 s, the injected fluid is suddenly driven into the fault, which leads to multiple events happening simultaneously: be stable (i.e., aseismic) or unstable (i.e., seismogenic). Usually a threshold of dynamic slip rate is used to distinguish these two slip regimes. This threshold can range from 5 mm/s to 0.1 m/s [START_REF] Dublanchet | Interactions and triggering in a 3-D rate-and-state asperity model[END_REF][START_REF] Gischig | Rupture propagation behavior and the largest possible earthquake induced by fluid injection into deep reservoirs[END_REF][START_REF] Mcclure | An investigation of stimulation mechanisms in Enhanced Geothermal Systems[END_REF]. When compared to these proposed thresholds of dynamic slip rates, the slip of fault F1 in this simulation is likely to be aseismic.

The contour plots of pore pressure at different times are presented in Figure 4.22. It can be seen that the fluid pressure inside the fracture HF1 and the fault F1 is nearly uniform and much higher than the pore pressure in the rock mass.

Effect of friction coefficient on fault slip

The influence of the friction coefficient on the slip of the fault F1 is investigated in this subsection. Three cases with different values of friction coefficient μ f were studied and are listed in Table 4.5. The case S2, μ f = 0.65 was chosen because it is typical of granitic rocks [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF], while case S3 μ f = 1.0, with a higher than normal friction coefficient, was

Simulation results and discussion:

The time evolution of fluid pressure (at the injection point and at the center of fault F1), accumulative slip and slip rate averaged over the length of the fault F1 is presented in Figure 4.23 and Figure 4.24 for cases S2 and S3, respectively. Similar patterns of evolution were observed in both cases:

1) The injection pressure increases immediately after injection commences and then decreases as the fracture HF1 starts to propagate

2) The fluid pressure in the fault F1 remains constant initially, decreases as the fracture HF1 approaches it, and then increases suddenly as HF1 intersects F1

3) The slip rate is very small for the majority of the injection time, including at the onset of HF1, while a high slip rate occurs only when the fracture HF1 intersects the fault F1.

Figure 4.25 shows the slip rate as function of the friction coefficient when the fracture HF1 is initiated and when it intersects the fault F1. It can be seen that the slip rate at the onset of the hydraulic fracture HF1 decreases as the friction coefficient increases, from 5.1 mm/s for case S1 (μ f = 0.35) to almost zero for case S3 (μ f = 1.0). This can be explained as follows. The fault F1 will slip at relatively high rates whenever the shear stress on the fault exceeds the shear strength of the fault. Otherwise, the fault will have elastic deformations, which results in lower slip rates. When the fault has a small friction coefficient (e.g., case S1 where μ f = 0.35), the shear strength of the fault is also small. Thus, the shear stress caused by the onset of HF1 can exceed the small shear strength of the fault, causing the fault to slip. As the friction coefficient increases, the shear strength also increases and becomes larger than the shear stress induced by the onset of HF1. This prevents the fault from slipping. In these cases the fault F1 only slips when it is intersected by HF1, which is when the fluid pressure in the fault increases suddenly, causing the shear strength to decrease significantly.

It is also noted from Figure 4.25 that the slip rate upon intersection of HF1 with F1 increases with the increase of friction coefficient, from 6.2 mm/s for case S1 to 19.2 mm/s for case S2

and to 23 mm/s for case S3. This can be interpreted as follows.

For a fault that has a small friction coefficient, the shear stress required to cause the fault to slip is also small. This means that the strain energy accumulated during the stress buildup process prior to the slip event, which is also the energy released during the slip event, is small.

Consequently, the slip rate is small when slip occurs. Also, when the friction coefficient is small, the fault may slip multiple times during injection, but at a relatively small slip rate. In other words, the accumulated strain energy is released more evenly. This phenomenon is observed for case S1 where the fault F1 slips twice during fluid injection and both times the slip rate is relatively small (5.1 mm/s and 6.2 mm/s).

For a fault that has a higher friction coefficient, more strain energy is accumulated prior to the slip event because the fault shear strength is higher. This accumulated strain energy is all released at one time when the fault is intersected by the hydraulic fracture. Consequently, the slip is more violent with higher slip rates. For instance, the slip rate of fault F1 when it is intersected by the hydraulic fracture HF1 is 19.2 mm/s for case S2 (μ f = 0.65 ) and 23 mm/s for case S3 (μ f = 1). These slip rates are significantly higher than the lower threshold of the dynamic slip rate, which is from 5 mm/s to 10 mm/s [START_REF] Dublanchet | Interactions and triggering in a 3-D rate-and-state asperity model[END_REF][START_REF] Mcclure | Investigation of injection-induced seismicity using a coupled fluid flow and rate/state friction model[END_REF]. Thus, these slip events are likely to be seismic.

Correlations between the friction coefficient of fractures/faults and their slip behavior were observed in lab experiments and in in-situ hydraulic stimulation tests. Through a series of lab tests on shale reservoir rocks, [START_REF] Kohli | Frictional properties of shale reservoir rocks[END_REF] found that: (1) the friction coefficient of a fault increases linearly with a decrease in the clay content;

(2) faults with low clay contents, i.e., high friction coefficients, exhibit a velocity-weakening frictional behavior and the fault slip is unstable, whereas faults with high clay contents, i.e., low friction coefficients,

Effect of injection rate on fault slip

The effect of the injection rate on the fault slip is investigated in this subsection. Four cases with different injection rates were studied, ranging from 0.25 L/s to 3 L/s per unit thickness.

The injection rates for these four cases, denoted by Q1 to Q4, are listed in Table 4.6. The fluid is injected over 100 s for cases Q1 to Q3. In case Q4, the higher injection rate (3 L/s) leads to fast propagation of the hydraulic fractures HF1 and HF2. Therefore, the injection is terminated after 53 s. The friction coefficient of the fault F1 was chosen as 0.65 for all cases; this friction 151 coefficient is typical of granitic rocks [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF]. The orientation angle of the fault F1 is θ = 22°. Other material properties were given in Table 4.4.

Table 4. [START_REF] Gu | Hydraulic Fracture Crossing Natural Fracture at Nonorthogonal Angles: A Criterion and Its Validation[END_REF] The injection rates for cases Q1 -Q4 The time evolution of the injection pressure and the fluid pressure at the center of the fault F1 is presented in Figure 4.27 for all four cases Q1 to Q4. It can be seen that as the injection rate increases, the fracture propagates faster. For instance, the intersection between fracture HF1

Case

and the fault F1 occurs after 77.5 s of injection when the injection rate is Q = 0.25 L/s. In case Q2 (Q = 0.5 L/s) the intersection occurs after 35.1 s, and is 18.5 s for case Q3 (Q = 1 L/s), and 8.1 s for case Q4 (Q = 3 L/s). The injection rate increases 12 times (3/0.25) between Q1 and Q4, while the propagation velocity of the fracture HF1 is increased by 77.5/8.1 = 9.6 times.

This shows that the injection rate plays an important role on the propagation velocity of the hydraulic fracture HF1, but is not the only factor. Poroelastic effects and permeability of the rock mass also affect the fracture propagation velocity. The breakdown pressure, i.e., the maximum injection pressure, also increases with the injection rate, from 65.8 MPa for case Q1 to 74.5 MPa for case Q4 (Figure 4.29). Again, this is reflective of the influence of poroelastic effects. If poroelastic effects were not considered, the breakdown pressure would be the same for any injection rate.

The simulated accumulative slip and slip rate of the fault F1 are presented in Figure 4.28. It can be seen that the slip rate is close to zero during most of the injection time. A high slip rate only occurs upon the intersection of HF1 with F1. The slip rate of the fault F1 when it is intersected by the fracture HF1 and the breakdown pressure as function of the injection rate are presented in Figure 4.29. It can be seen that upon intersection the slip rate increases with the injection rate. However, this effect of the injection rate on the slip rate is not significant as the slip rate only increases from 15.8 mm/s to 21.9 mm/s (less than 40% increase) when the injection rate is increased 12 times from 0.25 L/s to 3 L/s. The relatively small effect of the injection rate on the slip rate can be explained by the fluid pressure in the fault F1 as shown in Figure 4.27b. After intersection of HF1 with F1, the fluid pressure in the fault F1 is more or less similar for cases Q1, Q2, Q3, and Q4. The fluid pressure is the main factor that affects the shear strength, and thus the slip of the fault; small differences in the fluid pressure in the fault result in small differences in slip rate. Nonetheless, the slip rate in all four cases studied is higher than the 10 mm/s threshold of the dynamic slip rate, indicating that the slip is likely to be unstable (i.e., seismic). [START_REF] Aochi | Self-induced seismicity due to fluid circulation along faults[END_REF]) used a different model to study the fault slip;

the fault is modeled as permeable zone of finite thickness but the authors came to similar conclusion that the induced seismicity, which is reflected by the fault slip, is self-induced and less sensitive to the injection rate. Rather, the injection rate delays the fault slip as the pressure buildup process at the injection point takes more time when the injection rate is small.

The hydraulic aperture at the center of the fault F1 after injecting 25 L of fluid as function of the injection rate is shown in Figure 4.30. It should be noted that unlike the breakdown pressure and the maximum slip rate that both increase with the injection rate, the hydraulic aperture of the fault F1 after injection of 25 L of fluid is relatively similar for all cases. This indicates that more unstable fault slip may not necessarily mean more permeability enhancement.

The contour plots of pore pressure at different times for all four cases are presented in Figure 4.31. The effect of the injection rate on the velocity of fracture propagation is clearly illustrated. For instance, after 40 s injection, the fracture in case Q1 is only approximately 5 m long, while for case Q4 the fracture HF1 (10 m long) is completely generated and the fracture HF2 has begun to propagate. four cases, denoted by Z1 to Z4, are listed in Table 4.7. For all of these cases, the fluid is injected at 0.5 L/s in 100 s, while the friction coefficient of the fault F1 is μ f = 0.65. Other material properties are unchanged and were given in Table 4.4. Simulation results and discussion:

The time evolution of the accumulative slip and the slip rate of the fault F1 is presented in The fracture tip is at the same location, approximately 3.7 m from the injection point. This

shows that the fracture HF1 propagates at a similar velocity. However, it can also be seen in

Concluding remarks

The propagation of a hydraulic fracture and the slip of a pre-existing fault were studied. The effect of several controlling parameters on the fault slip was investigated. These include the friction coefficient of the fault, the injection rate, and the orientation of the fault with respect to the initial stresses. The friction coefficient and the fault orientation were found to have a strong influence on the fault slip: the fault slip rate increases significantly with an increase in the fault friction coefficient or the fault orientation angle. While the injection rate also affects the fault slip, its effect is of a lesser extent.

PROPAGATION OF MULTIPLE HYDRAULIC FRACTURES AND INTERACTION WITH EXISTING FAULTS

The propagation of multiple hydraulic fractures and their intersection with pre-existing faults are studied. Specifically, the effect of the fault friction coefficient and the injection rate on the fault slip is examined in more detail. The model setup is presented in the next subsection.

Then the results of the sensitivity studies for the fault friction coefficient and the injection rate are subsequently presented.

Model setup

The model studied is a 2D plane strain one that consists of four potential hydraulic fractures, denoted by HF1, HF2, HF3, and HF4, and three pre-existing faults, denoted by F1, F2, and F3 Four potential hydraulic fractures HF1 to HF4 are modeled by four layers of an intact cohesive material which have similar fracture mechanics properties (i.e., tensile strength and fracture energy). The three faults F1, F2, and F3 are modeled by three cohesive layers that are initially damaged. The shear strength of the faults is modeled using the Coulomb friction law with the same friction coefficient. The rock medium has low permeability and is considered to be linear poroelastic.

The stimulation of the system of fractures and faults is performed by injecting fluid at 0.5 L/s over 200 s through the left of the fracture HF1. When a hydraulic fracture intersects a preexisting fault, it is assumed that the hydraulic fracture diverts into the fault. Thus, after the stimulation is completed a step like fracture system is obtained.

The properties of the materials (i.e., rock mass, cohesive materials, and fracturing fluid), initial stresses, and initial pore pressure are similar to those used in the previous section 4.2 and presented again in Table 4.8 for convenience. For the finite element mesh, the rock mass is discretized by linear plane strain coupled pore fluid pressure/deformation continuum elements CPE4RP, while the cohesive layers are discretized by coupled fluid pressure/deformation cohesive elements COH2D4P. At the intersection of two cohesive layers, a mid-edge node that has the pore pressure as a degree of freedom is shared between those cohesive elements to support fluid flow continuity. The mesh is refined around the fractures and faults (smallest size is 0.05 m) and gradually becomes coarser towards the boundaries (largest size is 10 m). In total, the mesh has 21450 nodes and 20240 elements.

Sensitivity study of friction coefficient

The effect of the friction coefficient of the faults is studied in this subsection. All three faults F1, F2 and F3 have the same friction coefficient. However, the friction coefficient of these fault are changed to examine how it affects the slip of the faults. Three cases with different values of friction coefficient are studied and summarized in Table 4.9. The injection rate is fixed at 0.5 L/s for all of these cases. The decrease becomes smaller the farther the fault is from the injection point. For instance, the decrease in the hydraulic aperture at the injection point is approximately 0.33 mm when HF1 intersects F1, is 0.17 mm when HF2 intersects F2, and is 0.08 mm when HF3 intersects F3.

Case

Meanwhile, the hydraulic aperture of each fault increases suddenly when it is intersected by a hydraulic fracture and decreases upon subsequent intersections between the other faults and other hydraulic fractures (Figure 4.37b). For instance, the aperture of the fault F1 increases by nearly 0.44 mm when it is intersected by hydraulic fracture HF1, but decreases by 0.12 mm when the hydraulic fracture HF2 intersects fault F2, and decreases again by 0.06 mm when the hydraulic fracture HF3 intersects fault F3. The sudden increase in the hydraulic aperture of a fault when it is intersected by a hydraulic fracture is caused by the sudden increase in fluid pressure in the fault at the moment of intersection. Upon subsequent intersections between other hydraulic fractures and other faults, the injected fluid is driven into the additional faults, causing the aperture of the previously activated fault to decrease.

The contour plots of pore pressure at different times and the advancement of hydraulic fractures and pre-existing faults are presented in Figure 4.38. Again, it can be seen that the fluid pressure is almost constant along the hydraulic fractures and the pre-existing fractures, so that the whole structure extends at a pressure just below the minimum stress if the tensile strength is disregarded.

The slip rate of each fault, presented in Figure 4.39, remains small during most of the injection duration. High slip rates occur only when a hydraulic fracture intersects a fault. The slip rate of faults F1, F2, and F3 upon intersection with fractures HF1, HF2, and HF3 is 20.4 mm/s, 34.9 mm/s, and 6.1 mm/s, respectively. When compared to the threshold of the dynamic slip rate of 5 mm/s, all these intersection events are deemed to be seismic.

CHAPTER CONCLUSIONS

The propagation of hydraulic fractures and the slip of pre-existing faults in low-permeability poroelastic rock masses were studied. The hydraulic fractures are modeled using the cohesive zone model concept in combination with cohesive elements, while the shear strength of the pre-existing faults is modeled using the Coulomb friction law. It was found that the preexisting faults slip at higher rates once they are intersected by the hydraulic fractures. This is because during the intersection, the fluid pressure inside the fault increases suddenly, causing the shear strength of the fault to decrease significantly and thus allowing the fault to slip. In most cases studied, the fault slip rate upon intersection is higher than the threshold of the dynamic slip rate, which indicates that the slip event is likely to be seismic.

The effect of the most important parameters on the slip of the pre-existing faults was investigated, including the fault friction coefficient, the injection rate, and the fault orientation with respect to the initial stresses. Faults that are more critically oriented with respect to the initial stresses slip at higher rates when they intersect with hydraulic fractures. Meanwhile, increasing the injection rate also increases the fault slip rate. However, short faults were found to be more sensitive to an increase in the injection rate. The friction coefficient was also found rate. This suggests that it is possible achieve the goal of stimulating a rock reservoir (i.e., activating pre-existing faults and connecting them together) while minimizing the risk of inducing seismic fault slip by reducing the friction coefficient of the pre-existing faults before hydraulic stimulation begins. This could be achieved by appropriate geochemical treatments according to a pre-identification of key minerals within the faulted zones.

At present, only the fault slip rate is used to characterize the dynamic nature of a slip event.

However, it is worth noting that the fault slip rate does not provide direct information about the history of the induced seismic waves or the dominant frequency content, both of which are crucial in assessing human perception of the seismic waves and the damage potential to structures. In the next chapter, a modeling procedure to tackling these questions is proposed and a more in-depth analysis of fracture propagation/ fault slip-induced seismic waves is presented. This will provide a source model for induced micro earthquakes from small finite structures that experience shear and extend by developing wing cracks that grow in directions parallel to the most compressive stress [START_REF] Johnson | A source model for induced earthquakes at The Geysers geothermal reservoir[END_REF]. However, contrary to [START_REF] Johnson | A source model for induced earthquakes at The Geysers geothermal reservoir[END_REF], who built a dynamic source tensor using the appropriate Green's function, the same dynamic mechanical concepts will be used to simulate seismograms at any location outside the source region.

Chapter 5

Simulation of induced dynamic effects

INTRODUCTION TO INDUCED MICROSEISMICITY

Induced microseismicity refers to the generation of seismic waves due to the slip of preexisting faults or the propagation of brittle fractures [START_REF] Hardy | Acoustic Emission/Microseismic Activity: volume 1: Principles, Techniques and Geotechnical Applications[END_REF][START_REF] Swindlehurst | Acoustic emission -1 Introduction[END_REF].

When a brittle material is loaded to some extent, atomic bonds between the atoms are broken, allowing the fracture to form and propagate. The strain energy of the system is altered and some of the strain energy is converted into kinetic energy in the form of elastic waves, which can be transmitted through the medium and recorded by seismic sensors. A schematic representation of the generation of seismic waves from fracture propagation is shown in Figure 5.1. The generation of seismic waves due to fault slip is somewhat similar: when the stress on fault exceeds the fault strength, either due to the reduction of the fault shear strength or an increase of the shear stress, the fault slips and a part of the strain energy is released in the form of seismic waves [START_REF] Ellsworth | Injection-Induced Earthquakes[END_REF].

The fault slip that occurs during hydraulic stimulation of geothermal reservoirs is typically tied to the reduction of shear strength of the fault due to an increase in the fluid pressure within the fault. It should be noted, however, that the slip of faults can be seismic (i.e., generating detectable seismic waves) or aseismic (i.e., not generating detectable seismic waves). In in-situ experiments, direct measurement of the ground motions (e.g., particle velocity or acceleration) during the fault slip can distinguish whether the fault slip is seismic or aseismic. In numerical simulation studies, the slip rate of the fault is generally accepted as a parameter to assess the dynamic character of the slip event [START_REF] Cappa | On the relationship between fault permeability increases, induced Stress perturbation, and the growth of aseismic slip during fluid injection[END_REF][START_REF] Dublanchet | Interactions and triggering in a 3-D rate-and-state asperity model[END_REF]McClure andHorne, 2011, 2014).

Microseismic monitoring has been the only tool for assessing the development of enhanced geothermal reservoirs [START_REF] Baria | Phase 2B Final Report of the Camborne School of Mines Project (Part 5 Sesmics. Section 1 Microseismic results)[END_REF][START_REF] Baria | Microseismic monitoring of the world's largest potential HDR reservoir[END_REF][START_REF] Schill | Hydraulic performance history at the Soultz EGS reservoirs from stimulation and long-term circulation tests[END_REF]. The underlying idea behind the use of microseismic monitoring as a mapping tool is that the locations of the seismic events are considered to be a proxy for the presence of hydraulic fractures or the activation of pre-existing faults. Thus, the development of the seismic "cloud", also known as a swarm of events, provides information about the development and dimensions of the reservoir accessible to fluid. Microseismic source characteristics can also be inferred from the recorded waveforms [START_REF] Gibowicz | Source parameters of seismic events at the Underground Research Laboratory in Manitoba, Canada: Scaling relations for events with moment magnitude smaller than-2[END_REF][START_REF] Pearson | Parameters and a magnitude moment relationship from small earthquakes observed during hydraulic fracturing experiments in crystalline rocks[END_REF][START_REF] Wang | Relationships between moment magnitude and fault parameters: theoretical and semi-empirical relationships[END_REF]. In order to measure the "size" of a seismic event, various magnitude scales have been derived.

Classical magnitude scales, which are calculated from the displacement amplitude of ground motion, include the local scale M L for local seismic events and m B or M S for teleseismic events [START_REF] Gutenberg | Earthquake magnitude, intensity, energy, and acceleration: (Second paper)[END_REF][START_REF] Richter | An instrumental earthquake magnitude scale[END_REF]. In 1977, the moment magnitude M W was introduced by [START_REF] Kanamori | The energy release in great earthquakes[END_REF] and has since then been considered as the only well defined, non-saturating magnitude scale and estimator of seismic events [START_REF] Bormann | The moment magnitude Mw and the energy magnitude Me: common roots and differences[END_REF]. The moment magnitude is defined as:

  10 0 2 log 9.1 3 W MM  (5.1)
where M 0 is the seismic moment, which is related to the fundamental parameters of the faulting process including the average rigidity of the faulted rock μ, the rupture area S, and the average displacement on the fault D as:

0 M SD   (5.2)
It should be noted that the seismic moment is the change in strain energy of the rock mass before and after the fault slips (i.e., strain energy release). Equation (5.2) was proposed by [START_REF] Kanamori | The energy release in great earthquakes[END_REF] based on analysis of the stress around a shear fracture embedded in an elastic medium by [START_REF] Knopoff | Energy release in earthquakes[END_REF].

It is worth noting that not all of the strain energy release during fracture propagation or fault slip is converted into seismic waves. A large portion of the strain energy release goes into creating new surfaces and melting the rock. Thus, the use of the seismic moment (i.e., strain energy release) to estimate the seismic magnitude assumes that a proportion of the strain energy release is actually being released in the form of seismic waves. In fact, only a tiny fraction of the strain energy release, roughly about 5x10 -5 , is converted into seismic waves [START_REF] Bormann | The moment magnitude Mw and the energy magnitude Me: common roots and differences[END_REF][START_REF] Kanamori | The energy release in great earthquakes[END_REF].

Although the moment magnitude M W has some advantage over the classical magnitude scales, it is still a static measure of the seismic event, i.e., based on the difference in strain energy before and after the fault slips. Thus, M W does not provide direct information about the time history of the seismic waves, which is crucial in assessing the human perception of the seismic waves and the damage potential to structures.

The effect of the induced elastic waves on human perception and on structures on the ground surface can be quantified using the peak ground acceleration (PGA) or the peak ground velocity (PGV). By comparing horizontal peak ground motions to observed intensities for significant earthquakes in California, [START_REF] Wald | Relationships between peak ground acceleration, peak ground velocity, and modified Mercalli intensity in California[END_REF] proposed a regression relationship between PGA or PGV and the Modified Mercalli intensity. The Modified Mercalli intensity scale was introduced by [START_REF] Wood | Modified Mercalli intensity scale of 1931[END_REF] to quantify the shaking pattern and the potential damage of earthquakes. An abbreviated description of the levels of Modified Mercalli intensity and the corresponding PGA and PGV is given in Table 5.1.

MODELING OF INDUCED DYNAMIC EFFECTS

Introduction

Currently, three main numerical methods are used to model induced microseismicity: the finite element method in combination with damage mechanics, the discrete element method, and the combined finite-discrete element method.

In the finite element method, the fracturing process is modeled using damage mechanics [START_REF] Kaiser | Numerical simulation of damage accumulation and seismic energy release during brittle rock failure-Part II: Rib pillar collapse[END_REF][START_REF] Tang | Numerical simulation of progressive rock failure and associated seismicity[END_REF][START_REF] Tang | Numerical simulation of sumulative damage and seismic energy release during brittle rock failure-Part I: Fundamentals[END_REF]. The material adopts a damageable constitutive law and microseismicity is simulated based on the assumption that the number of induced seismic events is proportional to the number of elements that are damaged. The strain energy released in the damaged elements is considered analogous to seismic energy.

In discrete element models, the material is represented as an assemblage of particles that have the shape of a circular disc for 2D models or a sphere for 3D models [START_REF] Cundall | Numerical modelling of discontinua[END_REF][START_REF] Potyondy | A bonded-particle model for rock[END_REF]. These particles are bonded together and interact with each other at contact points where constitutive contact laws are required to describe their normal and shear interactions. When the shear or tensile stress exceeds the shear or tensile strength of a contact bond, the contact bond is broken, allowing a crack to initiate and propagate. The magnitude of an induced seismic event is calculated based on the change in kinetic energy of particles before and after one or several bonds are broken. The method was used by Hazzard andYoung (2000, 2004) to simulate triaxial compression tests and a mine-by tunnel experiment. The authors showed that the model was capable of reproducing the Gutenberg-Richter relationship between the magnitude and the number of induced seismic events with the numerical b-value situated within the range found in experiments. Note that the Gutenberg-Richter relationship is expressed as:

log N a bM  (5.3)
where M is a magnitude scale, N is the number of events having a magnitude larger or equal to M, a and b are two positive constants. The constant b (i.e., the b-value), typically close to 1, is a parameter describing the relative abundance of large to smaller events [START_REF] Okal | On the variation of b-values with earthquake size[END_REF].

In the combined finite-discrete element method, the intact material is typically assumed to be linear elastic and is modeled by finite elements. Cohesive elements are inserted between the finite elements to model the fracture initiation and propagation [START_REF] Munjiza | The Combined Finite-Discrete Element Method[END_REF]. Once the fracture has occurred, the blocks of material are created and treated as discrete blocks. The kinetic energy that is induced by the fracture propagation is considered to be the seismic energy. This method has been used by [START_REF] Lisjak | Numerical simulation of acoustic emission in brittle rocks by two-dimensional finite-discrete element analysis[END_REF] to simulate induced seismic effects in an unconfined compression test on granite. The authors reported that the magnitudes of the simulated events tended to display a power-law distribution, with b values in agreement with those reported in the literature for granitic rocks.

The three numerical methods briefly presented above have been successful to some extent in reproducing several features of the dynamic effects that are induced by fracture propagation (e.g., b-value in the range found in experiments). However, the following fundamental questions (FQ) have not been fully resolved in these simulations:

 FQ1: How are dynamic effects generated through fracture propagation even when the loading is quasi-static?

 FQ2: How much of the strain energy is converted into kinetic energy?

 FQ3: What are the most important parameters that influence the induced dynamic effects?

Some explanations can be found as answers to these questions. For question FQ1, the brittle nature of the fracture propagation can induce dynamic effects even when the loading rate is quasi-static. Many experiments have shown that microseismicity is induced by fracture propagation [START_REF] Keshavarz | Comparing the results of acoustic emission monitoring in Brazilian and uniaxial compression tests[END_REF][START_REF] Moradian | Detection of cracking levels in brittle rocks by parametric analysis of the acoustic emission signals[END_REF]. For question FQ2, [START_REF] Kanamori | The energy release in great earthquakes[END_REF] proposed that a fraction, 2x10 -5 , of the strain energy change is converted into the kinetic energy of induced elastic waves. For the question FQ3, several authors have reported some relationship between the velocity of fracture propagation or fault slip rate and the intensity of the induced dynamic effects [START_REF] Boudet | Experimental study of the instability of crack propagation in brittle materials[END_REF][START_REF] Fineberg | Instability in the propagation of fast cracks[END_REF]. These explanations are plausible but they still have not answered the above questions in a concrete and quantitative way. In this section, I will attempt to quantitatively resolve these questions.

Explanations using discrete models

From an energy balance point of view, the generation of elastic waves from the fracture of a material is a process of conversion of strain energy into kinetic energy. A physics-based approach that models the induced elastic waves must be able to model this energy conversion process. The only model (that I know of so far) that illustratively explains the strain energykinetic energy conversion was proposed by [START_REF] Pollock | Acoustic emission -2: Acoustic emission amplitudes[END_REF] and is presented in Figure 5.2. The model consists of a mass m suspended by two springs of stiffness K. The stiffness of the lower spring is assumed to instantaneously change by an amount δK. As the stiffness of the system suddenly changes, the mass will move toward and vibrate around a new equilibrium position.

The change in strain energy before and after the change in stiffness ΔE s and the peak kinetic energy of the mass E kin are, according to [START_REF] Pollock | Acoustic emission -2: Acoustic emission amplitudes[END_REF]: The model by Pollock shows that the amount of induced kinetic energy is proportional to the change in the stiffness of the system. For a continuum material, the stiffness change can originate from the propagation of fractures.

2 22 1 1 ( ) ; 88 s kin K E K x E x K       (5.
To more clearly illustrate the generation of dynamic effects from the damage regardless of the loading rate, I also use a mass-spring model as shown in Figure 5.3. The model consists of two masses M1 and M2, and two springs 1 and 2. Each spring is characterized by its stiffness and its tensile strength: K 1 and R 1 for spring 1, K 2 and R 2 for spring 2, with R 2 smaller than R 1 (R 2 < R 1 ). The system is loaded by controlling the displacement of the mass M2, as illustrated in  Strain energy of spring 1:

  22 22 11 1 cos(2 ) 24 p E K x K a t     (5.8)
where ω is the angular frequency; and a is the displacement of the mass M1 when spring 2 is broken and is equal to R 2 /K 1 . Thus the maximum kinetic energy E kmax of the mass M1 is:

2 2 2 max 1 1 11 22 k R E K a K 
(5.9)

The time evolution of the kinetic and the strain energy is presented in Figure 5.4. From equations (5.5) to (5.9) and Figure 5.4, the following remarks can be made:

(i) The kinetic energy of the mass M1 is proportional to the tensile strength of the spring 2 and is inversely proportional to the stiffness of the spring 1 (equation (5.9))

(ii) There is a conversion between kinetic and strain energy: when the kinetic energy is zero the strain energy is a maximum, and vice versa (Figure 5.4) (iii) The frequency of the kinetic energy and strain energy is double that of the frequency of the displacement or the velocity (equation (5.5) to (5.8) and Figure 5.4) In the above calculation the strain energy of spring 2 and the energy loss due to the breakdown of the spring 2 are not included. For a discrete system of springs and masses, these simplifications are acceptable. The calculation shows clearly that the intensity of the induced dynamic effects, which is reflected by the maximum kinetic energy, depends on two factors:

the strength and the stiffness of the system. Also, the kinetic energy of the mass M1 is shown to be independent of the rate of the applied displacement, which means that no matter how slowly the displacement is applied, the mass M1 will vibrate harmonically once spring 2 is broken and the kinetic energy of the mass M1 will remain the same. By using this simple spring-mass model, the questions FQ1 and FQ3 are clearly and quantitatively resolved.

Explanations using continuum model: analytical calculation

Here the cohesive zone model is used to model the fracture in a uniaxial tension test of a plate and show that elastic waves can be induced even when the loading rate is quasi-static. The plate studied is shown in Figure 5.5. The height of the plate is L. The plate consists of two parts; the first part is a zero-thickness layer of cohesive material whose cohesive constitutive law is characterized by tensile strength R T and fracture toughness G IC (Figure 5.6). The second part is made of a linear elastic material, which is characterized by Young's modulus E and Poisson's ratio ν. The loading is given by controlling the displacement rate of the top edge of the plate (Figure 5.5). As the applied displacement increases, the cohesive layer will be damaged and the separation rate of the cohesive layer is calculated. In this subsection, only the phase after damage initiation in the cohesive layer is considered, i.e., after the stress in the cohesive layer has reached the tensile strength and the two faces of the cohesive layer start to separate. It can also be seen that if the length of the plate is large enough, α can be larger than 1. In such cases c u becomes negative. However, as the applied displacement continues to increase after damage initiation of the cohesive layer has occurred, the separation of the cohesive layer will increase with time and thus c u must be positive. This leads to the realization that in cases where α > 1 the calculation above is no longer valid. In fact, in the above calculation it is assumed that when the stress in the cohesive layer is changed by an increment Δσ, the whole plate experiences the same stress increment Δσ. This is indeed the case as long as α < 1. When α > 1, if the stress in the cohesive layer is changed by Δσ, it will not affect the whole plate instantaneously. Instead, the stress change will "propagate" as elastic "stress" waves from regions near the cohesive layer to outlying regions without dependence on the applied displacement rate. The length L of the plate can be seen as a characteristic length of the material. A closer look at α shows that α is in fact the ratio between the change in strain energy before damage initiation, From equations (5.17) and (5.18), α is calculated as:

22.5 4.5 5

S F E E     (5.19)
Since α is larger than 1, elastic waves are expected to be induced after damage to the cohesive layer has occurred.

Simulation results:

A contour plot of the normal stress S22 in direction 2 during one period from the beginning of damage to the cohesive layer is presented in The time evolution of the different types of energy in the whole model (strain energy, kinetic energy, and fracture energy) is presented in Figure 5.11. The strain energy increases at first and then starts to decrease when cohesive damage occurs, while the kinetic energy remains close to zero before the cohesive damage and starts to increase when the cohesive damage occurs. The fact that the kinetic energy is close to zero before cohesive damage insures that the relatively high applied displacement rate used in the simulation has a negligible effect on the solution. After cohesive damage, both the strain energy and the kinetic energy oscillate harmonically with the same period of approximately 0.8 ms. These two types of energy are also converted back and forth into each other: when the kinetic energy is zero the strain energy is largest and vice versa. The maximum value of the kinetic energy after damage to the cohesive material is approximately 17.5 J. The energy dissipated due to the damage of the cohesive layer, i.e., the fracture energy, increases suddenly from 0 to 5 J upon the damage of the cohesive layer. This value of fracture energy is equal to the predicted value in equation (5.18).

It is worth noting that the oscillation period of the displacement and the velocity (T = 1.6 ms) is double that of the kinetic and strain energies (T = 0.8 ms). This observation, also obtained in discrete spring-mass model, is clearly illustrated in Figure 5.12.

Influence of the applied displacement rate

In the previous numerical simulation (subsection 5.2.4.1), the applied displacement rate is increased so that the loading duration can be reduced, which in turn reduces the computation time. To demonstrate that the increase of the applied displacement rate has a negligible effect on the final results, the simulation is re-run with a displacement rate of 0.9 mm/s, which is ten times smaller than in the previous case.

The simulation results are presented in Figure 5.13 for the time evolution of displacement and velocity of the point A and in Figure 5.14 for the time evolution of different types of energy.

The damage of the cohesive material occurs at around 77 ms. After that, trends for the time evolution are seen to be identical to the previous case. For instance, after damage to the cohesive material, the displacement and the velocity of the point A and kinetic and strain energies oscillate. The oscillation period of the kinetic and strain energies is approximately 0.8 s (Figure 5.14b), while the displacement and velocity of the point A oscillate with a period of approximately 1.6 s (Figure 5.13b). These periods are identical to those that were obtained in the previous case. The maximum velocity of the point A is 0.32 m/s and the maximum kinetic energy is 17.5 J, which are again identical to those obtained in the previous case. Thus, the simulation shows that the reduced applied displacement rate only expands the elastic Five cases with different tensile strengths are considered. These are listed in Table 5.3. Since the applied displacement rate has been proven to have negligible effect on the response of the system after the damage of the cohesive material, a displacement rate of 9 mm/s is used. For all of these five cases the loading rate is applied within 20 ms.

Simulation results:

The time evolution of the displacement and the velocity of the point A are presented in Figure 5.15 and Figure 5.16, respectively. As can be seen in Figure 5.15, for case 1 (α = 0.24) the velocity of the point A is virtually zero the whole time, both before and after damage to the cohesive material. As predicted by equation (5.16), the velocity of the point A after the damage of the cohesive material is:

11 9 11.8 mm/s 1 1 0.24 A uu      
This velocity is of the same magnitude as the applied displacement rate.

200

As α is increased from 2.4 in case 2 to 14.1 in case 5 (Figure 5.15), which corresponds to an increase in the tensile strength from 1.9 MPa to 4.6 MPa, the damage to the cohesive material, which corresponds to the moment when the velocity of the point A changes significantly, is delayed, as expected. The dynamic behavior of the system after cohesive damage, which is reflected by the time evolution of the point A in Figure 5.15, is similar for all cases: the velocity of the point A oscillates with similar period of approximately 1.6 ms. However, as α is increased from 2.4 to 14.1, the maximum velocity of the point A increases from 0.19 m/s (21 times the applied displacement rate) to 0.99 m/s (111 times the applied displacement rate).

The maximum velocity of the point A is plotted in Figure 5.17a as a function of the coefficient α. An increasing linear relationship between these two quantities, which is represented as a straight line passing through the origin, is observed. Figure 5.17b shows the maximum velocity as a function of the tensile strength. Again, it is noted that the maximum velocity tends to zero for a small tensile strength and increases rapidly as the tensile strength increases.

The above observations are consistent with the conclusions regarding the role of the ratio between the strain energy release and the fracture energy, which is represented through the coefficient α, that as the ratio between the strain energy release and the fracture energy increases the intensity of the induced elastic waves will also increase.

The above results can be interpolated to a Mode II fracture (shear-mode): as the shear strength of the Mode II fracture increases the intensity of the elastic waves that are induced when the fracture slips will also increase. And for a Mode II fracture and a given stress state, the shear strength of the fracture is solely proportional to the friction coefficient of the fracture. In other words, if the friction coefficient of the fracture can be reduced, which will lead to a reduction in the shear strength of the fracture, the intensity of the induced elastic waves will be reduced.

ESTIMATION OF ACCELERATION INDUCED BY FAULT SLIP

As mentioned previously, the fault slip rate is generally used as the main parameter to assess the dynamic character of the fault slip. If the slip rate is higher than a certain threshold, the slip event is considered to be seismic. Otherwise, the slip event is aseismic. The magnitude of the slip event can also be obtained by using the slip rate in combination with the seismic moment. are preformed to estimate the potential effects on humans and structural systems by using simplified 2D models.

Model setup

The model studied is similar to the model used in section 4.2, except that the dimensions of the model are now real dimensions, as shown in Figure 5.18. Specifically, the model is a 2D plane strain one that contains only one existing 30 m long fault F1, i.e., similar to the case given in section 4.2 of Chapter 4. The fault is located at 2.5 km from the ground level and is oriented at angle θ from direction 1. The material is assumed to be linear elastic and the whole domain is composed of only one material. This assumption is not altogether realistic as real formations can be made of multiple layers of different materials. However, the objective is only to demonstrate the methodology and the simulation procedure. Additionally, undrained conditions are also assumed since high slip rates occur only during a short time period. The so-called undrained Poisson's ratio is assumed to be ν u = 0.3 [START_REF] Jin | Fluid-induced seismicity and fully dynamic earthquake rupture in fractured poroelastic media : models, discretization and computation[END_REF]. The undrained Young's modulus can be found by using the condition that the shear modulus is similar in both the drained and undrained conditions:

1 1 u u EE      (5.20)
where E and ν are Young's modulus and Poisson's ratio of the rock mass in the drained condition. When E = 30 GPa and ν = 0.22, E u is 31.96 GPa. The saturated density of the rock mass is assumed to be 2700 kg/m 3 . The speed of dilatational waves and shear waves are calculated using equation (2.62). All the material parameters are listed in Table 5.4. For boundary conditions, the fault is assumed to be located in an infinite medium. The ground surface (edge EH in Figure 5.18) is the only free surface. The model used in this simulation has finite dimensions. Thus, to eliminate the reflection of waves from the outer boundaries of the finite model, "quiet" boundary conditions must be used. In ABAQUS, the quiet boundary is assigned using infinite elements CINPE4.

206

For the finite element mesh, the 2D domain is uniformly discretized with elements that have a characteristic size of 0.25 m.

We have the relationship between the wavelength λ, wave speed C p , and the wave frequency f as: where L e is the size of the finite element and N is the number of finite elements used to discretize one wavelength.

Combining equations (5.21) and (5.22) gives the following relationship: approximately 3193 Hz. It should be noted that seismic events in the field typically have frequencies from 1 to 5 kHz [START_REF] Hardy | Acoustic Emission/Microseismic Activity: volume 1: Principles, Techniques and Geotechnical Applications[END_REF]. Thus, the finite element mesh used is considered appropriate.

Dynamic simulations were run for all cases that have been studied in subsections 4.2.3 (Effect of friction coefficient on fault slip) and 4.2.4 (Effect of injection rate on fault slip), and are repeated again in Table 5.5 and Table 5.6 for convenience. The time histories of the displacements of all nodes on two surfaces of the pre-existing fault are prescribed using the results from the quasi-static simulations done in subsections 4.2.3 and 4.2.4. The simulation results for these cases are presented in the next sections. the second arrival is Δt ≈ 1.65 -1.15 = 0.5 s. This is consistent with the theoretical prediction, which is given by:

1 1 1 1 2500 0.49 2174 3767 SP t d s CC            
where d is the distance from the source to the observation point, which is roughly 2500 m. It can also be seen that point C is farther from the wave-emitting source (i.e., the pre-existing fault) than the point B and the point D is the most distant (see Figure 5.18). Therefore, the incident waves (both dilatational and shear waves) arrive at point B earlier than point C, and arrive last at point D.

The maximum acceleration in directions 1 and 2 at points B, C, and D are plotted against the distance from the points to the projection of the pre-existing fault on the top surface in Figure 5.25. Among these three points, the accelerations at B are highest in both directions with A 1max = 0.015 m/s 2 (maximum acceleration in direction 1) and A 2max = 0.007 m/s 2 (maximum acceleration in direction 1). The total acceleration at this point B is: Compared to intensity levels presented in Table 5.1, this peak acceleration (A = 0.17% g)

would not be felt by human beings and is not capable of causing any structural damage.

PGA as a function of friction coefficient:

The peak ground acceleration (PGA) in directions 1 and 2 and the total PGA at points B, C, and D for all three cases are presented in Table 5.7, Table 5.8, and Table 5.9, respectively.

These results are plotted against the friction coefficient of the pre-existing fault and are presented in Figure 5.26. In general, for all three points B, C, and D, the PGA increases with an increase in the friction coefficient.

For case S1 with a friction coefficient μ f = 0.35, the total PGA, measured in %g, is maximum at point C with a value of 0.167, which corresponds to the intensity level I (Table 5.1), i.e., the induced elastic waves may not be perceived by human beings and there is no potential for damage to structural systems. However, for the case S2 where friction coefficient μ f = 0.65, which is a normal value for a fault friction coefficient, the total PGA measured in %g at points B, C, and D is 0.44, 0.53, and 0.27, respectively. When compared to the Modified Mercalli intensity scale figures given in Table 5.1, these PGA values correspond to an intensity level II-III (0.17% g ≤ PGA ≤ 1.4% g), which again means that the induced elastic waves are not capable of causing any structural damage but can be "weakly" perceived by humans.

For case S3 with a higher-than-normal friction coefficient, the total PGA increases up to 1.68% g at point B, which corresponds to the intensity level IV. Human perception of the induced elastic waves is greater than case S2 (human perception of induced elastic waves from 214 "weak" in case S2 to "light" in case S3, Table 5.1). However, there is still no potential for structural damage. 5.10, Table 5.11, and Table 5.12, which show the PGA in all cases for points B, C, and D, respectively. These results are plotted against the injection rate in Figure 5.27. Note that for all four cases, the friction coefficient of the fault is 0.65.

As can be seen from Figure 5.27 the PGA increases with the increase of the injection rate.

However, the increase rate of the PGA is faster when the injection rate increases from 0.25 L/s to 1 L/s than when the injection increases from 1 L/s to 3 L/s.

The total PGA at all three points B, C, and D for all four cases is in the range from 0.17% g to 1.4% g, which means the induced elastic waves may be felt the human beings but only weakly and there is no damage potential for structural systems. Conclusions and perspectives

CONCLUSIONS

This research is a part of the major project GEOTREF whose aim is to improve the understanding of fractured geothermal reservoirs in order to reduce geothermal geological risk and to ensure sustainable exploitation of the reservoir during the production phase. The principal objective of this research is to build numerical models to simulate the fracture propagation and the slip of pre-existing faults during hydraulic stimulation of deep geothermal reservoirs and to propose a procedure to model the induced seismicity. The main results can be summarized in the following:

1. The effect of thermal cooling on fracture propagation:

Thermal cooling effects are shown to be able to cause fracture propagation. These thermal fractures can contribute to an increase in the reservoir permeability. Numerical models based on the extended finite element method or cohesive element concepts were performed and were

shown to be able to capture the main characteristics of thermal fracture propagation. The low thermal conductivity of the rock reduced heat transfer process; thus, the thermal cooling effect is more pronounced over the long term.

Modeling of the propagation of hydraulic fractures and the slip of pre-existing faults during the stimulation phase

The hydraulic fractures are modeled using the cohesive zone model concept, while the shear strength of the pre-existing faults is modeled using the Coulomb friction law. It was observed that the slip rate of the pre-existing faults increases suddenly when they are intersected by the newly created hydraulic fractures. The main reason for this is that during intersection the fluid pressure inside the faults increases suddenly due to fluid flow from the hydraulic fractures into the faults, causing a significant decrease in the shear strength of the faults, allowing the faults to slip at higher rate. In most cases studied, the fault slip rate upon intersection with the newly created hydraulic fractures is higher than the threshold of the dynamic slip rate, which indicates that the slip event is likely seismic.

The effect of the fault orientation, fault friction coefficient, and the injection rate on the fault slip

Through parametric studies, it was found that:

 Faults that are more critically oriented with respect to the initial stresses show higher slip rates when they are intersected with hydraulic fractures  Increasing the fluid injection rate also increases the fault slip rate. However, short faults are more sensitive to an increase in the injection rate than longer faults.

 The friction coefficient was also found to have a strong influence on the fault slip. As the fault friction coefficient decreases, the fault slip rate upon intersection with hydraulic fractures also decreases. If the friction coefficient is low enough, the fault slip rate can be reduced to levels below the threshold of the dynamic slip rate. This suggests that one can still achieve the goal of stimulating a rock reservoir (i.e., activation of the pre-existing faults and connecting them together) while minimizing the risk of inducing seismic fault slip by reducing the friction coefficient of the preexisting faults before the hydraulic stimulation begins.

The generation of elastic waves from fracture propagation and fault slip

 It was theoretically shown (see chapter 5) that the main mechanism leading to the generation of elastic waves when the material is damaged (e.g., fracture propagation or fault slip) is that, during the material damage process the amount of strain energy released is larger than the amount of energy that can be dissipated through the material damage (e.g., energy consumed to create new fracture surfaces). Thus, the excess strain energy is converted into kinetic energy in the form of elastic waves.

 The intensity of the induced elastic waves is dependent on the amount of the excess strain energy, or in other words, is dependent on the ratio between the strain energy release and the energy dissipated through damage.

 It is also shown that as long as the external loads do not introduce significant dynamic effects into the system (i.e., the external loads are still considered to be quasi-static), the dynamic response of the system that results from the induced material damage is independent of the loading rate.

The methodology for computing peak ground accelerations (PGAs) induced by fault slip

 A numerical modeling approach is proposed to compute the ground accelerations induced by fault slip. The sequential approach combines quasi-static coupled hydromechanical simulations with dynamic simulations.

 For the configurations studied in this research, the estimated PGAs may be felt by the public but will not cause damage to structural systems or the built environment. For instance, the maximum PGA for all studied cases is approximately 2 %g, which is smaller than the 4 %g limit specified by the Building code Eurocode 1998 above which seismic designs must be implemented (Clause 3.2.1 (5)P of the Eurocode 1998 -part 1).

RECOMMENDATION FOR FUTURE RESEARCH

 The main results of the numerical models (chapter 4), with respect to the relationship between the friction conditions of the faults, the injection rate, and fault orientation with its slip behavior, are qualitatively supported by in-situ observations. New  The behavior of fractured geothermal reservoirs during the stimulation phase is the main focus of this research. The reservoir behavior during the production phase still requires further research as the production phase is expected to last as long as 30 to 50 years. The models that have been built for the stimulation phase can be modified to study the behavior of the reservoirs during the production phase, e.g., adding physical couplings such as thermal processes.

 The thermal effects have been evaluated separately from fluid processes in chapter 3.

The cooling experiment using rock salt showed that a temperature decrease of 25 °C was able to create a 10 m long, 1 m deep, and 1 mm wide fracture. In deep geothermal reservoirs, temperature changes may be of the order of hundreds of °C. Therefore, it is expected that the thermal cooling effects in geothermal reservoirs will be more significant. The first well-known effect is the increase over time in the injectivity at the reinjection borehole, which is linked to the contraction of the rock mass. The second possible effect is the risk of thermal breakthrough, i.e., the injected fluid concentrates in a more direct flow path linking directly the two boreholes. The thermal breakthrough, if occurs, will be detrimental to the sustainability of any geothermal project. Several studies have also suggested that thermal contraction of the geothermal reservoirs is responsible for the settlements observed on the ground surface above the reservoir. Thus, the thermal cooling effects may require further research at a larger scale and over long time scales with respect to its impact on the reservoir permeability

Figure 1 . 1

 11 Figure 1.1 Estimates of the worldwide resource base for geothermal and other energy sources[START_REF] Brown | Mining the earth's heat: hot dry rock geothermal energy[END_REF]. Note: 1 quad = 293 x 10 9 kWh; world energy consumption in 2010 is 524 quads

Figure 1 . 2

 12 Figure 1.2 Original concept of an HDR geothermal system consisting of two wells and a vertical fracture created by hydraulic stimulation. Left: Front view showing typical dimensions for the wells, depth of the reservoir, and temperature. Right: 3D view showing the orientation of the fracture relative to the minimum stress. (Adapted from[START_REF] Brown | Mining the earth's heat: hot dry rock geothermal energy[END_REF].)

Figure 1 . 4

 14 Figure 1.4 Schematic of the EGS concept that contains one injection well and three production wells. The permeability of the rock mass is enhanced by the shearing and opening of existing joints/fractures. Source: DOE (2015).

Figure

  Figure 1.5.

Figure 1 . 5

 15 Figure 1.5 Elevation view of the Phase I (GT-2B and EE-1) and Phase II (EE-2 and EE-3) wellbore at the Fenton Hill HDR project[START_REF] Grigsby | Rock-water interactions in the Fenton Hill, new Mexico, hot dry rock geothermal systems I. fluid mixing and chemical geothermometry[END_REF]

Figure 1 . 6

 16 Figure 1.6 Front view of the seismic cloud of the EGS-system at Rosemanowes[START_REF] Jung | EGS -Goodbye or Back to the Future[END_REF]

  The two initial boreholes (GPK1, 3600 m deep and EPS1, 2200 m deep) were mainly for the purpose of investigating the nature and the structure of the rock mass. The other three boreholes (GPK2, GPK3 and GPK4) are located at depths 4400 to 5000 m. These boreholes form a triplet system in which GPK3 is used for injection of fluids, while GPK2 and GPK4 are dedicated to production. An elevation showing the relative locations of the four boreholes GPK1 to GPK4 is presented in Figure 1.7. All boreholes drilled had open sections of 500 -750 m in the lowest part that would be stimulated. Each borehole of the triplet system was stimulated by injecting large volumes of water at various flow rates. A volume of 22000 m 3 was injected at 50 L/s into borehole GPK2, which resulted in a wellhead pressure of 14.5 MPa at shut-in. Stimulation of the borehole GPK3 involved pumping of 34000 m 3 of water at a flow rate from 50 to 90

Figure 1 . 7

 17 Figure 1.7 (a) Schematic map of the Rhine graben and location of the Soultz-sous-Forêts EGS site (Calò et al., 2014), (b) Simplified deep geology, well trajectories and location of open sections and leakages in the Soultz-sous-Forêts reservoir[START_REF] Held | Economic evaluation of geothermal reservoir performance through modeling the complexity of the operating EGS in Soultz-sous-Forêts[END_REF]

Figure 1 . 8

 18 Figure1.8 Slip of fault due to fluid injection, propagation of induced seismic waves, and potential impact on surface structures and human perception. Adapted from[START_REF] Rutqvist | Modeling of induced seismicity and ground vibrations associated with geologic CO2 storage, and assessing their effects on surface structures and human perception[END_REF]

  (a) Propagation of hydraulic fractures (b) Intersection of hydraulic fractures with pre-existing faults

Figure 2 . 1

 21 Figure 2.1 Schematic representation of the hydraulic stimulation of geothermal reservoirs

  ); p and 0 p are the current and initial pore pressures, respectively;

Figure 2 . 3

 23 Figure2.3 Two distinct zones in the vicinity of a crack tip. The nucleation and growth of microseparations take place in the process zone, leading to the creation of voids. The propagation and coalescence of voids allow the crack to grow. After[START_REF] Broberg | Cracks and Fracture[END_REF]

  Figure 2.4, to illustrate the formation and propagation of a macrocrack. Several micro-cracks exist in the initial configuration (Figure 2.4-a). As the applied load is increased, more micro-cracks are nucleated in the process zone (Figure 2.4-b) and propagate (Figure 2.4-c). The coalescence of these micro-cracks allows the macro-crack to propagatee). A new process zone develops a head of the tip of the macro-crack in parallel with the propagation of the macro-crack. The size of the process zone depends on the nature of the material and its internal structure. For brittle rocks, the size of the process zone is relatively small in comparison with the dimensions of the cracked body and hence a large part of the body remains elastic.There are three basic modes of crack tip displacement, as illustrated in Figure2.5. Mode I is the opening mode where the principal load is applied normal to the crack plane and the crack tends to open. Mode II corresponds to the case where the crack faces tend to slide over each other due to in-plane shear. Mode III occurs when out-of-plane shear is the predominant applied load. A cracked body can be loaded in any one of these modes, or a be studied and modeled using either the theory of fracture mechanics or the cohesive zone model. The fundamental elements of these concepts and related issues are presented in the following sections.

Figure 2 . 4 Figure 2 . 5

 2425 Figure 2.4 Development of a process zones and extension of the macro-crack from the tip of a blunt notch due to applied load σ y . The applied load increases from drawings (a) to (e). (a) Existing microcracks. (b) Nucleation of micro-cracks. (c) Propagation of micro-cracks. (d) Coalescence of microcracks. (e) Extension of the macro-crack. Source: adapted from[START_REF] Atkinson | Fracture Mechanics of Rock[END_REF]

  .

Figure 2 . 7

 27 Figure 2.7 Definition of the coordinate axis ahead of a crack tip

  contributions of Mode I, II, and III, respectively.



  Minimum strain energy criterionThis criterion is based on the local density of the strain energy in the region around the crack tip. The crack is assumed to grow in a direction along which the strain energy density factor S reaches a minimum value. Thus the direction of crack growth is determined by:

Figure 2 . 8

 28 Figure 2.8 Arbitrary contour around the tip of a crack

Figure 2 . 9

 29 Figure 2.9 Schematic of cohesive zone

Figure 2 .

 2 Figure 2.10 Schematic of a linear cohesive law for Mode I fracture. The three parameters of the model are the tensile strength R T , the critical separation δ f , and the initial stiffness K. σ and δ are the current normal traction and current normal separation, respectively

Figure 2 . 11

 211 Figure 2.11 Different types of faulting. The vertical stress is the major principal stress in normal faulting, the minor principal stress in thrust faulting, and the intermediate principal stress in strike-slip faulting. Adapted from Jaeger et al. (2009)

Figure 2 .

 2 Figure 2.[START_REF] Bormann | The moment magnitude Mw and the energy magnitude Me: common roots and differences[END_REF] Force-displacement for a granite crack surface. T is the shear force, N is the normal force, μ and μ d are the static and dynamic friction coefficients[START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF]

Figure 2 .

 2 Figure 2.14 Geometry of three single fracture models: (a) KGD model, (b) PKN model, (c) Pennyshaped model (adapted from Savitski and Detournay (2002)). S h is the minimum compressive stress

  two coordinates that correspond to an energy dissipation mechanism and a component of global fluid balance, as Figure 2.15 illustrates.

Figure 2 .

 2 Figure 2.15 Diagram illustrating the four limiting propagation regimes of hydraulic fractures (after[START_REF] Bunger | Toughness-dominated hydraulic fracture with leak-off[END_REF]

Figure 2 .

 2 Figure 2.16 Potential scenarios of interaction between a hydraulic fracture and existing fracture

  led to a simple criterion for predicting fracture will propagate across a frictional interface orthogonal to the hydraulic fracture. The Renshaw and Pollard crossing criterion depends on the in-situ stresses, the tensile strength R T of the medium, and the friction coefficient μ of the interface, and is given

Figure 2 .

 2 Figure2.17 Plot of crossing criterion for different intersection angles between 15° and 90°. The region to the right of each curve represents the crossing condition for that intersection angle. Source:[START_REF] Gu | Hydraulic Fracture Crossing Natural Fracture at Nonorthogonal Angles: A Criterion and Its Validation[END_REF]

Figure 2 .

 2 Figure 2.18 Illustration of normal and tangential coordinates for a smooth crack (ABAQUS, 2016)

Figure 2 .

 2 Figure 2.19 A 4-node enriched element in three different configurations: undeformed before damage, undeformed after damage, and deformed after damage. Nodes numbered from 1 to 4 are real nodes and denoted by full circles. Nodes numbered from 1 p to 4 p are phantom nodes that are superposed on the real nodes and become active only after fracture begins. (After (Zielonka et al., 2014)) Undeformed configuration before damage

Figure 2 . 21

 221 Figure 2.20 Coupled pressure/displacement cohesive element. Corner nodes are denoted by black circles and have both displacement and fluid pressure degrees of freedom (DOF). Mid-edge nodes denoted by red triangles have only fluid pressure DOF (after (Zielonka et al., 2014)) Undeformed configuration

Figure 2 .RÉSUMÉ

 2 Figure 2.22 Intersection of two fractures (horizontal and vertical) modeled by two cohesive layers. The mid-edge node 100 is shared by all four cohesive elements 10, 20, 30 and 40.(modified from (ABAQUS, 2016))

  .1) a niche was purposely excavated to house the experiment. The dimensions of the niche are 17.5 m in length, 12 m in width, and 4.5 m in height. The salt block used for the cooling experiment was left unexcavated and has a thickness of 1.5 m (Figure3.1b). A section of the niche floor (dimensions 3.6 m x 3.6 m) was isolated by a chamber(Figure 3.1a, b and Figure 3.2a), in which the temperature can be regulated by a refrigeration system and fans.

Figure 3 . 1

 31 Figure 3.1 Geometry of the gallery and position of the cooling chamber (adapted from Hévin et al. (2016)). Dimensions are given in meters

Figure 3

 3 Figure 3.2 (a) View of the cooling chamber (Hévin et al., 2016); (b) Time evolution of air temperature inside the cooling chamber during the first cooling stage

Figure 3

 3 Figure 3.3 (a) Layout of thermocouples: thermocouples K and S measured temperatures on the floor; thermocouples T3, T4, T5 measured temperatures at depths 0.2 m, 0.4 m, and 0.8 from floor's surface, respectively. (b) View from inside the chamber with sensors on the floor and fans on the left[START_REF] Hévin | Perform a Thermo-mechanical Test in a Salt Mine, as part of SMRI's Research Program on High Frequency Cycling of Salt Storage Caverns[END_REF]

  Figure 3.4 Sub-domain used for heat transfer and fracture propagation simulations. Dimensions are given in meters[START_REF] Ngo | Numerical modeling of thermally-induced fractures in a large rock salt mass[END_REF]

Figure 3 . 5

 35 Figure 3.5 Convection zones for heat transfer simulation with different values of the heat transfer coefficient (measured in W.m -2 .C -1 ) as follows: h = 2 for the 2 slots, h = 20 for the floor under the fan box, h = 3.2x 2 -1.6x + 19 for the floor inside the cooling chamber, h = 0.2 for the insulated zone, and h = 10 for the uncovered zone

Figure 3 . 6 Figure 3 . 7 Figure 3 . 8

 363738 Figure3.[START_REF] Gu | Hydraulic Fracture Crossing Natural Fracture at Nonorthogonal Angles: A Criterion and Its Validation[END_REF] Comparison between simulation and measurement for temperatures along (a) the main profile and (b) the transverse profile. Legends followed by an "S" represent computed results while those without "S" represent the measurements. The cooling chamber is located in 0 ≤ x ≤ 3.6 m and 0 ≤ y ≤ 3.6 m.

Figure 3 . 9 Figure 3 .

 393 Figure 3.9 Extension and closure of Mode I fracture. The energy released when a crack extends by a distance δa is the same as the energy required to close the crack by a distance δa. Black dots represent the nodes of finite elements. Source: adapted from ABAQUS (2016)

  was first reproduced using a full model with dimensions of 180 m x 180 m x 180 m (Figure 3.11) and consisting of three layers, which are, from the top: a 40 m thick sandstone layer, an 80 m thick rock salt layer, and a 60 m thick marl layer. All lateral faces and the bottom face of the domain are constrained. The simulation has 3 steps: In the first step, the geostatic stress is

Figure 3 . 11

 311 Figure 3.11 Geometric domain used for reproduction of the initial stress state

Figure 3 .

 3 Figure3.12c, the stress σ x was more relaxed than σ y , decreasing from -2.0 MPa before excavation of the slots to about -1.4 MPa afterwards. This is due to the direction of the slots, which are normal to the x-axis.

Figure 3 .

 3 Figure 3.13 shows the stresses in the directions x and y (σ x on the left, σ y in the middle) and the temperature (on the right) on the floor of the cooling chamber at different times during the first cooling stage. Note that stress σ x inside the chamber, which is compressive prior the cooling, changes to tension 1h after the start of cooling. It is also greater than σ y and increases upon further cooling. At about 4h30, when σ x exceeds the tensile strength of the rock salt, two cracks are initiated in two locations at a distance of about 1.1 m from the fan box (i.e. x = 1.1 + 0.5 = 1.6 m, with the width of the fan box is 0.5 m), as shown in Figure 3.13b. The locations of the two cracks are immediately in front of the fans where temperature is the most reduced.

Figure 3 . 13 Figure 3 . 14

 313314 Figure 3.13 Numerical results (σ x on the left, σ y in the middle, temperature on the right) on the floor of the cooling chamber at different times during the cooling stage

  Figure 3.16. The 2D rock formation has dimensions of width 10 m and length 10 m. A wellbore of radius 0.15 m is located at the center of the 2D formation. Coldwater is injected into the rock formation through the wellbore and potential thermal fractures may initiate from the wellbore. The dimensions of the 2D formation were chosen such that they are still reasonably larger than the length of the expected fractures. The modeling consists of two simulations: a heat transfer simulation and a fracture propagation simulation.

Figure 3 .

 3 Figure 3.16 Geometry for studying fracture propagation around a wellbore

Figure 3 .

 3 Figure 3.17 Layout of cohesive layers around the wellbore and the boundary condition.

Figure 3 .

 3 Figure3.19 shows the temperature along a radius at different times. It should be noted that for one specific curve in Figure3.19 the temperature gradients, which are a measure of the slope of the curve, are largest near the wellbore and gradually decrease farther away from the wellbore. This indicates that the potential induced tensile stresses are also largest near the wellbore. It is also observed that the temperature gradients in the region near the wellbore decrease over time and the temperature field becomes more homogeneous.

Figure 3 .Figure 3 .Figure 3 .

 333 Figure 3.20 represents the time evolution of the temperature at four points A, B, C, and Dwhich are located at 0.25 m, 0.7 m, 1.5 m, and 3.7 m from the center of the wellbore, respectively. The temperature decreases sharply in the hours immediately after the beginning of the simulation, and continues to decrease subsequently but at a much slower rate. The temperature at points closer to the wellbore also decreases more than at points farther away from the wellbore.

Figure 3 . 21

 321 Figure 3.21 Contour plots of damage variable D at different times (when S h = 12 MPa). Fractured zones correspond to D = 1.

Figure 3 .Figure 3 . 24

 3324 Figure 3.22 Deformed shape (amplified 20 times) with fracture pattern at day 10 for different values of the minimum stress S h . The maximum stress is S H = 36 MPa

3. 3 Figure 3 . 25 Figure 3 . 26 0 Figure 3 . 27

 33253260327 Figure 3.25 Results of heat transfer simulation for different temperatures of the injected water. (a) and (b) are contour plots of temperature at day 10. (c) is the temperature along a radius at day 10

Figure 4 . 1

 41 Figure 4.1 KGD model: (a) 3D domain with a vertical fracture, (b) 2D plane strain model used for simulation. a is the initial fracture length and s is the curvilinear coordinate along the fracture.

Figure 4 . 2

 42 Figure 4.2 Finite element mesh for the KGD model: (a) Meshing scheme, (b) Actual finite element mesh with 7076 nodes and 6603 elements, and boundary conditions. Note: u 1 and u 2 are displacements in direction 1 and 2, respectively; p is the pore pressure

Biot

Figure 4 . 3 Figure 4 . 4 Figure 4 . 5

 434445 figures, the simulation results are compared with analytical solutions. Good agreement between the simulation results and the analytical solutions was found for all quantities. The final fracture length is approximately 20 m, which is less than one half of the width of the model as expected (the width of the model is W = 45 m).

Figure

  Figure 4.6 Simulated spatial distribution of pore pressure at different times

Figure 4 . 7

 47 Figure 4.7 Simulated spatial distribution of effective normal stress in direction 2 at different times (Deformation in direction 2 was scaled up 5000 times)

Figure 4 . 8

 48 Figure 4.8 Damage energy and viscous regularization dissipation (for numerical viscosity μ = 10 -4 s -1 )

Figure 4 . 9

 49 Figure 4.9 Simulation results with a numerical viscosity μ = 10 -2 s -1 : (a) Damage energy and dissipated energy associated with viscous regularization, (b) Hydraulic aperture at injection point

Figure 4 .

 4 Figure 4.10 2D axisymmetric model used for simulation of penny-shaped fracture model

  Figure 4.11, with good agreement for all three quantities. The final fracture radius is approximately 6.8 m, which is much smaller than the radius of the model (R = 45 m), which ensures that the model boundaries have a negligible effect on the simulation results.

Figure 4 .

 4 Figure 4.12 displays the variation of the fluid pressure in the fracture and the effective stress normal to the fracture plane along the fracture radius. As in the case of the KGD model, the fluid pressure is nearly uniform along the fracture radius due to the negligible viscous dissipation. The effective stress normal to the fracture plane is zero where the fracture has been created, increases in the cohesive zone and peaks at the tip of the fracture with a peak value approximately equal to the tensile strength of 1.25 MPa (Figure 4.12b).

Figure 4 .Figure 4 . 11 Figure 4 . 12

 4411412 Figure 4.13 and Figure 4.14, respectively. These figures again confirm the uniform distribution of the fluid pressure in the fracture and the stress concentration in the region near the tip of the fracture.

Figure 4 . 13

 413 Figure 4.13 Simulated spatial distribution of pore pressure at different times

Figure 4 .

 4 Figure 4.14 Simulated spatial distribution of effective stress in direction 2 at different times. Deformation in direction 2 was scaled up 4000 times.

Figure 4 .

 4 Figure 4.15 Time evolution of damage energy and viscous regularization dissipation

Figure 4 .

 4 Figure 4.16 Single fault model. The model contains two hydraulic fractures HF1 and HF2 and one pre-existing fault F1. The sketch is not to scale.

Figure 4 .

 4 Figure 4.17 Finite element mesh and boundary conditions with an orientation angle θ = 22 . u 1 and u 2 are displacements in directions 1 and 2, p is pore pressure. The mesh has 16469 nodes and 15380 elements

  results for case S1 with μ f = 0.35, Q = 0.5 L/s, and θ = 22° The time evolution of the injection pressure, the fluid pressure at the center of the fault F1, and the hydraulic aperture at the injection point and at the center of fault F1 is presented in Figure 4.18. It can be seen that the injection pressure increases sharply from 23.7 MPa to 66.3 MPa after injection begins. Due to this fluid pressure increase, the fracture HF1 is initiated, accompanied by a sudden decrease in the injection pressure. As HF1 propagates the injection pressure continues to decrease and stabilizes at around 55 MPa (Figure 4.18a), which is approximately equal to the initial minor total stress plus the tensile strength of rock mass. The initial minor effective stress S h and the initial pore pressure p 0 are S h = 29 MPa and p 0 = 23.7

  (i) the injection pressure decreases suddenly by approximately 1.5 MPa (Figure 4.18a), which leads to a decrease of the hydraulic aperture at the injection point by 0.5 mm (Figure 4.18b) (ii) the fluid pressure in the fault increases suddenly from the initial pressure to 53.3 MPa (Figure 4.18a), which leads to the opening of the fault F1. The hydraulic aperture at the center of fault F1 increases suddenly by roughly 0.3 mm from 0.4 mm to 0.7 mm (Figure 4.18b). After the intersection between HF1 and F1, the fluid pressures at the injection point and in the fault are identical. This is because the fluid has a very small dynamic viscosity (i.e., dynamic viscosity of the fluid is 0.001 Pa.s) and the hydraulic aperture of fracture HF1 and fault F1 are relatively large. In fact, the fluid pressure is virtually uniform along both the fracture HF1 and the fault F1, as Figure 4.19 illustrates. It should be noted that the pressure at the center of the fault F1 decreases slightly immediately before it intersects with the hydraulic fracture HF1 (Figure 4.18a). This pressure decrease is because the zone ahead of the tip of HF1 is in volumetric expansion due to the stress concentrations, as shown in Figure 4.20a. As a consequence, the fluid in the nearby regions is sucked into this expansion zone (Figure 4.20b), which in turn causes a decrease in the fluid pressure in the nearby regions, including within the fault F1. If the volumetric fluid flow rate into this expansion zone is smaller than the volumetric expansion rate, the fluid in this zone may even undergo cavitation; in such cases the capillary effect needs to be considered. The accumulative slip and the slip rate averaged over the length of the fault F1 are presented in Figure 4.21. The slip rate remains relatively close to zero during most of the injection time except in two instances. The first instance is at the onset of the hydraulic fracture HF1. The onset of HF1 is accompanied by a rapid decrease in the injection pressure. This induces stress disturbances in the whole model. If these stress disturbances exceed the shear strength of the fault F1, which is the case in this simulation, the fault will slip at a relatively high rate. The second instance is when HF1 intersects F1, during which the fault shear strength is lost by the rapid increase in fluid pressure within the fault. The slip rate is approximately 5.1 mm/s at the onset of HF1 and 6.2 mm/s when HF1 intersects F1. It should be noted that the fault slip can

Figure 4 . 18 Figure 4 .

 4184 Figure 4.18 Time evolution of (a) injection pressure and fluid pressure at the center of the fault F1 and (b) hydraulic aperture at the injection point and at the center of the fault F1

Figure 4 .Figure 4 .(

 44 Figure 4.20 (a) Strain in direction 2 and (b) Fluid velocity vector at 4.4 s after injection

Figure 4 .

 4 Figure 4.22 Contour plots of the pore pressure at different times

  the purpose of parametric study. The injection rate is Q = 0.5 L/s, fault orientation angle is θ = 22°. All the other material properties were given in Table4.4.

Figure 4 .Figure 4 . 24 Figure 4 .

 44244 Figure 4.23 Results for case S2 (friction coefficient μ f = 0.65): (a) time evolution of injection pressure and fluid pressure at the center of fault F1, (b) time evolution of accumulative slip and slip rate of fault F1

Figure 4 .

 4 Figure 4.29 Slip rate of fault F1 when intersected by hydraulic fracture HF1 and breakdown pressure as a function of the injection rate

Figure 4 .

 4 Figure 4.32. The slip rate of the fault F1 when it is intersected by the hydraulic fracture HF1, as a function of the fault orientation angle, is shown in Figure 4.33. It can be seen from Figure 4.32b and Figure 4.33 that the maximum slip rate when the hydraulic fracture HF1 intersects with the fault F1 increases as the orientation angle increases. For instance, the slip rate increases from 4.1 mm/s for θ = 10 to 46.6 mm/s for θ = 45. It can also be seen from Figure4.32b that the fracture propagation is faster as the orientation angle increases. For instance, the fracture HF1 intersects the fault F1 at a time of 50.9 s, 38.1 s, 35.4 s, and 27.1 s for θ = 10°, 22°, 35°, and 45°, respectively. This is counterintuitive because the orientation angle was expected to only affect the slip of the fault F1. The propagation velocity of the fracture HF1 is similar in all cases Z1 to Z4 as long as the fracture tip is far from the fault F1. This is illustrated on Figure4.34, which shows the effective stress normal to the plane of the fracture HF1 (called stress S22) along the fracture HF1 for cases Z1 to Z4 after 10 s of fluid injection.

Figure 4 .

 4 Figure 4.34 that the stress S22 near the intersection point between HF1 and F1 is disturbed differently for different orientation angles of the fault F1. S22 is reduced more (i.e., becoming

Figure 4 .Figure 4 . 33 Figure 4 . 34

 4433434 Figure 4.32 Time evolution of (a) the accumulative slip and (b) the slip rate of the fault F1 for cases Z1 to Z4 with different fault orientation angles ranging from 10° to 45°

Figure 4 .

 4 Figure 4.35 Contour plots of pore pressure at different times for 4 cases Z1 to Z4 with different fault orientation angles ranging from 10° to 45°

(

  Figure 4.36). The model is 120 m wide and 200 m high. The length of the fracture HF1 is 10 m, while HF2, HF3, and HF4 are 3 m long. All of the three faults are inclined at 22° with respect to the major principal stress S H and are 10 m long. The initial stress state consists of major and minor effective far-field stresses S h and S H that act in directions 1 and 2, respectively. All hydraulic fractures HF1 to HF4 are perpendicular to the minor initial stress S h .

Figure 4 .

 4 Figure 4.36 Multiple fault model. The model contains four hydraulic fractures HF1 to HF4 and three pre-existing faults F1, F2, and F3. S h and S H are the minor and major initial stresses. The unit of measurement is meters. The sketch is not to scale.

  Figure 4.37a that the injection pressure increases immediately after the start of fluid injection and then decreases as the fracture HF1 propagates. The fluid pressure in each of the three faults remains at the initial pore pressure before their intersection with the hydraulic fractures occurs and increases sharply upon the intersection (Figure 4.37a). For instance, the hydraulic

Figure 4 . 37

 437 Figure 4.37 Simulation results of case MS2 (friction coefficient μ f = 0.65): time evolution of (a) injection pressure and fluid pressure at the center of faults F1, F2, F3 and (b) hydraulic aperture at the injection point and at the center of faults F1, F2, F3

FigureFigure 4 . 39 Figure 4 . 41 Figure 4 .

 4394414 Figure 4.38 Contour plots of pore pressure at different times

  to have a strong influence on the fault slip. As the fault friction coefficient decreased, the fault slip rate upon intersection with hydraulic fracture also decreased. If the friction coefficient is low enough, the fault slip rate can be reduced to levels below the threshold of the dynamic slip

FurthermoreFigure 5

 5 Figure 5.1 (a) Illustration of generation of elastic waves from fracture propagation, (b) Example of seismic signal (e.g., acceleration) recorded by sensors

  4)where x is initial extension of the spring system due to the self weight of the mass m. From equation (5.4) the ratio between the kinetic energy and the change in strain energy can

Figure 5 . 2

 52 Figure 5.2 A spring-mass model. Source:[START_REF] Pollock | Acoustic emission -2: Acoustic emission amplitudes[END_REF]

Figure 5 . 3 .

 53 Figure 5.3. As the mass M2 is slowly pulled downwards, internal forces of the same magnitude are generated in both the springs. At some point the internal force in the spring 2 reaches its tensile strength R 2 , and spring 2 breaks. Assume that the damage of spring 2 is

Figure 5 . 3

 53 Figure 5.3 New spring-mass model. K1 and R1 are the stiffness and strength of spring 1, K2 and R2 are the stiffness and strength of spring 2.

Figure 5 . 4

 54 Figure 5.4 Time evolution of kinetic energy and strain energy. T is the period of the vibration of the mass M1

Figure 5 . 5 Figure 5 . 6 where

 5556 Figure 5.5 Plate model made of two materials. The cohesive layer has zero thickness. The sketch is not to scale

  , and the energy loss due fracture G IC .

Figure 5 . 7

 57 Figure 5.7 Dimensions of the plate, boundary conditions and loading

Figure 5 . 8

 58 Figure 5.8 Normal stress S22 in direction 2 during one period after the damage to the cohesive layer. The damage begins at time t = 8 ms

  evolution of displacement and velocity of the point A (Figure5.7) is presented in Figure5.10. It can be seen that the displacement increases linearly at first, and then starts to decrease as the cohesive layer is damaged. But instead of decreasing to zero, it oscillates around the zero-position. A similar evolution pattern is observed for the velocity: the velocity is relatively small before the damage of the cohesive layer, but starts to increase once the cohesive damage occurs and oscillates afterward. The maximum velocity of the point A is 0.32 m/s, which is 35.5 times the applied displacement rate. Both the displacement and the velocity oscillate with the same period of approximately 1.6 ms.

Figure 5 .Figure 5 . 11 Figure 5 . 12

 5511512 Figure 5.10 Time evolution of displacement and velocity of point A

  As can be seen in the above expression α can be varied by changing either the tensile strength R T , the length of the specimen L, the fracture toughness G IC , or the Young's modulus E. Here, the tensile strength is change. The change in tensile strength for a Mode I fracture is equivalent to the change in shear strength for a Mode II fracture. Thus, conclusions obtained in this simulation are also valid for Mode II fractures.

Figure 5 .Figure 5 . 16 Figure 5 . 17

 5516517 Figure 5.15 Time evolution of the velocity of the point A in direction 2 for the 5 cases studied

  However, the most important question that remains unanswered is: what is the effect of the induced seismic waves on human perception and structural systems at a specific location? In this section, dynamic simulations of the propagation of seismic waves induced by the fault slip

Figure 5 .

 5 Figure 5.18 Model for dynamic simulation

  accurately model the propagation of elastic waves, one wavelength should be discretized into a number of finite elements, i.e.:

  = 3992 m/s, L e = 0.25 m, and N = 5, f is 3193 Hz. This means that the finite element mesh is capable of accurately modeling the propagation of waves that have frequencies up to

  Detailed results for case S1 are presented in this subsection, while those for cases S2 and S3 are presented in Appendix 4.

Figure 5 .

 5 Figure 5.19 and Figure 5.20 show the contour plots of accelerations in directions 1 and 2 during the arrival of the dilatational waves to the top surface of the model. Note that only the upper part of the model containing the top surface is presented. The dilatational waves arrive at the top surface of the model at time t = 1.15 s (Figure 5.19a and Figure 5.20a). After that, surface waves are induced and propagate along the top surface (Figure 5.19b,c and Figure 5.20b,c); a part the incident dilatational waves are reflected off the top surface back into the model and interfere with the incident waves (Figure 5.19d and Figure 5.20d).

Figure 5 .

 5 Figure 5.21 and Figure 5.22 give the contour plots of accelerations in directions 1 and 2 during the arrival of the shear waves to the top surface of the model. The shear waves arrive at the top

  the gravitational acceleration which is approximately equal to 9.8 m/s 2 .

Figure 5 . 23 Figure 5 . 24 Figure 5 . 25

 523524525 Figure 5.23 Time evolution of acceleration in direction 1 at three points B, C, and D. The arrows indicate the arrivals of elastic waves: left arrow shows the arrival of dilatational waves and dilatationalwave-induced surface waves; right arrow shows the arrival of shear waves and shear-wave-induced surface waves

Figure 5 .

 5 Figure 5.26 PGA in directions 1 and 2 and total PGA at points B, C, and D as a function of the friction coefficient

  27 PGA in directions 1 and 2 and total PGA at points B, C, and D as function of the injection

  laboratory experiments under controlled conditions analogous to shear-faulting reservoirs need to be done in order to quantitatively validate and calibrate the numerical models. Using the methodological developments in chapter 5 it is possible to analyze the responses recorded by micro-seismic sensors used in laboratory experiments, which are usually in the form of electrical signals (piezoelectric sensors), and calibrate the waveforms with geometrical and mechanical parameters of the shearing and/or opening zone. A similar analysis can be carried out for in-situ situations to find out possible relationships between the responses at seismometers and the source mechanism and mechanical properties of the source.

Figure AN- 2 Figure AN- 6 Figure AN- 7 Figure AN- 8 Figure AN- 9 Figure AN- 11 Figure AN- 12 Figure AN- 13 Figure AN- 20

 2678911121320 Figure AN-2 Time evolution of accumulative slip and slip rate for faults F1, F2, and F3 (case MS1)

Table 1 .1

 1 Required properties for an EGS reservoir. Source:[START_REF] Rybach | The future of geothermal energy" and its challenges[END_REF]

	Property	Value
	Fluid production rate	50 -100 L/s
	Wellhead temperature	150 -200 °C
	Total effective heat exchange surface	> 2 x 10 6 m 2
	Rock volume	> 2 x 108 m 3
	Flow impedance	< 0.1 MPa/(L/s)
	Water loss	< 10%

Table 1 .

 1 2 Modified Mercalli intensity scale and corresponding peak ground acceleration and peak ground velocity. Source:[START_REF] Wald | Relationships between peak ground acceleration, peak ground velocity, and modified Mercalli intensity in California[END_REF],[START_REF] Wood | Modified Mercalli intensity scale of 1931[END_REF]

	Intensity Peak acceleration (%g)	Peak velocity (cm/s)	Perceived shaking	Potential damage
	I	< 0.17	< 0.1	Not felt	None
	II-III	0.17 -1.4	0.1 -1.1	Weak	None
	IV	1.4 -3.9	1.1 -3.4	Light	None
	V	3.9 -9.2	3.4 -8.1	Moderate	Very light
	VI	9.2 -18	8.1 -16	Strong	Light
	VII	18 -34	16 -31	Very strong	Moderate
	VIII	34 -65	31 -60	Severe	Moderate to heavy
	IX	65 -124	60 -116	Violent	Heavy
	X+	> 124	> 116	Extreme	Very heavy

  Fracture geometry and fluid flow in the fracture. w is the fracture aperture, q is the longitudinal flow, g is the leakoff flow rate

	Substituting the Poiseuille equation (2.18) into the continuity equation leads to the field
	equation that governs fluid flow in the fracture, as:
	3 12 ww g t s   f p s        	(2.20)
	Figure 2.2
				3
	k	t		12 w 	(2.19)

Figure 2.6 Plate with a crack of length 2a subjected to a remote stress σ

			IC G		2 	s	(2.34)
	The fracture stress, which is the stress at which the crack starts to propagate, can be obtained
	by combining equations (2.30), (2.32), and (2.34), as:
		f			2  E a 	s	(2.35)
	0 The fracture criterion, as in equation (2.35), is obtained through the introduction of the s dW dE d     (2.29) dA dA dA concept of the energy release rate (equation (2.27)), which is the net change in potential
	Substituting equations (2.27) and (2.28) into equation (2.29) gives the criterion for fracture energy of the whole system that accompanies an increment of crack extension. Thus, the
	propagation as: energy release rate characterizes the global behavior without regard for the local stress
	I GG  IC concentrations, which is in contrast to the stress analysis approach presented below.	(2.30)
	Griffith used the results from the stress analysis of Inglis and showed that for a plate
	containing a crack subjected to Mode I loading as shown Figure 2.6, the potential energy is
	given by:
	0    		22 aB E	(2.31)
	where 0  is the potential energy of the uncracked body, which is independent of the crack
	length; B is the plate thickness. Combining equation (2.27) and equation (2.31) leads to the
	expression of the energy release rate G I as:
	I G		2 E 	a	(2.32)
	Since the formation of one crack of length 2a creates two surfaces of the same area A = 2aB,
	the total surface energy s W is given by:
	4 ss W aB 	(2.33)
	where γ s is the unit surface energy, in J/m 2 .
	Combining equation (2.28) and equation (2.33) leads to the expression of the fracture energy
	G IC for the Mode I fracture as:

Table 2 .1

 2 Functions f ij (θ) for Mode I and II fracture[START_REF] Anderson | Fracture Mechanics: Fundamentals and Applications[END_REF]

  88) are the traditional nodal shape functions associated the finite element mesh, u I is the nodal displacement vector, H(x) is the Heaviside function, a I is the vector of nodal enriched degrees of freedom associated with discontinuous displacement across the crack surface,

	where	()
	/ S h b is the vector of nodal enriched degrees of freedom associated with crack tip I  Ratio S H singularity, and () Fx  are asymptotic crack tip functions that are defined in equation (2.89):
			 ( , )	sin , cos , sin sin , sin cos 2 2 2 2
				Coefficient of friction

I N x

Table 3 .3

 3 Physical and thermal properties necessary for the heat transfer simulation[START_REF] Ghassemi | Analysis of Fracture Propagation under Thermal Stress in Geothermal Reservoirs[END_REF]

	Property	Unit	Value
	Density ρ	kg/m 3	2650
	Heat capacity C p	J/kg °C	790
	Heat conductivity λ	W/m °C	10.7

Table 4 .2

 4 Material properties used for the KGD model

	Property	Value	Source / Comment
	Poroelastic medium
	Young's modulus
	Poisson's ratio
	Biot's coefficient

Table 4

 4

	.5 The friction coefficient values, cases S1, S2, and S3
	Case	Friction coefficient
	S1	0.35 (studied in section 4.2.2)
	S2	0.65
	S3	1.0

Table 4 . 7

 47 Fault orientation angle for the four cases studied

	Case	Orientation angle θ
	Z1	10°
	Z2	22° (Case S2 in Section 4.2.3)
	Z3	35°
	Z4	45°

Table 4 .9

 4 Friction coefficient for cases MS1, MS2, and MS3

RÉSUMÉ

  Ce chapitre est consacré à l'analyse et à la simulation de l'effet dynamique induit par la propagation de la fracture ou le glissement de la faille (c'est-à-dire les ondes élastiques induites). La toute première question que l'on se posera lorsqu'il s'agira de modéliser les ondes élastiques dynamiques induites par la propagation de fractures ou le glissement de failles est la 'inertie reste encore trop compliquée. Ainsi, une procédure de modélisation séquentielle est adoptée pour modéliser les ondes élastiques induites par le glissement de faille en combinant deux simulations distinctes. La première simulation est une simulation de diffusion contrainte-transitoire quasi-statique dans laquelle la propagation de la fracture et le glissement de faille sont modélisés. Cependant, les forces d'inertie sont considérées comme négligeables dans la première simulation. Cette première simulation a été réalisée au chapitre 4, sections 4.2. La deuxième simulation est une simulation dynamique tenant compte des forces d'inertie et visant à modéliser la propagation des ondes élastiques induites par le glissement de la faille. À cette fin, les historiques temporels des déplacements des deux faces de la faille, obtenus à partir de la première simulation, sont considérés comme des charges externes pour la deuxième simulation. Les accélérations calculées à la surface du sol sont ensuite utilisées pour évaluer l'intensité du séisme en relation avec la a perception humaine des ondes sismiques et le potentiel d'endommagement des infrastructures.

	suivante: comment les ondes élastiques dynamiques sont-elles induites par la propagation de
	fractures ou le glissement de failles, même lorsque les charges externes sont quasi-statiques?
	Cette question est plus difficile qu'il n'y paraît, même si la génération d'ondes élastiques
	dynamiques est un fait incontesté qui a été observé dans de nombreuses expériences quasi-
	statiques telles que l'essai de compression uniaxiale et l'essai brésilien. Pour apporter une
	réponse quantitative à ce phénomène, un modèle discret composé de ressorts et de masses
	élastiques est proposé, puis étendu à un modèle continu en utilisant le concept de matériau
	cohésif. La réponse à la question posée est progressivement dévoilée, de même que les
	principaux facteurs qui influent sur l'intensité des ondes élastiques induites.
	Ensuite, je vais essayer de réaliser des simulations dynamiques de la propagation des ondes
	sismiques induites par le glissement de faille. Même s'il est démontré que le concept de
	matériau cohésif est capable de modéliser le processus de conversion d'énergie (énergie de
	déformation -énergie cinétique) et donc la génération d'ondes élastiques à partir de
	dommages matériels, une simulation hydromécanique entièrement couplée prenant en compte
	les forces d

Table 5 .

 5 1 Modified Mercalli intensity scale and corresponding peak ground acceleration and peak ground velocity. Source:[START_REF] Wald | Relationships between peak ground acceleration, peak ground velocity, and modified Mercalli intensity in California[END_REF],[START_REF] Wood | Modified Mercalli intensity scale of 1931[END_REF]

	Intensity	Peak acceleration	Peak velocity	Perceived	Potential
		(%g)	(cm/s)	shaking	damage
	I	< 0.17	< 0.1	Not felt	None
	II-III	0.17 -1.4	0.1 -1.1	Weak	None
	IV	1.4 -3.9	1.1 -3.4	Light	None
	V	3.9 -9.2	3.4 -8.1	Moderate	Very light
	VI	9.2 -18	8.1 -16	Strong	Light
	VII	18 -34	16 -31	Very strong Moderate
	VIII	34 -65	31 -60	Severe	Moderate to heavy
	IX	65 -124	60 -116	Violent	Heavy
	X+	> 124	> 116	Extreme	Very heavy

Table 5 .3

 5 Tensile strength and the corresponding α for simulation cases 1 -5

	Case	Tensile strength R T (MPa) α
	1	0.6	0.24
	2	1.9	2.4
	3	2.6	4.5 (studied in subsection 5.2.4.1)
	4	3.5	8.2
	5	4.6	14.1

Table 5 .4

 5 Material parameters for the dynamic simulation

	Property	Value
	Undrained Young's modulus	E u = 31.96 GPa
	Undrained Poisson's ratio	ν u = 0.3
	Dilatational wave speed	C p = 3992 m/s
	Shear wave speed	C s = 2133 m/s
	Saturated density	ρ = 2600 kg/m 3

Table 5 .5

 5 Friction coefficient for cases S1, S2, and S3

	Case	Friction coefficient
	S1	0.35 (studied in section 4.2.2)
	S2	0.65
	S3	1.0

Table 5 . 6

 56 Injection rate for cases Q1 to Q4

Table 5 .

 5 7 PGA at point B for cases S1, S2, and S3

	Case	Friction coefficient	PGA in direction 1 (m/s 2 )	PGA in direction 2 (m/s 2 )	Total PGA (m/s 2 )	Total PGA (%g)
	S1	0.35	0.013	0.004	0.013	0.14
	S2	0.65	0.042	0.011	0.043	0.44
	S3	1	0.159	0.042	0.164	1.68

Table 5 .

 5 8 PGA at point C for cases S1, S2, and S3

	Case	Friction coefficient	PGA in direction 1 (m/s 2 )	PGA in direction 2 (m/s 2 )	Total PGA (m/s 2 )	Total PGA (%g)
	S1	0.35	0.015	0.007	0.016	0.17
	S2	0.65	0.046	0.025	0.052	0.53
	S3	1	0.084	0.045	0.095	0.97
			Table 5.9 PGA at point D for cases S1, S2, and S3
	Case	Friction coefficient	PGA in direction 1 (m/s 2 )	PGA in direction 2 (m/s 2 )	Total PGA (m/s 2 )	Total PGA (%g)
	S1	0.35	0.006	0.006	0.009	0.09
	S2	0.65	0.019	0.018	0.026	0.27
	S3	1	0.037	0.030	0.048	0.49

Table 5 .

 5 10 PGA at point B for cases Q1, Q2, Q3, and Q4

	Case	Injection rate (L/s)	PGA in direction 1 (m/s 2 )	PGA in direction 2 (m/s 2 )	Total PGA (m/s 2 )	Total PGA (%g)
	Q1	0.25	0.030	0.008	0.031	0.32
	Q2	0.55	0.042	0.011	0.043	0.44
	Q3	1	0.100	0.030	0.105	1.07
	Q4	3	0.156	0.044	0.162	1.65

Table 5 .

 5 11 PGA at point C for cases Q1, Q2, Q3, and Q4

	Case	Injection rate (L/s)	PGA in direction 1 (m/s 2 )	PGA in direction 2 (m/s 2 )	Total PGA (m/s 2 )	Total PGA (%g)
	Q1	0.25	0.037	0.022	0.043	0.44
	Q2	0.55	0.046	0.025	0.052	0.53
	Q3	1	0.067	0.035	0.076	0.77
	Q4	3	0.074	0.056	0.093	0.95

Table 5 .

 5 12 PGA at point D for cases Q1, Q2, Q3, and Q4

	Case	Injection rate (L/s)	PGA in direction 1 (m/s 2 )	PGA in direction 2 (m/s 2 )	Total PGA (m/s 2 )	Total PGA (%g)
	Q1	0.25	0.019	0.016	0.024	0.25
	Q2	0.55	0.019	0.018	0.026	0.27
	Q3	1	0.032	0.024	0.040	0.41
	Q4	3	0.054	0.026	0.060	0.61

Acknowledgement I am grateful to the many people who participated in different ways in the success of my doctoral research. First of all, I express my sincere gratitude to my advisors: Prof. Frédéric Pellet and Dr. Dominique Bruel, who entrusted me with this research topic. They gave me a lot of encouragement and guidance during my work at the Center for Geosciences of MINES ParisTech. Their advice was very valuable. I thank Prof. Patrick Selvadurai from McGill University and Prof. Jean Sulem from Ecole des Ponts ParisTech for their acceptance of being the reviewers of my research work. I greatly appreciated their advice and all the discussions I had with them. I thank Prof. Hervré Chauris from MINES ParisTech, who accepted to be president of the jury, as well as Prof. Eva Schill from University of Darmstadt and Dr. Michel Garcia from KIDOVA who accepted to participate in the evaluation of this work.

I would like also to thank Prof. Michel Tijani for all the helpful discussions and the tough questions he gave me as well.

My sincere thanks also go to my friends and colleagues: Zezhong, Hao, Xiangdong, Paule, Mejda, Sara, Hafsa, Aurélien, Qinglin and Shuaitao.

Finally, my last thoughts are with my family and all my loved ones for their unconditional support in difficult times. To all, thank you very much.

Fontainebleau, Juin 2019 NGO Dac Thuong

12

SUMMARY

Thermal effect on the fracture of rock is studied in this chapter. Numerical simulations are used to model the initiation and propagation of a fracture due to cycles of temperature changes. First, the fracturing of rock salt in an in-situ cooling experiment is modeled using the extended finite element method [START_REF] Ngo | Numerical modeling of thermally-induced fractures in a large rock salt mass[END_REF]. This method allows a representation of the fracture geometry independent of the finite element mesh while still being able to predict the orientation of fracture propagation based on the stress state around the fracture tip.

The numerical results show good agreement with the experimental data.

Second, the fracture propagation from a wellbore due to fluid injection is studied. The problem is simplified and only thermal effects are considered. The simulations show that temperature changes are able to cause the rock to fracture, even under high compressive confinement. The initial stress has been shown to have a significant effect on the growth of the thermal fractures.

SUMMARY

This chapter is dedicated to the study of the propagation of hydraulic fractures and the slip of pre-existing faults in low-permeability igneous rocks due to fluid injection. The problem involves temperature changes, fluid flow both in the porous space and the hydraulic fractures/faults, and deformation of the rock mass. Temperature changes may exert some effects but are not coupled to the analysis and are not considered here because of the complexity that would be added. To account for the thermal effect, a sequential simulation procedure as presented in Chapter 3, which is used as a first approximation. In this chapter, the focus is exclusively on the hydro-mechanical effects. All the numerical simulations in this chapter were done using the ABAQUS finite element package. The cohesive zone model in combination with cohesive elements is exclusively used to model the hydraulic fractures and the pre-existing faults. Although the XFEM is able to do the same job, it requires much more effort to obtain an accurate solution and convergence issues often arise.

Firstly, the propagation of a single hydraulic fracture in a low-permeability poroelastic medium is simulated. Both the KGD plane strain model and the axisymmetric penny-shaped model are studied. The simulation results are then compared with the theoretical solutions presented in subsection 2.5.2 of Chapter 2. The main objective of this section is to demonstrate the capability of the cohesive element technique and to gain insight into all numerical aspects (e.g, viscous regularization technique, finite element mesh creation, and numerical parameters for the cohesive element) that require attention in order to obtain an accurate solution to the coupled hydro-mechanical problem.

Secondly, the propagation of a hydraulic fracture in rock medium that contains a pre-existing fault is studied. The slip behavior of the fault once it is intersected by the hydraulic fracture is carefully investigated. This scenario is a simplified version of the geological conditions that might be encountered in igneous basements targeted for the development of deep geothermal reservoirs. Some of the most important aspects that could influence the fault slip are examined. These include the fault friction coefficient, the rate at which fluid is injected, and the orientation of the fault with respect to the initial stresses. Several general conclusions on the role of these aspects are drawn.

SIMULATION OF THE PROPAGATION OF A SINGLE HYDRAULIC FRACTURE

In this section, two single hydraulic fracture models, the KGD plane strain and the axisymmetric penny-shaped models, are studied. The numerical results are then compared with the analytical solutions presented in subsection 2.5.2 of Chapter 2. The cohesive zone model in combination with cohesive elements is used to model the fracture propagation. In order to make the numerical results comparable to the analytical solutions, the model dimensions and several material properties are chosen such that all the assumptions made to obtain the analytical solutions are satisfied. For instance, the dimensions of the model are much larger than the fracture dimensions (e.g., length and aperture) to limit the effect of model boundaries. The cohesive properties are selected such that the cohesive zone is small relative to the fracture length. This ensures that the fracture propagation modeled using the cohesive zone model is equivalent to that modeled in the analytical solutions using the theory of linear elastic fracture mechanics. The permeability is also defined to minimize the poroelastic effects since these effects are not considered in the analytical solutions. This section only analyzes and presents results for the toughness/storage-dominated regime of propagation (near vertex K in Figure 2.15). Models for the viscosity/storage-dominated regime (near vertex M) can be readily obtained based on the models for vertex K by adjusting the viscosity of the fluid.

Plane strain KGD model

Geometry and finite element mesh

The model studied consists of a vertical fracture that propagates in a poroelastic medium as illustrated in Figure 4.1a. As the fracture is under plane strain conditions, a 2D plane strain model, as shown Figure 4.1b, can be used for the simulation without compromising the accuracy.

The 2D domain has dimensions of width W and height B. These dimensions are chosen such that the maximum fracture length is expected to be less than one half of the width W. The numerical model also requires the presence of an initial crack in order to incorporate the physics for the fluid flow in the fracture. The model dimensions and the length of the initial crack are summarized in Table 4.1.

fluid flow continuity (see Section 2.7.2). The mesh is refined around the fractures and the fault (smallest size is 0.05 m) and gradually becomes coarser towards the boundaries (largest size is 10 m). A typical mesh, obtained with an orientation angle θ = 22°, is shown in Figure 4.17.

The mesh has 16469 nodes and 15380 elements. Injection rate Q = 0.5 L/s per unit thickness

Subjected to sensitivity study show a velocity-strengthening frictional behavior and the fault slip is stable. An experimental study by [START_REF] Ikari | On the relation between fault strength and frictional stability[END_REF] also showed that weak faults (faults with a friction coefficient μ f < 0.5) only exhibit increased friction with slip rate, which suppresses instable slip, while faults with higher friction coefficients can exhibit unstable slip behavior. Similar observations were obtained in 2000 during the hydraulic stimulation of two wellbores GKP1 and GKP3 at Soultz-sous-Forêts. Based on an analysis of the induced seismic events, [START_REF] Meller | The significance of hydrothermal alteration zones for the mechanical behavior of a geothermal reservoir[END_REF] showed that large-magnitude seismic events occurred on fractures that have no significant clay filling, i.e., fractures with high friction coefficients. Meanwhile, clay-rich fractures, i.e., fractures with lower friction coefficients, tend to produce smaller seismic events. They also concluded that clay inside the fractures reduces the friction on fractures, thus preventing large stress buildup and promoting aseismic slip.

The time evolution of the hydraulic aperture at the center of the fault F1 is presented in Figure 4.26 for all three cases S1, S2, and S3. It can be seen that in all cases the hydraulic aperture increases by approximately the same amount when intersection between HF1 and F1 occurs, even though the fault F1 slips at a much higher rate in cases 2 and 3 than in case 1.

(a) (b)

Effect of fault orientation on fault slip

The effect of the fault orientation with respect to initial stresses is studied in this subsection.

Four cases with a different orientation angle θ ranging from 10° to 45° are examined. These

Sensitivity study of injection rate

The effect of the injection rate is studied in this subsection. The injection rate is varied from 0.5 L/s to 2 L/s and the slip rate of the faults is examined. Three cases with different injection rates are listed in Table 4.10. The friction coefficient of all of the three faults is 0.65. In conclusion, if the strain energy that is released during the fracture propagation is larger than the energy needed to create new fracture surfaces (i.e., the fracture energy), elastic waves will be generated regardless of the loading rate. The intensity of the induced elastic waves is proportional to the ratio between the strain energy release and the fracture energy.

Explanations using a continuum model: numerical simulation

The model presented in subsection 5.2.3 is numerically investigated. Again, the objective is to model the elastic waves induced by the fracture and investigate the factors that influence the intensity of the induced elastic waves. The following points will be demonstrated using the numerical simulations:

(i) Elastic waves are induced if α, i.e., the ratio between the change in strain energy E S before and after fracture occurs and the energy loss due to the fracture E F , is larger than 1 (subsection 5.2.4.1)

(ii) The loading rate does not affect the induced elastic waves (subsection 5.2.4.2)

(iii) The intensity of the induced elastic waves also increases if α increases (subsection 5.2.4.3)

Reference simulation when α > 1

The plate studied has a width W = 0.2 m and length L = 1 m. The top edge of the plate is pulled upwards at a constant and relatively small rate such that the loading does not introduce any dynamic effects. In this simulation, a displacement rate of 9 mm/s, which corresponds to a strain rate 0.009 s -1 , is chosen. This loading rate may seem too high, but as shown later this displacement rate actually does not cause any significant dynamic effect on the system.

The three other edges of the plate are fixed. An illustration of the model with boundary conditions and loading is given in Figure 5.7. With these boundary conditions and loading, the problem is equivalent to a 1D problem and only the displacement, velocity, and acceleration in direction 2 are non-zero. Thus, in the following the words displacement, velocity, and acceleration are used to describe the displacement, velocity, and acceleration in direction 2.

Parameters for the cohesive material and the linear elastic material are presented in Table 5.2.

deformation phase and delays the damage of the cohesive layer. After the cohesive layer is damaged, the behavior of the model is identical to the previous case.

Influence of tensile strength on the intensity of induced elastic waves

In this section, the influence of the ratio α on the intensity of the induced elastic waves, which is reflected by the maximum velocity of the point A, is investigated. The expression for α, as presented in equation (5.15), is repeated here for convenience:

CHAPTER CONCLUSIONS

In this chapter 5, the generation of elastic waves from the material damage is studied using a discrete model (mass-spring model) and a cohesive material model. It is shown that when a material is damaged the main reason leading to the generation of elastic waves is the excess strain energy that is released during material damage, which is greater than the energy that can be dissipated through the material damage (e.g., energy consumed to create new surfaces). The intensity of the induced elastic waves is dependent on the amount of the excess strain energy, or, put another way, is dependent on the ratio between the strain energy release and the energy dissipated through damage. It is also shown that as long as the loading does not introduce significant dynamic effects to the system (i.e., the loading is still considered as quasi-static), the dynamic response of the system that is induced by the material damage is independent of the loading rate.

The propagation of elastic waves induced by the fault slip is simulated by a modeling procedure. The effect of the induced elastic waves on human perception and structural systems is quantified using the peak ground acceleration. The simulation results show that lower fault friction coefficients or lower injection rates result in smaller peak ground accelerations. For the studied configuration, it has been shown that the induced elastic waves may be "weakly"

or "lightly" felt by human beings but are not capable of posing any concern to the structural systems on the ground.

Appendix 1: List of publications

ABSTRACT

The development and the exploitation of deep geothermal reservoirs are usually accompanied with induced seismicityan unwanted side effect. This research is focused on using numerical simulations to investigate the propagation of hydraulic fractures and the reactivation of pre-existing faults during the hydraulic stimulation of the reservoirs in an effort to better understand the fractured reservoir behavior and to reduce the potential risk of induced seismicity.

The induced seismicity is studied first from the standpoint of using the law of energy conservation in order to explain the mechanism of generating elastic waves from rock failure. Then an approximate approach is proposed to calculate the peak ground accelerations (PGAs) that are induced by the fault slip. The computed PGAs on ground surface are used to assess the human perception of the seismic waves and the damage potential to structures.

MOTS CLÉS

Géothermie, stimulation hydraulique, modélisation numérique, sismicité induite, mécanique de la rupture, FEM RÉSUMÉ Le développement et l'exploitation de réservoirs géothermiques profonds s'accompagnent généralement d'une sismicité induite -un effet secondaire indésirable. Cette recherche est axée sur l'utilisation de simulations numériques pour étudier la propagation des fractures hydrauliques et la réactivation de failles préexistantes lors de la stimulation hydraulique des réservoirs afin de mieux comprendre le comportement du réservoir fracturé et de réduire le risque potentiel de sismicité induite.

La sismicité induite est d'abord étudiée du point de vue de l'utilisation de la loi de conservation de l'énergie afin d'expliquer le mécanisme de génération d'ondes élastiques à partir d'une rupture de roche. Ensuite, une approche approximative est proposée pour calculer les accélérations de pointe de pointe (PGA) induites par le glissement de faille.

Les PGA calculés à la surface du sol servent à évaluer la perception humaine des ondes sismiques et le potentiel de dégradation des structures.

KEYWORDS

Geothermal energy, hydraulic stimulation, numerical modeling, induced seismicity, fracture mechanics, FEM