Qinglin Deng

Jérôme, Nolween, Quentin, Bérénice Catalina

Kamel Emmanuel

Xin

Lei, Chong, Haitao, Xiang, Yuzhou, Ke Xuan Jiang

Multi-scale hydro-mechanical behavior of a fracture: implications for EGS stimulation

Keywords: fractures/failles, auto-affinité, rugosité, écoulement de fluide, raideur, fermeture Numerical modeling, fractures/faults, self-affinity, roughness, fluid flow, stiffness, fracture closure, sealing, hydraulic stimulation, hydro-mechanical coupling Fracture roughness, Fluid flow, Hydraulic diffusivity, Fracture closure, Anisotropy self-propped fracture, self-affine fracture roughness, asperity yield, fracture sealing, hydraulic stimulation, hydro-mechanical coupling, pore pressure diffusion, poro-elastic stressing, fault reactivation, induced seismicity

Mots-clés Modélisation numérique,

Since I first came to France in 2018, I have met many people who have brought me great joy and helped me out of difficult situations. I still remember Mingxiang picking me up at CDG airport and accompanying me to Fontainebleau, even also as a newcomer. And Peng, who was the first person to visit Paris with me. Thank you, these two guys, who were also my roommates during French training in Shanghai. Then I met Shuaitao, Hao, Kaiwen, Dingan, Yang & Zimin, Tiantian etc. Thank you for all those great time and help, and special thanks to Shuaitao for hosting me for long time. I would also like to thank my previous colleagues at the Geoscience center in Mines Paristech -PSL university: Fédéric, Dominique, Dac, Sara, Hafsa and Houssam, for any support from you.

I am deeply grateful to Prof. François Cornet for helping me to continue my PhD studies in Strasbourg. Sadly he passed away suddenly and unexpectedly in 2020. And those who have helped me a lot during the very difficult times: Yu, Jinzhou, Zhixiong, Yexia and Fujiao, are highly appreciated.

VII

Résumé

L'exploitation de l'énergie souterraine telle que les ressources géothermiques profondes implique souvent un milieu discontinu, par exemple des failles et des fractures. Dans le système géothermique amélioré (EGS), ces structures imparfaites servent de pièces faibles au sens mécanique mais d'éléments forts au sens hydraulique. Dans le cas des premiers, ils ont généralement une résistance inférieure à celle de la roche intacte et peuvent provoquer une instabilité structurelle globale, telle que la fermeture d'une fracture ou un glissement de faille. Ils peuvent affecter de manière significative la circulation des fluides ou même induire des tremblements de terre. Dans le cas de ces derniers, profitant des surfaces rugueuses, ils sont susceptibles de former des vides de contact partiels et de canaliser par la suite des flux généralement beaucoup plus perméables que la matrice rocheuse environnante, et peuvent donc servir de conduit principal pour les fluides souterrains.

Pour éviter ces problèmes indésirables et obtenir une bonne circulation des fluides dans le réservoir, il est impératif de mieux comprendre le comportement hydromécanique de ces fractures rugueuses à des échelles variées du laboratoire au domaine d'intérêt. Dans l'EGS, la cible d'extraction de chaleur a subi un processus évolutif progressif. Initialement, la tâche a consisté à utiliser des injections massives de fluides pour créer de nouvelles fractures ou réseaux de fractures. Ensuite, on est progressivement passé à de multiples stratégies de stimulation (par exemple, cisaillement hydraulique, stimulation chimique et thermique) des fractures préexistantes pour améliorer leur perméabilité. À l'heure actuelle, la méthode la plus populaire consiste à cibler directement une grande faille avec une perméabilité naturelle élevée et un échange de chaleur sur une grande surface. Fait intéressant, de nombreuses mesures de surfaces rocheuses, depuis des échantillons de carottes jusqu'à des failles régionales sur plusieurs échelles de longueur, montrent des propriétés d'échelle auto-affines. Cela nous donne la possibilité de reproduire des fractures géologiques dans un ordinateur pour effectuer des simulations hydrauliques et/ou mécaniques à un coût relativement faible.

Ainsi, notre approche est basée sur le développement d'un schéma numérique qui nous per-XV RÉSUMÉ met de simuler l'écoulement de fluide et la déformation mécanique d'une fracture ou d'une roche fracturée. Ceci est fait principalement dans le cadre de la méthode des éléments finis, à l'aide d'un simulateur open source MOOSE/GOLEM dédié à la résolution de problèmes liés aux processus T-H-M (Thermique-Hydraulique-Mécanique) couplés dans les réservoirs fracturés (Thèmes 1 et 2). En outre, nous étudions également la fracturation hydraulique dans le massif rocheux en considérant l'intersection de la fracture hydraulique et des fractures naturelles à l'aide du package ABAQUS (Thème 3).

Thème 1 : Diffusion hydraulique le long d'une fracture grossière. La diffusion de la pression est d'une grande importance dans le développement des réservoirs géothermiques souterrains et considérée comme l'un des mécanismes de la sismicité induite. L'objectif ici est de quantifier une diffusivité hydraulique effective d'une fracture grossière partiellement ouverte (échelle de ∼ 500 m) lors de la fermeture progressive en considérant à la fois les effets temporels et spatiaux, en particulier pour une grande faille unique dominée par l'écoulement du fluide. Afin de résoudre l'équation de diffusion de pression transitoire, nous utilisons une représentation 2D de l'ouverture de fracture 3D. La diffusivité hydraulique effective est obtenue en comparant les résultats numériques avec une solution analytique du modèle de plaques parallèles. L'influence de la rugosité de la fracture sur la diffusion hydraulique est soulignée et les résultats sont comparables à ceux rapportés lors de certains séismes. La méthodologie peut être étendue pour quantifier la diffusivité hydraulique effective des réseaux de fractures en incluant la rugosité de chaque fracture.

Thème 2 : Impact de la cimentation sur la rigidité et la perméabilité d'une fracture de grande échelle. L'objectif est de comprendre comment la précipitation minérale peut influencer la perméabilité de la fracture, comment on peut prédire le degré de scellement et quelles implications cela aura sur la stimulation chimique. Pour ce faire, nous développons un modèle 3D hydro-mécanique couplé où une grande fracture rugueuse unique est enchâssée dans un réservoir granitique imperméable (256 m × 256 m × 256 m). Nous imitons le scellement minéral en ajoutant progressivement une couche d'épaisseur constante aux parois de la fracture avant un processus de fermeture mécanique. En simplifiant la rhéologie de contact, nous supposons ici un régime rigide-plastique pour la fermeture mécanique. A chaque étape de fermeture de la fracture (mécanique ou scellement), nous calculons la rigidité de la fracture en appliquant une petite perturbation de contrainte sur le sommet du modèle et calculons la perméabilité en imposant une petite chute de pression (10 Pa/256 m) dans les deux x-et y-directions jusqu'à ce que les seuils de percolation du débit soient atteints. Pour relier les deux paramètres, nous dérivons une loi basée sur des études antérieures avec laquelle les résultats numériques sont cohérents. Cette loi peut être utilisée pour évaluer l'état hydraulique/mécanique du réservoir fracturé. De plus, l'anisotropie de perméabilité observée peut XVI RÉSUMÉ guider l'utilisation de la stimulation chimique pour améliorer la transmissivité du flux. Afin de considérer l'effet de la température ou de la dissolution des minéraux, un modèle T-H-M-C (Thermique-Hydraulique-Mécanique-Chimique) entièrement couplé reste à développer. Thème 3 : Stimulation hydraulique autour du puits d'injection en tenant compte de la fracturation hydraulique et du cisaillement hydraulique. L'objectif est d'étudier les réponses hydromécaniques telles que la propagation de la fracture, la diffusion de la pression, la contrainte poroélastique et le glissement de la fracture préexistante lors de l'injection de fluide et après la fermeture. Le chemin de stimulation est préréglé comprenant deux segments de fracture hydraulique et deux fractures préexistantes. Les fractures sont modélisées par des éléments cohésifs couplés pressiondéplacement spéciaux, qui utilisent un concept de zone cohésive pour décrire la zone partiellement endommagée près de la pointe de la fissure et entre la zone fissurée et le matériau intact. La propagation de la rupture se produit lorsque les éléments cohésifs sont complètement endommagés en fonction de l'énergie critique de rupture. Le modèle résout plusieurs processus couplés : la déformation poro-élastique de la roche, l'écoulement de fluide interstitiel, l'écoulement de fracture et la propagation de fracture. Le glissement de cisaillement sismique de la fracture préexistante suit la loi de frottement de Coulomb et est utilisé pour calculer la magnitude du moment. Des études paramétriques sur la perméabilité de la matrice, le coefficient de frottement, l'orientation de la fracture et le schéma d'injection montrent que la sismicité induite est liée à la fois à l'état initial et à l'intervention humaine. La diffusion de la pression du fluide ainsi que le changement de contrainte poro-élastique sont responsables de la sismicité après fermeture. Le modèle reste à étendre pour simuler l'injection de fluide dans des réseaux de fractures complexes.

Nos études sont dédiées à l'application des techniques numériques dans les fractures et les réservoirs fracturés, en particulier l'EGS. Les simulations offrent de nouvelles perspectives pour comprendre les comportements hydrauliques/mécaniques d'une fracture rugueuse, la relation entre l'écoulement de fluide et la déformation mécanique d'une fracture rugueuse, et les réponses couplées H-M des roches fracturées. L'analyse ouvre de nouvelles perspectives pour relier le modèle à des sujets d'intérêt dans le développement de l'EGS, tels que la diffusion de la pression, l'amélioration de la perméabilité, la fermeture des fractures, la réactivation des failles et la sismicité induite. [START_REF] Sausse | 3d model of fracture zones at soultz-sous-foręts based on geological data, image logs, induced microseismicity and vertical seismic profiles[END_REF]; The blue and red arrows indicate injection and production well(s), respectively 1.4 3D view of the induced seismicity in Soultz-Sous-Forêts, modified from [START_REF] Sausse | 3d model of fracture zones at soultz-sous-foręts based on geological data, image logs, induced microseismicity and vertical seismic profiles[END_REF]; The circles indicate seismic events . . . . . . . . . . . . . . . . . . .

XVII

List of Figures

1.5 Fluid injection-induced seismicity related to pore pressure diffusion from GPK1, Soultz-Sous-Forêts, 1993, modified from [START_REF] Shapiro | Large-scale in situ permeability tensor of rocks from induced microseismicity[END_REF] . . . . . . . . . . . 1.6 Vertical tensile fracture in the GPK1 borehole at Soultz-Sous-Forêts which extends over a depth range of about 3 m, modified from [START_REF] Brudy | Drilling-induced tensile wall-fractures: implications for determination of in-situ stress orientation and magnitude[END_REF] . . . . . dicular (H ⊥ = 0.83 ± 0.03) to the slip direction; derived from [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. 31 2.3 Log-log graph of Fourier power spectra C(q) of linear profiles as a function of the wavenumber q for our synthetic self-affine surface (large blue dots) and for a group of natural fault surfaces (small dots) from [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. The synthetic selfaffine surface is generated with H = 0.8 and h rms = 0.1 m at the 512 m scale to match field observations . . XXII 3.9 a) Fitting of the fracture normal stiffness as a function of fracture closure during the mechanical stage for F1 using the stress-linear stiffness approach. b) Linear regression of the stiffness-stress data of the stress-linear stiffness approach. . . . . .

3.10

Comparison of different approaches to obtain the normal stress against fracture closure during the mechanical closure phase for fault F1: the red curve indicates the stress-linear method in the Bandis model; the blue curve indicates the elastic-limit method; and the green curve indicates the elasto-plastic (EP) model using the solution from the boundary element method as in [START_REF] Röttger | contact. engineering-create, analyze and publish digital surface twins from topography measurements across many scales[END_REF]. . . . . . . . . .

3.11

Log-linear plot of the fracture normal stiffness against the fault volume for all studied faults. The pink curve denotes the prediction from Eq. (3.16): )) for the 6 studied faults, which is fitted by k

κ ≈ κ 0 exp (-V f /V c f ).
≈ (B/A 2 ) • (V f -V perc f ) 2
, where A is the total block area and B = 1.24 × 10 -4 is a geometrical prefactor. . . . . . . . . . . . . .

(a)

Power spectral density C(q) of the fault topography F1 along both directions (X and Y). The dashed lines indicate the domain of measured surfaces from the field and laboratories [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF] (top line is C(q) = 2 × 10 -2 q -3 , medium line is C(q) = 2 × 10 -4 q -3 and bottom line is C(q) = 2 × 10 -6 q -3 ); (b) RMS h rms evolution with scale L along X and Y for the same synthetic fault. . . . . . . . . . 2022)) or for electricity generation at higher temperatures (> 150 °C, e.g., [START_REF] Muffler | Methods for regional assessment of geothermal resources[END_REF]; [START_REF] Genter | Overview of the current activities of the european egs soultz project: from exploration to electricity production[END_REF]; [START_REF] Bahadori | A review of geothermal energy resources in australia: current status and prospects[END_REF]; [START_REF] Hirschberg | Energy from the Earth: Deep Geothermal as a Resource for the Future?[END_REF]; [START_REF] Paulillo | Geothermal energy in the uk: The life-cycle environmental impacts of electricity production from the united downs deep geothermal power project[END_REF]). A prerequisite for the success of such systems is maintaining good heat fluid circulation. Deep geothermal energy recovery typically employs a doublet well system involving an injection well and a production well (Fig. 1.1). EGS sites can be categorized into hydrothermal systems and petrothermal systems. The former are suitable for existing aquifers with adequate fractures, water flow, and sufficient temperatures, and therefore require generally minimal stimulation. However, such condition is not always met, in particular at a great depth where the host rock is typically impermeable. Petrothermal systems are thus designed to improve reservoir transmissivity and productivity, which is formerly referred to as Hot Dry Rock (HDR) that originated from Los Alamos National Labs [START_REF] Mortensen | Hot dry rock: a new geothermal energy source[END_REF]. The involved hydraulic fracturing technique often requires high pressures to overcome the sum of the minimum principal stress and the rock tensile strength and create new fractures. It is therefore environmentally and economically unfriendly but still has gathered plenty of important results, e.g., the Fenton Hill HDR project [START_REF] Aki | Interpretation of seismic data from hydraulic fracturing experiments at the fenton hill, new mexico, hot dry rock geothermal site[END_REF][START_REF] Norbeck | Field observations at the fenton hill enhanced geothermal system test site support mixed-mechanism stimulation[END_REF].

In between these two systems, Enhanced Geothermal Systems (EGS) are proposed to increase the permeability of the host rock with seismicity control (Fig. 1.1). With the increasing exploration of EGS, the concept of hydraulic fracturing has been gradually shifting to hydraulic shearing, aiming to improve the reservoir permeability by reactivation of pre-existing fractures, with the resulting overpressure capable of bringing about fracture shear sliding and irreversible dilation [START_REF] Mcclure | An investigation of stimulation mechanisms in enhanced geothermal systems[END_REF][START_REF] Blöcher | Hydraulic history and current state of the deep geothermal reservoir groß schönebeck[END_REF][START_REF] Bijay | Geothermal reservoir stimulation through hydro-shearing: An experimental study under conditions close to enhanced geothermal systems[END_REF]. Such shearing mechanism does not necessarily require high injection pressures, especially for those natural fractures critically stressed (i.e. sets of fractures are favorably oriented to fail in shear under the present-day stress field) or optimally oriented, but still has the ability to improve the permeability by 2-3 orders of magnitude [START_REF] Häring | Characterisation of the basel 1 enhanced geothermal system[END_REF][START_REF] Mcclure | An investigation of stimulation mechanisms in enhanced geothermal systems[END_REF]. In fact, previous studies also showed that these two 'end-member' mechanisms (fracture opening and shearing mechanisms) both contribute to the hydraulic injection resulting in a mixed stimulation mechanism [START_REF] Mcclure | An investigation of stimulation mechanisms in enhanced geothermal systems[END_REF]. Specifically, in the area close to the wellbore, the hydraulic fracturing mechanism might be of more importance to link the wellbore with the natural fracture network in the target formation. While in areas a few tens of meters or more away from the wellbore, hydraulic shearing would dominate the permeability enhancement process [START_REF] Cornet | The hydromechanical behaviour of a fracture: an in situ experimental case study[END_REF][START_REF] Cacace | Modelling of fractured carbonate reservoirs: outline of a novel technique via a case study from the molasse basin, southern bavaria, germany[END_REF][START_REF] Amann | The seismo-hydromechanical behavior during deep geothermal reservoir stimulations: open questions tackled in a decameter-scale in situ stimulation experiment[END_REF]. In addition to hydraulic stimulation, other techniques, such as thermal and chemical stimulation have been proposed. Thermal stimulation, using cold water injection, is mainly based on tem-perature difference (or cooling effect) that enables to induce stress changes, leading to new cracks generation or natural fracture network stimulation, thereby increasing the initial permeability. Its effect is similar to that of the hydraulic stimulation, see, e.g., the Rittershoffen project [START_REF] Baujard | Hydrothermal characterization of wells grt-1 and grt-2 in rittershoffen, france: Implications on the understanding of natural flow systems in the rhine graben[END_REF], the Raft River geothermal site [START_REF] Bradford | Thermal stimulation and injectivity testing at raft river, id egs site[END_REF]. Chemical stimulation is applicable to reservoirs rich in secondary minerals that seal rock matrix or fractures, e.g., in Soultz-sous-Forêts [START_REF] Nami | Chemical stimulation operations for reservoir development of the deep crystalline hdr/egs system at soultz-sous-forêts (france)[END_REF][START_REF] Portier | Chemical stimulation techniques for geothermal wells: experiments on the three-well egs system at soultz-sous-forêts, france[END_REF]. It consists in injecting fluids with chemical additives into the reservoir in order to dissolve/remove these secondary minerals and thereby increase fluid pathways and permeability [START_REF] Griffiths | Geothermal implications for fracture-filling hydrothermal precipitation[END_REF][START_REF] Watanabe | Novel chemical stimulation for geothermal reservoirs by chelating agent driven selective mineral dissolution in fractured rocks[END_REF]. Meanwhile, the target of EGS has been also evolving from a fracture network to a large regional fault where a higher flow rate and larger heat exchange area are expected (e.g., [START_REF] Baujard | Hydrothermal characterization of wells grt-1 and grt-2 in rittershoffen, france: Implications on the understanding of natural flow systems in the rhine graben[END_REF]; [START_REF] Schmittbuhl | Induced and triggered seismicity below the city of strasbourg, france from november 2019 to january 2021[END_REF]). This is because EGS not only needs to improve the reservoir permeability but also needs to avoid heat loss and maintain its longevity/sustainability, and this requires plenty of flow paths designed to provide an abundance of rock surface area for subsurface fluid flow [START_REF] Kluge | Sustainability of engineered fractured systems: An experimental study on hydromechanical properties[END_REF].

In summary, the objectives of EGS can be achieved by various methods depending on the geological situation, such as rock types, rock structures, tectonic conditions, in-situ stress fields, etc. However, no matter which approach is used and what challenges are faced, the fundamental and important elements of EGS cannot be skipped, which are, fractures and faults.

Fractures and faults

In rock mechanics, fractures are open spaces between opposite rock surfaces due to fracturing. They can be classified into faults and joints [START_REF] Jaeger | Effect of fluid pressure heterogeneity on injection-induced fracture activation[END_REF]. The major difference between the two is lying on whether the relative displacement parallel to the nominal fracture plane has occurred or not. In structural geology, faults are unique and complex objects consisting of a fine-grained fault core and a surrounding damage zone [START_REF] Caine | Fault zone architecture and permeability structure[END_REF][START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF][START_REF] Cacace | Modelling of fractured carbonate reservoirs: outline of a novel technique via a case study from the molasse basin, southern bavaria, germany[END_REF]. For the purpose of this study, we do not differentiate joints from fractures but use the global term 'fracture' instead. In addition, we simplify the fault structure by neglecting the damage zone. The fault damage zone can affect fluid flow properties but depends on the fracture size, orientation, spatial distribution, etc [START_REF] Odling | Permeability scaling properties of fault damage zones in siliclastic rocks[END_REF][START_REF] Nguyen | Integrating damage zone heterogeneities based on stochastic realizations of fracture networks for fault stability analysis[END_REF]. In this thesis, we use synthetically generated fault geometry and focus on the influence of the fault surfaces' roughness. Therefore the impact of the damage zone is beyond the scope of this thesis.

As mechanically weak parts of the rock mass, fractures and faults control the stability of subsurface structures but simultaneously serve as major conduits for fluids due to their common higher permeability, although they can also be barriers in less frequent cases [START_REF] Council | Rock fractures and fluid flow: contemporary understanding and applications[END_REF][START_REF] Fisher | Fluid-flow properties of faults in sandstone: The importance of temperature history[END_REF][START_REF] Bense | Faults as conduit-barrier systems to fluid flow in siliciclastic sedimentary aquifers[END_REF][START_REF] Agosta | From fractures to flow: A field-based quantitative analysis of an outcropping carbonate reservoir[END_REF]. Because of their importance, the related study has been a hot topic for the past several decades [START_REF] Zoback | Laboratory hydraulic fracturing experiments in intact and pre-fractured rock[END_REF][START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF][START_REF] Tsang | Flow channeling in a single fracture as a two-dimensional strongly heterogeneous permeable medium[END_REF][START_REF] Barton | Fluid flow along potentially active faults in crystalline rock[END_REF][START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF][START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF][START_REF] Watanabe | Diversity of channeling flow in heterogeneous aperture distribution inferred from integrated experimental-numerical analysis on flow through shear fracture in granite[END_REF][START_REF] Detournay | Mechanics of hydraulic fractures[END_REF][START_REF] Vogler | Experiments and simulations of fully hydro-mechanically coupled response of rough fractures exposed to highpressure fluid injection[END_REF]Kluge et al., 2021a). Although the understanding of faults and fractures has improved considerably, there are many aspects that still remain unknown or unclear. This is primarily because fractures are deeply buried below the Earth's surface. Quantifying precisely fracture aperture or fracture roughness at depth and varied length scales thereby becomes difficult under current techniques. In practice, obtaining fracture morphology at depth is typically based upon image logging interpretation, e.g., seismic, acoustic, and electric logs [START_REF] Luthi | Fracture apertures from electrical borehole scans[END_REF][START_REF] Leary | Deep borehole log evidence for fractal distribution of fractures in crystalline rock[END_REF][START_REF] Prioul | Fracture characterization at multiple scales using borehole images, sonic logs, and walkaround vertical seismic profile[END_REF][START_REF] Ponziani | Experimental validation of fracture aperture determination from borehole electric microresistivity measurements[END_REF][START_REF] Lai | Fracture characterization and detection in the deep cambrian dolostones in the tarim basin, china: Insights from borehole image and sonic logs[END_REF], or samples taken from boreholes [START_REF] Wang | Stereological interpretation of rock fracture traces on borehole walls and other cylindrical surfaces[END_REF][START_REF] Rogers | Volumetric fracture intensity measurement for improved rock mass characterisation and fragmentation assessment in block caving operations[END_REF]. However, these methods can only provide limited information about fractures/fracture networks, for example, to have the strike and dip of a large-scale fault, or reveal fracture traces on wellbore walls. Obtaining details of the geometric information of the fracture, such as the thickness variation, effective aperture variability, contact zone, etc, is much more difficult. Moreover, despite the fact that numerous experimental results on rock samples have been obtained, using the laboratory scale to represent the field scale remains a debate. In this case, the use of analogues from outcrops is an alternative to reproduce fracture surfaces at varied scales.

By using different scanner devices (e.g., profilometer, LiDAR, white light interferometer) and statistical analysis approaches (e.g., Fourier Power Spectrum Method, Standard Deviation of the Correlation Function Method), many researchers have conducted surface measurements on natural joints, artificial fractures, exhumed faults profiles, and earthquake rupture surfaces, both in the laboratory and in the field, and created a huge amount of data on the rock fracture surfaces. These studies have revealed a universal geometric property of fracture surfaces over several decades of lengths from µm to km scale, i.e., self-affine property [START_REF] Brown | Broad bandwidth study of the topography of natural rock surfaces[END_REF][START_REF] Power | Roughness of natural fault surfaces[END_REF]Schmittbuhl et al., 1993a;[START_REF] Power | Topography of natural and artificial fractures in granitic rocks: Implications for studies of rock friction and fluid migration[END_REF][START_REF] Renard | High resolution 3d laser scanner measurements of a strike-slip fault quantify its morphological anisotropy at all scales[END_REF][START_REF] Bistacchi | Fault roughness at seismogenic depths from lidar and photogrammetric analysis[END_REF][START_REF] Brodsky | Faults smooth gradually as a function of slip[END_REF][START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. In general, a statistical invariance in terms of roughness amplitude and roughness exponent exists. When we look at such surfaces in Fourier space, the power spectrum vs spatial frequency (or wavenumber) shows a linear trend in a log-log plot. The slope equal to -2H-1 or -2H-2 (for 1D profiles and 2D surface respectively), with H the roughness exponent [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF][START_REF] Brodsky | Constraints from fault roughness on the scaledependent strength of rocks[END_REF]. The roughness exponent, also called Hurst exponent, was found to be about 0.8 for artificial tensile fractures or faults perpendicular to the sliding direction, and about 0.6 for faults along the sliding direction [START_REF] Nigon | Multiscale characterization of joint surface roughness[END_REF]Schmittbuhl et al., 1995a;[START_REF] Candela | Characterization of fault roughness at various scales: Implications of three-dimensional high resolution topography measurements[END_REF]. In addition, no roll-off wavelength was found over multiple length scales [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]Schmittbuhl et al., 1993a). The self-affinity of fracture surfaces is crucial since it appears to be one of the best models to account for the multi-scale property of the fracture geometry. Having the fracture geometrical model provides us the possibility to generate fractures synthetically on different scales.

Measurements of fracture aperture distribution are also associated with measurements of fracture surfaces. The aperture is defined as the forming space between the opposite fracture surfaces. By measuring the two surfaces' topography separately, the fracture aperture can then be constructed, see, for example, [START_REF] Watanabe | Determination of aperture structure and fluid flow in a rock fracture by high-resolution numerical modeling on the basis of a flow-through experiment under confining pressure[END_REF]; Neuville et al. (2012a);[START_REF] Hofmann | Transmissivity of aligned and displaced tensile fractures in granitic rocks during cyclic loading[END_REF]; Kluge et al. (2021b).

However, one issue raised by measuring crack surfaces is how to ensure the accuracy of matching points on two surfaces facing each other. To avoid this issue, several approaches have been developed to directly measure the fracture aperture distribution, including injection resin into a fracture [START_REF] Hakami | Aperture measurements and flow experiments on a single natural fracture[END_REF], casting techniques, i.e. making a replica of the fracture's void space [START_REF] Yeo | Effect of shear displacement on the aperture and permeability of a rock fracture[END_REF][START_REF] Renshaw | Measuring fracture apertures: A comparison of methods[END_REF], non-invasive techniques, such as sonic method [START_REF] Hornby | Fracture evaluation using reflected stoneleywave arrivals[END_REF], electrical method [START_REF] Luthi | Fracture apertures from electrical borehole scans[END_REF], optical method [START_REF] Detwiler | Measurement of fracture aperture fields using transmitted light: An evaluation of measurement errors and their influence on simulations of flow and transport through a single fracture[END_REF][START_REF] Ameli | High-resolution fracture aperture mapping using optical profilometry[END_REF], nuclear magnetic resonance imaging [START_REF] Renshaw | Measuring fracture apertures: A comparison of methods[END_REF], X-ray computed tomography (Pyrak- [START_REF] Pyrak-Nolte | Volumetric imaging of aperture distributions in connected fracture networks[END_REF][START_REF] Keller | High resolution, non-destructive measurement and characterization of fracture apertures[END_REF][START_REF] Ramandi | Micro-ct image calibration to improve fracture aperture measurement[END_REF]. According to the measuring results, the fracture aperture field can be described by varied distributions. The most common two distributions are normal distribution [START_REF] Hakami | Aperture measurements and flow experiments on a single natural fracture[END_REF][START_REF] Pyrak-Nolte | Volumetric imaging of aperture distributions in connected fracture networks[END_REF][START_REF] Yeo | Effect of shear displacement on the aperture and permeability of a rock fracture[END_REF]Neuville et al., 2012a) and lognormal distribution [START_REF] Snow | The frequency and apertures of fractures in rock[END_REF][START_REF] Pyrak-Nolte | Volumetric imaging of aperture distributions in connected fracture networks[END_REF][START_REF] Keller | High resolution, non-destructive measurement and characterization of fracture apertures[END_REF][START_REF] Xu | Aperture measurements and seepage properties of typical single natural fractures[END_REF], while a few follows gamma distribution [START_REF] Tsang | Channel model of flow through fractured media[END_REF][START_REF] Xu | Aperture measurements and seepage properties of typical single natural fractures[END_REF] and other distributions [START_REF] Ramandi | Micro-ct image calibration to improve fracture aperture measurement[END_REF][START_REF] Xu | Aperture measurements and seepage properties of typical single natural fractures[END_REF].

However, fracture void geometry cannot only be described by the aperture distributions but is also influenced by the spatial correlations. Fractures with the same aperture frequency distributions may have different spatial correlations between the apertures on the fracture surface [START_REF] Hakami | Aperture measurements and flow experiments on a single natural fracture[END_REF]. The spatial correlations can be quantified by a semi-variogram, which shows sill (the plateau level) and range (the lag distance at which the sill is reached). The range is called the correlation length. The correlation length of the aperture was found typically much less than the fracture length at sample scale (Brown et al., 1986;[START_REF] Hakami | Aperture measurements and flow experiments on a single natural fracture[END_REF][START_REF] Pyrak-Nolte | Volumetric imaging of aperture distributions in connected fracture networks[END_REF][START_REF] Keller | High resolution, non-destructive measurement and characterization of fracture apertures[END_REF], although the correlation length of the fracture aperture in the field is difficult to determine. Therefore, unlike fracture surfaces, the fracture aperture is self-affine only at scales smaller than the correlation length, beyond which the Fourier spectrum of the aperture is constant, namely without self-affinity [START_REF] Brown | Broad bandwidth study of the topography of natural rock surfaces[END_REF][START_REF] Méheust | Scale effects related to flow in rough fractures[END_REF].

The spatial correlations of fracture aperture are affected by a variety of factors such as shear displacement and normal stress, e.g., [START_REF] Plouraboué | Aperture of rough cracks[END_REF][START_REF] Yeo | Effect of shear displacement on the aperture and permeability of a rock fracture[END_REF][START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF] and have a significant impact on the hydraulic and mechanical properties of a fracture, e.g., [START_REF] Yeo | Effect of shear displacement on the aperture and permeability of a rock fracture[END_REF]; Pyrak- [START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF]; [START_REF] Méheust | Scale effects related to flow in rough fractures[END_REF]; [START_REF] Petrovitch | Scaling of fluid flow versus fracture stiffness[END_REF]; Wang et al. (2016); [START_REF] Li | On the relationship between normal stiffness and permeability of rock fractures[END_REF]; [START_REF] Lenci | A lubrication-based solver for shearthinning flow in rough fractures[END_REF]. Among these studies, Pyrak-Nolte and [START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF] found that fractures with long-range correlations have only a few flow channels and the contact area is limited in a few regions. The correlation length decreases with increasing normal stress, and the flow decreases rapidly with changes in the fracture's normal stiffness. While fractures with short-range correlations are more likely to have multiple flow paths and uniformly distributed contact areas, and the flow changes slowly with changes in the fracture's normal stiffness. In addition, [START_REF] Méheust | Scale effects related to flow in rough fractures[END_REF] systematically studied the influence of the mismatch scale (i.e. the correlation length) on the hydraulic behavior of a fracture. They showed that at the fixed mechanical aperture, the deviation from the cubic law decreases as the correlation length decreases, and the loss of the correlations at the fracture scale has a strong impact on the flow behavior. Therefore, the correlations at the fracture scale are of great importance such as understanding the flow channels caused by large apertures, see also [START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF], which will also be the main focus of this thesis. The effect of correlations smaller than the fracture scale will be studied in the future.

Over the past half-century, a great deal of research has been dedicated to understanding the hydraulic and mechanical behavior of fractures and faults, mainly including field tests, laboratory experiments, numerical simulations, as well as some theoretical analysis.

In many studies, due to the speciality and requirement of the problem, fractures are typically simplified to the parallel plate model, and accordingly fluid flow through them is governed by the cubic law [START_REF] Lei | Modelling fluid injection-induced fracture activation, damage growth, seismicity occurrence and connectivity change in naturally fractured rocks[END_REF][START_REF] Carrier | Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model[END_REF][START_REF] Detournay | Propagation regimes of fluid-driven fractures in impermeable rocks[END_REF][START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF], which describes that the flow rate is proportional to the cube of the fracture (mean) aperture. However, early experiments showed that this simplification cannot stand when the fracture is subjected to relatively high normal stress (e.g., 20 MPa [START_REF] Witherspoon | Validity of cubic law for fluid flow in a deformable rock fracture[END_REF]), this is due to the fact that the fracture surfaces are more or less rough, as discussed above. Under normal stress, these rough surfaces can thus create flow channels and contact areas and therefore deviate from the parallel plate model, as also investigated by numerous authors hereafter, e.g., [START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF]; [START_REF] Yeo | Effect of shear displacement on the aperture and permeability of a rock fracture[END_REF]; [START_REF] Brown | Experimental observation of fluid flow channels in a single fracture[END_REF]; [START_REF] Méheust | Flow enhancement of a rough fracture[END_REF]. In addition, a lot of theoretical models have been developed for the fluid flow through a single rough fracture. For example, [START_REF] Tsang | Hydromechanical behavior of a deformable rock fracture subject to normal stress[END_REF] replaced the aperture with the weighted average aperture in the cubic law. They also proposed a channel conceptual model to study the flow through fractured rock [START_REF] Tsang | Channel model of flow through fractured media[END_REF]. Afterward, Robert Zimmerman and his co-authors mathematically studied the flow problem in a single fracture. These include, but are not limited to: obtaining analytical expressions for the fracture permeability of a fracture with cylindrical asperities [START_REF] Kumar | Permeability of a fracture with cylindrical asperities[END_REF]; examining the validity of the lubrication approximation [START_REF] Zimmerman | Lubrication theory analysis of the permeability of rough-walled fractures[END_REF] and revising the hydraulic aperture by considering the contact area and the mean aperture [START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF]; discussing conditions for simplification of Navier-Stokes equations to the Reynolds equations (i.e. Local cubic law, LCL) [START_REF] Zimmerman | Fluid flow in rock fractures: From the navier-stokes equations to the cubic law[END_REF]. The Reynolds equations are valid typically at low Reynolds number (Re < 1) and with small variations of apertures, see also [START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF].

Field tests are the most intuitive and direct way to reveal the in situ fracture properties. It has the most explicit engineering purpose and can provide first-hand knowledge to guide the production of subsurface fluids. For example, tracer tests are commonly used in the field to characterize the fluid flow properties (e.g., permeability) [START_REF] Saar | Geothermal heat as a tracer of large-scale groundwater flow and as a means to determine permeability fields[END_REF]. In addition, [START_REF] Rutqvist | A cyclic hydraulic jacking test to determine the in situ stress normal to a fracture[END_REF] used a cyclic hydraulic test to determine the in situ stress normal to a fracture. Moreover, through testing a single natural granitic fracture at ∼ 1 m scale, [START_REF] Cornet | The hydromechanical behaviour of a fracture: an in situ experimental case study[END_REF] found that hydraulic aperture and mechanical aperture are comparable only above a fracture opening threshold and channeling may control the fluid flow; and (Guglielmi et al., 2015b) observed a significant permeability enhancement of a 20-m critically stressed fault in shale, which was only related to a small shear slip. However, for these field tests, it is often hard to control their experimental conditions due to lots of uncertainties in the subsurface thereby they are typical of high-cost risk if the results are undesired. Moreover, because of the complexity of the reservoir, most of the results are highly site-dependent, which means that relying on them for other projects may be difficult or unavailable.

Laboratory experiments can capture some phenomena that are difficult to observe in the field. As an example, with the help of microscope and CT scanner, micromorphological changes of the fracture surface prior to and after the hydraulic experiments can be determined [START_REF] Kluge | Sustainability of engineered fractured systems: An experimental study on hydromechanical properties[END_REF]Ji et al., 2020;[START_REF] Cheng | Hydromechanical investigations on the self-propping potential of fractures in tight sandstones[END_REF]. Because the conditions are relatively readily controllable in the laboratory, it is generally accepted that measurements from the laboratory are able to facilitate the understanding of field observations (e.g., [START_REF] Rummel | Application of laboratory fracture mechanics data to hydraulic fracturing field tests[END_REF]; [START_REF] Zimmermann | Hydraulic stimulation of a deep sandstone reservoir to develop an enhanced geothermal system: Laboratory and field experiments[END_REF]) and provide parameters for numerical models (e.g., [START_REF] Vogler | Experiments and simulations of fully hydro-mechanically coupled response of rough fractures exposed to highpressure fluid injection[END_REF][START_REF] Vallier | Large scale hydro-thermal circulation in the deep geothermal reservoir of soultz-sous-forêts (france)[END_REF]). Among these experiments, several common purposes for a single fracture are dedicated to discovering its fluid flow behavior: normal closure behavior, shear dilation behavior, and hydro-mechanical coupling behavior. To name just a few, fluid flow through a single fracture has been shown to deviate from the cubic law due to roughness [START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF][START_REF] Méheust | Flow enhancement of a rough fracture[END_REF][START_REF] Kluge | Permeability of displaced fractures[END_REF]; [START_REF] Goodman | A model for the mechanics of jointed rock[END_REF] and [START_REF] Bandis | Fundamentals of rock joint deformation[END_REF] proposed empirical laws for describing the fracture closure under normal stress based on experimental results. Although laboratory experiments have yielded significant results, upscaling them to the engineering scale remains controversial [START_REF] Viswanathan | From fluid flow to coupled processes in fractured rock: Recent advances and new frontiers[END_REF].

As a complementary technique, numerical models play an essential role in understanding fracture behavior. On the one hand, they can serve as validation for field and laboratory measurements. For example, numerical studies on the flow through a rough fracture revealed the flow channeling behavior [START_REF] Blöcher | Permeability of matrix-fracture systems under mechanical loading-constraints from laboratory experiments and 3-d numerical modelling[END_REF][START_REF] Brown | Applicability of the reynolds equation for modeling fluid flow between rough surfaces[END_REF][START_REF] Skjetne | High-velocity flow in a rough fracture[END_REF][START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF][START_REF] Schmittbuhl | Fracture morphology and viscous transport[END_REF], which was also found in the laboratory; [START_REF] Kling | Numerical simulations and validation of contact mechanics in a granodiorite fracture[END_REF]; [START_REF] Li | On the relationship between normal stiffness and permeability of rock fractures[END_REF] showed that the granodiorite fracture normal closure behavior can be well described by an elastoplastic contact model. On the other hand, numerical studies are capable of providing new insights to overcome the constraints imposed by observations. For instance, on the basis of the self-affine geometry of a fracture, numerical simulations can be performed without the limitation of the number of fractures, thus are able to obtain generalized laws [START_REF] Petrovitch | Scaling of fluid flow versus fracture stiffness[END_REF][START_REF] Wang | Development of an empirical model relating permeability and specific stiffness for rough fractures from numerical deformation experiments[END_REF][START_REF] Pyrak-Nolte | Approaching a universal scaling relationship between fracture stiffness and fluid flow[END_REF][START_REF] Morris | Numerical investigation of alternative fracture stiffness measures and their respective scaling behaviours[END_REF][START_REF] Li | On the relationship between normal stiffness and permeability of rock fractures[END_REF]. In addition, with the help modern high-performance computer, the numerical technique has undergone rapid development, e.g., numerical experiments on large-scale fractures or on 3D fractures become available [START_REF] Molinero | Large-scale modeling of reactive solute transport in fracture zones of granitic bedrocks[END_REF][START_REF] Blöcher | 3d numerical modeling of hydrothermal processes during the lifetime of a deep geothermal reservoir[END_REF]Wang et al., 2016;[START_REF] Blöcher | Evaluating micro-seismic events triggered by reservoir operations at the geothermal site of groß schönebeck (germany)[END_REF][START_REF] Chen | Heat transfer in a 3d rough rock fracture with heterogeneous apertures[END_REF]. If economic cost is also taken into account, numerical simulation is no doubt a superior alternative to field tests and laboratory experiments.

Objective

Based on the previous studies and the above analysis, a simple process of EGS development and the related techniques/issues involved can be depicted in Fig. 1.2. As can be seen, EGS involves a great deal of work to be done in terms of various technologies and research methods.

The challenges to be overcome behind EGS stem from two main areas. On the one hand, how to maintain a sustainable heat source and sufficient subsurface flow rates at a level of commercial interest during the designed lifetime? i.e., the sustainability of the entire system [START_REF] Kumari | Sustainable development of enhanced geothermal systems based on geotechnical research-a review[END_REF][START_REF] Kluge | Sustainability of engineered fractured systems: An experimental study on hydromechanical properties[END_REF]. This requires considerable effort in investigating rock and fracture properties, rock-fluid interactions, and injection strategies. Furthermore, more attention should be paid to the early stage work, e.g., an EGS site targeting a large-scale fault maybe with higher perme-Figure 1.2: Simplified flow chart on harnessing geothermal energy by enhancing geothermal reservoirs properties and the associated techniques/issues, modified from [START_REF] Kluge | Sustainability of engineered fractured systems: An experimental study on hydromechanical properties[END_REF] ability and larger heat exchange area. On the other hand, how to control the margin effects within an acceptable limit? e.g., wellbore collapse, fluid leak-off, and induced seismicity. Among these potential risks, induced seismicity is of most interest because it is likely to reach intensities felt by people and even threaten life and properties. Some EGS projects that were eventually stopped due to induced earthquakes include the Basel EGS project with a maximum magnitude of 3.4 in 20063.4 in (Häring et al., 2008)), Pohang EGS project with a maximum magnitude of 5.4 in 2017 [START_REF] Kim | Assessing whether the 2017 mw 5.4 pohang earthquake in south korea was an induced event[END_REF] and Geoven EGS project with a maximum magnitude 3.6 in 2020 [START_REF] Schmittbuhl | Induced and triggered seismicity below the city of strasbourg, france from november 2019 to january 2021[END_REF]. Prevention or prediction of induced earthquakes requires a better understanding of the associated physical and mechanical mechanisms, for example, coupled T-H-M-C processes, pressure diffusion, stress changes, and fault slip. It should also be mentioned that these two aspects have commonly interacted with each other, e.g., the appropriate injection strategies should not only improve the poor hydraulic connections between boreholes but also need to avoid unwanted induced earthquakes.

In this study, focusing tightly on the topic of fractures (either natural or newly created/induced), we use numerical simulations to explore their hydraulic and mechanical properties. The aim of the numerical study is to further the understanding of pressure diffusion, mechanical response, and the interrelated hydro-mechanical relationship of a fracture at the reservoir scale. Specifically, we consider the following three aspects. 1) To investigate the transient pressure diffusion process inside a fracture, determine the hydraulic diffusivity and analyze the influence of roughness; This corresponds to Chapter 2.

2) To investigate the relationship between the fluid flow (permeability) and mechanical responses (stiffness) of a single rough fracture embedded in a porous rock matrix and to reveal how it may impact the chemical stimulation; This corresponds to Chapter 3.

3) To investigate the hydraulic fracturing and hydraulic shearing mechanisms during hydraulic stimulation and to analyze the influence of rock properties, fracture properties, and injection schemes on the hydro-mechanical response. This corresponds to Chapter 4.

For aspects 1) and 2), it is worth noting that we focus mainly on a single fracture study. Typically, the hydro-mechanical properties of fractures can be studied either on a single fracture or fracture networks. However, we first consider a single fracture as it is the more fundamental element of fractured rocks. In addition, as can be seen from Fig. 1.3, four wellbores are connected with a large fault zone in the well-known EGS project Soultz-Sous-Forêts. The figure shows that the large fault zone controls the fluid circulation of the EGS system under the triaxial stress conditions at a depth of about 5 km and a temperature of 200 • C. Therefore, it is critical to understand how fluid flow through this single fault, and how the mechanical part affects its hydraulic behavior (such as under normal stress, see Chapter 3). Fig. 1.4 shows that induced seismicity occurs mainly along the fault zones in Soultz-Sous-Forêts. This can be related to the pore pressure diffusion as shown in Fig. 1.5. The hydraulic diffusivity D can be estimated from the distance from the injection wellbore r as a function of the time t. Therefore, it is also essential to investigate the pore pressure diffusion along the fault zone. This will be the topic of Chapter 2.

For aspect 3), as mentioned above, the hydraulic stimulation process involves hydraulic shearing and hydraulic fracturing mechanisms, in particular around the injection wellbore. They may co-Figure 1.3: 3D view of the main fault zone in Soultz-Sous-Forêts, modified from [START_REF] Sausse | 3d model of fracture zones at soultz-sous-foręts based on geological data, image logs, induced microseismicity and vertical seismic profiles[END_REF]; The blue and red arrows indicate injection and production well(s), respectively exist and interact with each other [START_REF] Rutledge | Faulting induced by forced fluid injection and fluid flow forced by faulting: An interpretation of hydraulic-fracture microseismicity, carthage cotton valley gas field, texas[END_REF][START_REF] Mcclure | An investigation of stimulation mechanisms in enhanced geothermal systems[END_REF]. The induced tensile fracture was observed on the wellbore in Soultz-Sous-Forêts (e.g. on the GPK1, Fig. 1.6), indicating the importance of the hydraulic fracturing around the injection wellbore. The interactions of the natural fracture and tensile fracture depend largely on the initial stress state, injection pressure, fracture distribution and orientation, etc. Strengthening the understanding of this coupled process is significant for optimizing the stimulation strategies. As a consequence, it is necessary and valuable to exert efforts on this topic (this will be addressed in Chapter 4).

Thesis outline

Following the objectives of advancing the understanding of hydro-mechanical properties of a single fracture and the hydraulic stimulation process in the vicinity of the injection borehole, the thesis is organized based on three articles, which are presented in the following Chapters 2, 3 and 4, respectively.

Figure 1.5: Fluid injection-induced seismicity related to pore pressure diffusion from GPK1, Soultz-Sous-Forêts, 1993, modified from [START_REF] Shapiro | Large-scale in situ permeability tensor of rocks from induced microseismicity[END_REF] Chapter 2 introduces a methodology to determine the hydraulic diffusivity of a rough fracture, which is step-wise closed from an initial state with a few contact points up to the flow percolation threshold, under the hypothesis that the contact is perfectly plastic (a penetration model), therefore the fracture storage effect is neglected and the matrix is treated as a rigid body. The fracture geometry is described by the self-affine model. To solve the transient pressure fluid flow equation, the 3D fracture aperture is represented by a 2D plane with the aperture variations treated as a parameter. The highlight here is that the temporal and spatial effects are simultaneously considered when calculating the fracture hydraulic diffusivity. This deviates from the results estimated from the effective hydraulic aperture. It also illustrates that the fracture roughness either increases or decreases the hydraulic diffusivity and results in anisotropic flow behavior.

Chapter 3 presents a model for computing the permeability and normal stiffness of a rough fracture. The generation of fracture geometry/aperture is similar to that in Chapter 2. Here the 3D fracture is embedded in the porous rock matrix. Fracture closure here consists of discrete steps, including a mechanical process and a mineral precipitation process. At each step of closure, the model solves the Navier-stokes flow in the fracture and Darcy's flow in the matrix in order to determine the permeability, and the normal stress (assuming much larger than the fluid pressure, i.e., Figure 1.6: Vertical tensile fracture in the GPK1 borehole at Soultz-Sous-Forêts which extends over a depth range of about 3 m, modified from [START_REF] Brudy | Drilling-induced tensile wall-fractures: implications for determination of in-situ stress orientation and magnitude[END_REF] drained configuration) is imposed on the top the entire block in order to access the normal stiffness. The inherent relationship between permeability and stiffness is built and its impact on the chemical stimulation is analyzed. Chapter 4 introduces a numerical model to simulate the hydraulic stimulation in the vicinity of the injection wellbore. Both the hydraulic fracturing and hydraulic shearing mechanisms are considered. For simplicity and controllability, only a single stimulation path is modeled. The influence of the rock permeability, the fracture properties, and the injection protocols are studied parametrically. It also shows that the pressure diffusion and the poroelastic stressing both contribute to the fracture slip after the injection well is shut in.

Finally, Chapter 5 provides the general conclusions of the thesis, as well as some analysis and suggestions for potential future research directions.

Some key theoretical concepts

Below some theoretical backgrounds for each chapter are given. Note that there are some repetitions in the use of some symbols in different chapters, the reader should differentiate them in each chapter.

Derivation of the linear diffusion equation

This section involves the linear diffusion equation used in Chapter 2.

Considering one-dimensional linear flow in the parallel plate model (along x-direction) with an aperture h, the cubic law relating the volume flow rate per unit height q to the cubic of the aperture h can be written as follows [START_REF] Murphy | Semianalytical solutions for fluid flow in rock joints with pressure-dependent openings[END_REF]:

q = - h 3 12µ ∂p ∂x , (1.1)
where µ is the viscosity and p is the fluid pressure.

Considering also the impermeable rock matrix, the conservation of mass yields:

∂(ρq) ∂x = - ∂(ρh) ∂t (1.2)
where ρ and t are fluid density and time, respectively.

The fluid flow equation can be written as:

∂ ∂x h 3 12µ ∂p ∂x = -ρ ∂h ∂p + h ∂ρ ∂p ∂p ∂t (1.3)
In general, the density and viscosity have very small variations, therefore the equation above can be rewritten as:

1 12µh ∂ ∂x h 3 ∂p ∂x = c ∂p ∂t (1.4)
where the total compressiblity c is the sum of the fracture compressiblity c j and the fluid compressibility c f .

Considering that the fluid modulus K f is the reciprocal of the fluid compressiblity c f , Eq.(1.4) can be recast as:

h K f ∂p ∂t + ∇ • (hq) = 0, (1.5)
Eq.(1.5) is exactly the diffusion equation Eq.(2.13) used in Chapter 2.

Assuming very small pressure changes so that fracture deformation is negligible, Eq.(1.4) can be simplified to:

∂p ∂t = D ∂ 2 p ∂x 2 (1.6)
where the hydraulic diffusivity D = h 2 /(12µc).

Eq.(1.6) is then the linear diffusion equation Eq.(2.1) used in Chapter 2.

Fracture closure and contact mechanisms

This section involves the elasto-plastic model used in Chapter 3.

Pioneering work on fracture closure under normal stress was mainly based on laboratory experiments, followed by some empirical or semi-empirical fracture closure laws, e.g., the empirical hyperbolic law [START_REF] Goodman | A model for the mechanics of jointed rock[END_REF], the semi-logarithmic law (Barton-Bandis model) [START_REF] Bandis | Fundamentals of rock joint deformation[END_REF][START_REF] Barton | Strength, deformation and conductivity coupling of rock joints[END_REF]. These laws are convenient to use particularly in fitting the experimental data and therefore inherited and developed by numerous studies (e.g., [START_REF] Evans | The effects of fracture normal compliance on the long term circulation behavior of a hot dry rock reservoir: a parameter study using the new fully-coupled code'fracture[END_REF]; [START_REF] Lee | Empirical parameters for non-linear fracture stiffness from numerical experiments of fracture closure[END_REF]; [START_REF] Wang | A constitutive model for rock interfaces and joints[END_REF]; [START_REF] Zangerl | Normal stiffness of fractures in granitic rock: A compilation of laboratory and in-situ experiments[END_REF]) despite their disadvantages, such as the limitation of specific types of rocks and the lack of physical meaning of required fitting variables. Meanwhile, benefiting from the development of numerical techniques, many physicalbased models have been proposed to describe the fracture deformation behavior under normal stress. There are broadly several models. The first is the bearing-area model (e.g., interpenetration model [START_REF] Brown | Closure of rock joints[END_REF][START_REF] Watanabe | Numerical estimation of aperture structure and flow wetted field in rock fracture[END_REF]), which assumes rigid contact surfaces and the interpenetration of asperities such that the overlap parts are merged. The second is the independent asperity model (e.g., Greenwood and Williamson model (Greenwood and Wiliamson, 1966)) in which the rough surface is approximated by a set of non-interacting asperities. The third is the re-normalization model based on Persson's theory [START_REF] Persson | Elastoplastic contact between randomly rough surfaces[END_REF] and the approach proposed by [START_REF] Joe | Effect of roughness on the adhesive tractions between contacting bodies[END_REF], the solution of the contact problem at a given scale is deduced from the larger scale. The last type is the 'brute-force' model (e.g., [START_REF] Stanley | An FFT-based method for rough surface contact[END_REF]; [START_REF] Pastewka | Contact between rough surfaces and a criterion for macroscopic adhesion[END_REF]) solving the elastic problem based on a specific topography or continuum contact mechanics. Among these contact models, many researchers have tried to apply them to the field of rock fracture contact. Tab. 1.1 lists the most common contact models used in rock engineering. In this work, we use the elasto-plastic contact model to validate the linear stiffness-stress trend of our numerical computation. The model is integrated into a web app 'Contact.engineering' (http://contact.engineering) aimed to solve rough surface contact problems. Only the degrees of freedom of the surface are considered and the bulk is treated as an elastic half-space. In our case, since we deal with the contact of two rough surfaces, the difference between the top and bottom surfaces (i.e., the initial aperture) is used. We present below the basic theory behind this model, more details can be found in [START_REF] Röttger | contact. engineering-create, analyze and publish digital surface twins from topography measurements across many scales[END_REF].

The model uses the Boundary Element Method (BEM) to solve the elastic deformation of the contact problem. This is done in reciprocal space with the Fast Fourier Transform (FFT) technique. In linear elasticity, the displacement (normalized by RMS) and the pressure (normalized by the effective Young's modulus E * ) are related as follows:

u(x, y) = A G G(x -x ′ , y -y ′ )p(x ′ , y ′ ) dx ′ dy ′ = G(x, y) ⊛ p(x, y) (1.7)
where G(x, y) is the linear elastic surface Green's function and A G is its supported region.

The convolution G ⊛ P can be converted by Fourier transform into:

ũ(q, s) = G(q, s)p(q, s) = A G u(x, y)e -iqx-isy dxdy (1.8)
For periodic systems, the linear-elastic Green's function in reciprocal Fourier space is:

G(⃗ q) = 2/(E * |⃗ q|) (1.9)
where ⃗ q = (q, s) is the wavevector. E * is the effective young's modulus of the contact bodies:

E * = 1 (1 -ν 2 1 )/E 1 + (1 -ν 2 2 )/E 2 (1.10)
where E 1 , E 2 are Young's modulus of the two contact solids and ν 1 , ν 2 are their Poisson ratios.

For nonperiodic systems, which is also our case, the Green's function is:

πE * G(x, y) = (x + a) ln (y + b) + (y + b) 2 + (x + a) 2 (y -b) + (y -b) 2 + (x + a) 2 +(y + b) ln (x + a) + (y + b) 2 + (x + a) 2 (y -b) + (y + b) 2 + (x -a) 2 +(x -a) ln (y -b) + (y -b) 2 + (x -a) 2 (y + b) + (y + b) 2 + (x -a) 2 +(y -b) ln (x -a) + (y -b) 2 + (x -a) 2 (x + a) + (y -b) 2 + (x + a) 2 (1.11)
where a = ∆x/2 and b = ∆y/2. Typically, the surface topography data is provided on a discrete grid with a measured size L x ×L y and pixels N x ×N y , i.e., the measured resolution is ∆x = L x /N x and ∆y = L y /N y .

In order to describe the plastic part, a parameter called hardness H a [START_REF] Brown | Closure of rock joints[END_REF]) is introduced. It defines the maximum stress allowed in the contact.

The following three sections involve basic fracture mechanisms, cohesive zone method, and fault reactivation used in Chapter 4.

Basic concepts in fracture mechanics

When materials undergo fracture, the rate at which energy is transformed is called the energy release rate G, defined as the total potential energy loss per crack surface area. It is directly related to the stress intensity factor K in different mode as:

G I = K 2 I E ′ , G II = K 2 II E ′ , G III = K 2 III 2µ , (1.12)
where subscripts I, II, III indicate mode-I (tension), mode-II (in-plane shear) and mode-III (out-ofplane shear) loading, respectively. E is Young modulus with E ′ = E for plane stress and E ′ = E/(1 -ν 2 ) for plane strain. µ is the shear modulus.

Under mixed-mode loading conditions, the total energy release rate G can be superimposed as:

G = K 2 I E ′ + K 2 II E ′ + K 2 III 2µ (1.13)
Accordingly, the crack growth criterion (fracture criterion) can be denoted as:

G = G C (1.14)
where G C is the critical energy release rate (or the equivalent fracture energy) depending on the material itself. It corresponds to the fracture toughness (i.e. the critical stress intensity factor) as the form of Eq. (1.12) (e.g., K 2 IC = G IC E ′ is the Mode-I fracture toughness under Mode-I loading).

There are several commonly used analytical forms to estimate the equivalent fracture energy G C depending on different contributions from the three loading modes, as described below.

(1) The power law [START_REF] Goyal | Irreversible constitutive law for modeling the delamination process using interfacial surface discontinuities[END_REF]:

G I G IC α + G II G IIC β + G III G IIIC γ = 1 (1.15)
where α, β, γ are constants.

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

G C = G IC + (G IIC -G IC ) G II + G III G I + G II + G III α (1.16)
(3) The Reeder's law [START_REF] Reeder | 3-d mixed mode delamination fracture criteria-an experimentalist's perspective[END_REF]:

G C = G IC + (G IIC -G IC ) + (G IIIC -G IIC ) G III G II + G III G II + G III G α (1.17)
Figure 1.7: Sketch of cohesive zone model

Cohesive zone method

The pioneering work on the cohesive zone method was [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].

Along the direction of crack propagation, the cohesive zone lies between the cracked zone and intact material. It was used to describe a partially damaged zone around the crack tip observed in fracture experiments (Fig. 1.7).

In this transition zone, the stress still can be transferred through it. The separation between the fracture walls increases from zero at the crack tip to a critical value δ f when reaching the cracked zone, whereas the traction decreases from tensile strength τ c to zero. During the process of separation, the damage variable D varies linearly between zero and one and is commonly used to track the separation of the cohesive zone, i.e., from damage initiation when the traction reaches the tensile strength to fully damaged as rising to traction-free state. The evolution can be depicted by a traction-separation law as shown in Fig. 4.2. For fractures with mixed-mode damage, the effective separation δ m is defined as:

δ m = < δ n > 2 +δ 2 s + δ 2 t (1.18)
where δ n , δ s and δ t are separations along the normal and the two shear directions (Fig. 4.2), respectively. The damage variable D then can be denoted as:

D m = δ f m (δ m -δ 0 m ) δ m (δ f m -δ 0 m ) (1.19)
where δ f m represents the critical separation in the mixed mode, can be written as:

δ f m = 2G c τ e (1.20)
where G c is the equivalent fracture energy, τ e is the equivalent strength:

τ e = τ 2 n + τ 2 s + τ 2 t (1.21)
For the numerical simulation of the cohesive zone model, cohesive elements are used to model the fracture. Typically, the fracture paths should be defined beforehand, and cohesive elements are inserted along those paths that are surrounded by continuum elements. That is, here the cohesive elements and continuum elements share nodes on the intersection lines. Additionally, when dealing with fluid flow problems, intermediate nodes with pressure degrees of freedom will be added to pure mechanical cohesive elements (i.e., coupled deformation/pressure cohesive elements, see Fig. 1.8). This sort of cohesive element is able to handle both tangential flow (i.e., flow along fractures) and normal flow (i.e., leak-off). The fracture aperture w in the fluid flow equation Eq. ( 4.11) is defined as the difference between the deformed thickness and undeformed thickness:

w = h -H (1.22)
where the initial thickness H can be arbitrarily given.

In addition, cohesive elements are capable of modeling the intersection of multiple fractures.

The key point is to merge the nodes at the intersection to ensure flow continuity. As shown in Fig. 1.9, a horizontal and a vertical fracture are modeled respectively by a layer of cohesive elements which share one node 10 at their intersection.

Figure 1.9: Intersection of a horizontal fracture and a vertical fracture. Four cohesive elements 1, 2, 3, 4 share the same node 10 at the intersection.

Fault reactivation

Frictional forces act along the opposite surfaces of fractures and faults in the rock, which can lead to their sliding and therefore potentially generate earthquakes. Friction thus plays a key role in the study of induced earthquakes.

1) Amontons' law and coulomb friction law

The pioneering work on friction can be traced back to Amontons' research. Considering that a solid block is put on a planar surface with a normal force N on the top and a tangential shearing force T parallel to the contact surface. When T reaches a critical value, the block will slide, which yields the relationship between T and N as follows:

T = µN (1.23)
Dividing the contact area on both sides of eq. (1.23) yields:

τ = µσ (1.24)
Eq. (1.24) represents a state to initiate the sliding, therefore µ is called static friction coefficient; however, maintaining the sliding requires less shear stress, thus corresponding to a friction coefficient µ d < µ referred to dynamic friction coefficient.

Coulomb proposed the following equation for rock friction:

τ = µσ + C 0 (1.25)
where C 0 is called the cohesion of the rock surface.

2) Rate and state friction (RSF) model -the effect of time and velocity

Experiments showed that the friction coefficient also depends on time and velocity [START_REF] Dieterich | Modeling of rock friction: 1. experimental results and constitutive equations[END_REF][START_REF] Scholz | The role of asperity indentation and ploughing in rock friction-i: Asperity creep and stick-slip[END_REF][START_REF] Tullis | Rock friction constitutive behavior from laboratory experiments and its implications for an earthquake prediction field monitoring program[END_REF], which yields the rate and state friction law as follows:

µ = µ 0 + aln V V * + bln ψV * D c (1.26)
where µ 0 and V * are reference values for friction and velocity, a, b are material constants, D c is the critical slip distance, and ψ is referred to as the state variable and varies with time and slip.

There are several laws to describe the state variable, among them the aging law is suitable for the interseismic phase, aftershocks, afterslip, or slow earthquakes [START_REF] Dieterich | A constitutive law for rate of earthquake production and its application to earthquake clustering[END_REF][START_REF] Rubin | Earthquake nucleation on (aging) rate and state faults[END_REF], as follows:

dψ dt = 1 - V ψ D c (1.27)
In addition, the slipping law may work better for the nucleation phase or large changes in slip rate [START_REF] Bayart | Evolution of fault friction following large velocity jumps[END_REF][START_REF] Ampuero | Earthquake nucleation on rate and state faults-aging and slip laws[END_REF], as follows:

dψ dt = - V ψ D c ln V ψ D c (1.28)
For simplicity, we consider a constant friction coefficient in this study. This is based on the following reasons.

First, [START_REF] Byerlee | Friction of rocks[END_REF] found the following empirical relationships based on numerous experimental results:

   τ = 0.85σ f or σ < 200 MPa, τ = 0.6σ + 0.5 f or σ > 200 MPa (1.29)
Eq. (1.29) is often called Byerlee's law. It indicates a constant friction coefficient of 0.85 for most subsurface engineering cases. However, natural fractures/faults typically contain minerals that can greatly reduce friction.

Second, despite numerous researchers have successfully applied the rate and state friction model to the fault slip modeling during fluid injection, e.g., [START_REF] Mcclure | Investigation of injection-induced seismicity using a coupled fluid flow and rate/state friction model[END_REF][START_REF] Norbeck | Injection-triggered seismicity: An investigation of porothermoelastic effects using a rate-and-state earthquake model[END_REF], others argue differently, e.g., [START_REF] Scuderi | Frictional stability and earthquake triggering during fluid pressure stimulation of an experimental fault[END_REF] concluded that the rate and state effects may be of the second order importance during fluid injection induced fault slip experiments. In addition, they showed that an increase in fluid pressure can cause dynamic unstable slip even when the fault has velocity-strengthening behavior (it should favor aseismic creep).

Chapter 2

Hydraulic diffusivity of a partially open rough fracture

Résumé

Nous étudions l'impact de la fermeture d'une fracture ayant des surfaces rugueuses, sur la diffusivité hydraulique de la fracture, qui contrôle l'évolution spatio-temporelle des perturbations de la pression interstitielle dans les formations géologiques, en particulier celles composées d'une matrice imperméable et de fractures naturelles hautement perméables. Nous construisons des distributions d'ouvertures de fractures synthétiques à l'échelle du réservoir (∼ 500 m) à partir d'un modèle autoaffine avec des exposants de Hurst isotropes dérivés d'observations de terrain de surfaces de failles. Pour quantifier la diffusivité hydraulique des fractures rugueuses, nous effectuons des simulations par éléments finis de l'écoulement transitoire du fluide dans une seule fracture. Nous utilisons une représentation surfacique de l'ouverture de la fracture suivant l'approximation de lubrification de Reynolds. Nous vérifions que notre approximation est valide pour un écoulement en régime permanent et un faible nombre de Reynolds (Re ≪1) à partir de la comparaison avec un modèle d'ouverture de fracture représenté en volume résolu par les équations de Navier-Stokes pour les fluides incompressibles (INS). Ensuite, la diffusivité hydraulique effective de la fracture rugueuse est estimée en ajustant le champ de pression calculé à la solution d'un modèle équivalent de plaque parallèle. Les résultats montrent que le champ d'ouverture de corrélation à longue distance (jusqu'à l'échelle de la faille) dû à l'auto-affinité affecte significativement la diffusion de la pression hydraulique, ce qui se manifeste par une forte variabilité de la distribution de la pression avec l'orientation de la chute de pression imposée. Sur la base d'une rhéologie rigide-plastique, lors de la fermeture progressive de la fracture depuis le contact initial jusqu'au seuil de percolation du flux, une diminution de la diffusivité hydraulique de plus de 7 ordres de grandeur dans une direction le long de la fracture et de plus de 4 ordres de grandeur dans la direction perpendiculaire est obtenue. Nos résultats ont de fortes implications pour l'interprétation de certaines données de diffusivité hydraulique mesurées ainsi que pour l'utilisation de la diffusivité hydraulique dans l'interprétation de la distribution spatiale des événements sismiques induits par les fluides dans les réservoirs faillés.

Abstract

We investigate the impact of closing a fracture with rough surfaces on the fracture hydraulic diffusivity, which controls the spatiotemporal evolution of pore-pressure perturbations in geological formations, particularly those composed of an impermeable matrix and highly permeable natural fractures. We build distributions of synthetic fracture apertures at a reservoir scale (∼ 500 m) from a self-affine model with isotropic Hurst exponents derived from field observations of fault surfaces. To quantify the hydraulic diffusivity of rough fractures, we conduct finite element simulations of transient fluid flow in a single fracture. We use a surface representation of the fracture aperture following the Reynolds lubrication approximation. We verify that our approximation is valid for a steady-state flow and a low Reynolds number (Re ≪1) from the comparison with a volume-represented fracture aperture model solved by the Navier-Stokes equations for incompressible fluids (INS). Subsequently, the effective hydraulic diffusivity of the rough fracture is estimated by fitting the computed pressure field with the solution of an equivalent parallel plate model. The results show that the long-range correlation aperture field (up to the fault scale) due to self-affinity significantly affects hydraulic pressure diffusion, which is manifested as a strong variability in the pressure distribution with the orientation of the imposed pressure drop. Based on a rigid-plastic rheology, when closing the fracture stepwise from the initial contact to the flow percolation threshold, a decrease in the hydraulic diffusivity over 7 orders of magnitude in one direction along the fracture and over 4 orders of magnitude in the perpendicular direction is obtained. Our results have strong implications for the interpretation of some measured hydraulic diffusivity data as well as for the use of hydraulic diffusivity in interpreting the spatial distribution of fluid-induced seismic events in faulted reservoirs.

Introduction

Fluid injection into deep boreholes is often accompanied by a cluster of microseismic events [START_REF] Ellsworth | Injection-induced earthquakes[END_REF][START_REF] Shapiro | Fluid-induced seismicity[END_REF][START_REF] Cornet | Seismic and aseismic motions generated by fluid injections[END_REF][START_REF] Orlecka-Sikora | An open data infrastructure for the study of anthropogenic hazards linked to georesource exploitation[END_REF][START_REF] Cauchie | Seismic asperity size evolution during fluid injection: case study of the 1993 soultz-sous-forêts injection[END_REF]. Pore pressure diffusion through poroelastic rock is thought to be one of the primary mechanisms of fluid injection-induced seismicity since the increase in fluid pressure reduces the effective normal stress on pre-existing interfaces/faults and brings the optimally oriented interfaces/faults close to rupture [START_REF] Rice | Fault stress states, pore pressure distributions, and the weakness of the san andreas fault[END_REF][START_REF] Shapiro | Large-scale in situ permeability tensor of rocks from induced microseismicity[END_REF][START_REF] Parotidis | Back front of seismicity induced after termination of borehole fluid injection[END_REF][START_REF] Barth | Probability of earthquake occurrence and magnitude estimation in the post shut-in phase of geothermal projects[END_REF][START_REF] Blöcher | Evaluating micro-seismic events triggered by reservoir operations at the geothermal site of groß schönebeck (germany)[END_REF]. The temporal and spatial evolution of the induced microseismic events is then controlled by the hydraulic diffusivity of the rock-fracture system in the context of poroelasticity [START_REF] Shapiro | Estimating the crust permeability from fluid-injectioninduced seismic emission at the ktb site[END_REF][START_REF] Jin | Fully coupled nonlinear fluid flow and poroelasticity in arbitrarily fractured porous media: A hybrid-dimensional computational model[END_REF][START_REF] Segall | Injection-induced seismicity: Poroelastic and earthquake nucleation effects[END_REF]. The pore pressure is governed by a diffusion equation that contains the hydraulic diffusivity as the central parameter. Indeed, the linear pressure diffusion relates temporal and spatial derivatives of the fluid pressure p with the proportionality factor hydraulic diffusivity D, according to [START_REF] Jaeger | Effect of fluid pressure heterogeneity on injection-induced fracture activation[END_REF][START_REF] Wang | Theory of linear poroelasticity with applications to geomechanics and hydrogeology[END_REF][START_REF] Rozhko | Role of seepage forces on seismicity triggering[END_REF]:

∂p ∂t = D∇ 2 p. (2.1)
The hydraulic diffusivity delineates how the fluid pressure diffuses in the porous medium [START_REF] Rice | Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents[END_REF] and is an indicator of flow and transport connectivity [START_REF] Knudby | On the use of apparent hydraulic diffusivity as an indicator of connectivity[END_REF]. Indeed, it corresponds to the ratio of transport (permeability) and storage (specific storage capacity) properties that in turn depend on rock (geometrical characteristics of conduits, deformation characteristics) and fluid (viscosity, compressibility) properties. Assuming that the Biot coefficient (defined as the ratio of the fluid volume change to the bulk volume change in the drained condition) is 1, the hydraulic diffusivity is defined for fractured rock at a macroscopic scale as [START_REF] Rempe | Influence of effective stress and pore fluid pressure on fault strength and slip localization in carbonate slip zones[END_REF][START_REF] Renner | Modeling of fluid transport in geothermal research[END_REF]:

D m = k m µ(ϕc f + c pp ) , (2.2)
where k m is the intrinsic matrix permeability, µ is the dynamic viscosity of the fluid, ϕ is the porosity, and c f and c pp are the fluid compressibility and the pore space's compressibility, respectively.

Direct measurement of hydraulic diffusivity can be conducted in field tests [START_REF] Renner | Periodic pumping tests[END_REF][START_REF] Talwani | Pore pressure diffusion and the mechanism of reservoir-induced seismicity[END_REF][START_REF] Doan | In situ measurement of the hydraulic diffusivity of the active chelungpu fault, taiwan[END_REF][START_REF] Xue | Continuous permeability measurements record healing inside the wenchuan earthquake fault zone[END_REF] and laboratory experiments, such as via the pressure oscillation method [START_REF] Song | Analysis of oscillatory fluid flow through rock samples[END_REF][START_REF] Kranz | Hydraulic diffusivity measurements on laboratory rock samples using an oscillating pore pressure method[END_REF] and the pulsedecay test [START_REF] Brace | Permeability of granite under high pressure[END_REF]Hsieh et al., 1981;[START_REF] Wang | Theory of linear poroelasticity with applications to geomechanics and hydrogeology[END_REF][START_REF] Nicolas | Pore pressure pulse migration in microcracked andesite recorded with fibre optic sensors[END_REF]. In hydrogeology, Eq. (2.2) is commonly used to estimate the effective permeability of a reservoir when the hydraulic diffusivity is determined from hydraulic tests. This may work well when the rock matrix is highly permeable. In the case where fractures dominate the fluid flow, it would be better to consider the contribution of fractures separately [START_REF] Ortiz | Two-dimensional numerical investigations on the termination of bilinear flow in fractures[END_REF]. This is evident by different observed values of hydraulic diffusivity. For example, the hydraulic diffusivity of an intact rock sample is typically smaller than the hydraulic diffusivity derived from field tests where fractures exist, e.g., the values of sandstone range between 10 -6 and 10 -5 m 2 /s in laboratory measurements [START_REF] Song | Analysis of oscillatory fluid flow through rock samples[END_REF], compared with between 10 -1 and 10 0 m 2 /s in a field test [START_REF] Renner | Periodic pumping tests[END_REF].

Our interest here lies in the response of the fracture geometry to the linear fluid pressure diffusion (Eq. 2.1) when a single fracture acts as the preferential flow pathway. The linear diffusion equation can be derived from the conservation of mass in the fracture, but it requires that the pressure changes inside the fracture be sufficiently small so that the fracture deformation can be ignored, i.e., the fracture aperture remains reasonably constant [START_REF] Murphy | Semianalytical solutions for fluid flow in rock joints with pressure-dependent openings[END_REF]. This assumption is also made in other numerical modeling strategies [START_REF] Ortiz | Two-dimensional numerical investigations on the termination of bilinear flow in fractures[END_REF][START_REF] Vinci | The imprint of hydro-mechanics of fractures in periodic pumping tests[END_REF]. Previous studies also show that a small pressure disturbance (< 0.1 MPa) is able to trigger seismicity [START_REF] Keranen | Sharp increase in central oklahoma seismicity since 2008 induced by massive wastewater injection[END_REF][START_REF] Dempsey | Response of induced seismicity to injection rate reduction: Models of delay, decay, quiescence, recovery, and oklahoma[END_REF][START_REF] Goebel | The 2016 mw5. 1 fairview, oklahoma earthquakes: Evidence for long-range poroelastic triggering at> 40 km from fluid disposal wells[END_REF], particularly for stress-critical faults and during a post-induced seismicity or aftershock stage [START_REF] Schmittbuhl | Induced and triggered seismicity below the city of strasbourg, france from november 2019 to january 2021[END_REF][START_REF] Noir | Fluid flow triggered migration of events in the 1989 dobi earthquake sequence of central afar[END_REF][START_REF] Nur | Aftershocks caused by pore fluid flow?[END_REF]. It is also relevant for some EGS reservoirs or fractured media where the fluid pressure is close to hydrostatic conditions (drained conditions) and the fluid volume is connected over a long distance to the surface, so that the fluid pressure is significantly lower than the solid counterpart (trace of the stress tensor). In the framework of linear pressure diffusion, the fracture hydraulic diffusivity can be approximated as [START_REF] Murphy | Semianalytical solutions for fluid flow in rock joints with pressure-dependent openings[END_REF]:

D f = h 2 12µ(c f + h -1 c j ) , (2.3)
where h is the fracture aperture and c j is the fracture compressibility.

Eq. (2.3) is suitable for a parallel plate fracture where the permeability is calculated by the cubic law. For real rough fractures/faults, the deviation from the cubic law may be considerable [START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF][START_REF] Almakari | Effect of the injection scenario on the rate and magnitude content of injection-induced seismicity: Case of a heterogeneous fault[END_REF]Ji et al., 2020). More specifically, compared to a single parallel plate fracture with an identical mean aperture, roughness can either enhance or inhibit fluid flow [START_REF] Méheust | Flow enhancement of a rough fracture[END_REF][START_REF] Schmittbuhl | Fracture morphology and viscous transport[END_REF][START_REF] Guo | Thermal drawdown-induced flow channeling in a single fracture in egs[END_REF]. Moreover, the permeability of a self-affine fracture shows a certain degree of anisotropy when the orientation of the imposed pressure drop is changed [START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF].

To numerically compute the permeability of a rough or partially open fracture, the common approach is to apply Darcy's law under laminar flow conditions, although it should only be considered as a qualitative measurement [START_REF] Blöcher | Permeability of matrix-fracture systems under mechanical loading-constraints from laboratory experiments and 3-d numerical modelling[END_REF]. This method leads to replacement of the fracture aperture h by the hydraulic aperture d h in Eq. (2.3). The hydraulic aperture of a rough fracture is then classically introduced as an effective measure of the hydraulic performance using the directional total flux and pressure difference [START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF][START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF][START_REF] Murphy | Semianalytical solutions for fluid flow in rock joints with pressure-dependent openings[END_REF]Neuville et al., 2010b) and is therefore different from the arithmetic mean fracture aperture. Once the permeability is known, it is possible to estimate the diffusivity from Eq. (2.3). However, this estimation may seem somehow rough for a partially open rough fracture since the permeability calculation may not precisely describe the pressure propagation, i.e., the temporal and spatial pressure evolution.

In this work, we propose a new approach to determine the fracture hydraulic diffusivity. Instead of exploring the relationship between the diffusivity and those intrinsic hydraulic parameters, here, we focus on a direct quantification of the spatiotemporal evolution of the pressure inside the fracture. First, by solving the linear diffusion equation for a rough fracture with given initial and boundary conditions, we obtain the pressure profiles as a function of time and space. The pressure profiles are then compared to those derived from an analytical solution in which hydraulic diffusivity serves as an unknown. Next, we use the least square regression to search for a diffusivity that best matches the two pressure profiles. This diffusivity can then be seen as an effective hydraulic diffusivity of a rough fracture. Finally, we compare our results with the hydraulic diffusivity values obtained from Eq. (2.3).

Based on the proposed approach, the main objective of this study is to quantify the impact of fracture roughness and fracture closure on hydraulic diffusivity.

We use a rigid deformation to mimic fracture closure by interpenetrating the two fracture halves stepwise into each other. By interpenetrating the two fracture halves, we generate a partial overlap of the two volumes. During laboratory experiments (Kluge et al., 2021b), brittle deformation in these contact/overlap areas producing fines was observed. This plastic deformation is mimicked by removing the overlapping volume and setting the aperture to zero. Therefore, we call the simulated fracture closure process a "rigid-plastic" deformation. In the simulation, we close the fracture in a stepwise manner, forming different contact areas, and we compute the effective hydraulic diffusivity at each step. This paper is organized as follows. In Sect. 2.2, the generation of a self-affine fracture aperture distribution based on field observations is introduced. In Sect. 2.3, the mathematical formulation for fluid flow through a single fracture is given for the two cases of interest for this study, that is, fluid flow at low Reynolds numbers under either steady state or transient conditions. This is followed by Sect. 2.4, where the modeling results and the analysis of the influence of roughness and fracture closure on the effective hydraulic diffusivity are presented. In Sect. 2.5, the anisotropy of the effective hydraulic diffusivity is further discussed, and the results are compared to those available in the literature. The paper ends with a brief conclusion in Sect. 2.6. In addition, appendices regarding the fluid velocity calculation and the effective hydraulic diffusivity estimation are provided in detail at the end of the paper. A brief workflow of the numerical procedure is given in Fig. 2.1.

Fracture aperture generation

In our modeling, the aperture field of a partially open fracture between two opposite rough surfaces is built in four steps: (1) generate a single fracture surface from a self-affine surface generator; (2) build a self-affine aperture by mirroring the generated surface; (3) stepwise close the fracture based on an imposed normal displacement field and by assuming a perfectly rigid plastic rheology of the asperities; and (4) finite element (FE) mesh the fracture of the generated geometries.

Step 1: generation of a self-affine fracture surface

It has been shown that fresh surfaces of lab-scale samples can be well described by self-affinity (Schmittbuhl et al., 1993a;Zimmerman et al., 2004;Neuville et al., 2012b). The roughness of a two dimensional (2D) self-affine profile is statistically invariant under the following scaling transformation (Schmittbuhl et al., 1995a):

δx → λδx δz → λ H δz , (2.4)
where δx is the coordinate along the profile, δz is the vertical direction amplitude, λ is a scaling factor that can take any positive value, and H (0 < H < 1) is called the roughness exponent or Hurst exponent. The magnified profile is statistically similar to the original profile, as shown in Fig. 2.2. [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF].

Scanning measurements of the surface roughness of a set of faults (Fig. 2.3) also reveal a self-affine behavior over nine decades of length scales (i.e., from 50 µm to 50 km) [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF][START_REF] Renard | High resolution 3d laser scanner measurements of a strike-slip fault quantify its morphological anisotropy at all scales[END_REF][START_REF] Candela | Characterization of fault roughness at various scales: Implications of three-dimensional high resolution topography measurements[END_REF]. The data gap in Fig. 2.3 can be supplemented by the surface roughness measurement of the fault from the Gole Larghe Fault Zone in the range of 0.5 mm -500 m [START_REF] Bistacchi | Fault roughness at seismogenic depths from lidar and photogrammetric analysis[END_REF]. At scales ranging from 50 µm to 10 m, when surfaces in contact experience significant slip, they exhibit anisotropy with the Hurst exponent H ∥ ∼ 0.6 in the slip direction and H ⊥ ∼ 0.8 perpendicular to it [START_REF] Candela | Characterization of fault roughness at various scales: Implications of three-dimensional high resolution topography measurements[END_REF]. For rupture traces on scales of 200 m to 50 km, the self-affinity is isotropic and consistent with the slip-perpendicular behavior of the smaller-scale measurements (i.e., H ∼ 0.8) [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF].

Figure 2.3: Log-log graph of Fourier power spectra C(q) of linear profiles as a function of the wavenumber q for our synthetic self-affine surface (large blue dots) and for a group of natural fault surfaces (small dots) from [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. The synthetic self-affine surface is generated with H = 0.8 and h rms = 0.1 m at the 512 m scale to match field observations

To characterize the roughness of the self-affine surface from topographic measurements, the Fourier power spectrum C(q) (i.e., the square of the modulus of the Fourier transform) is introduced as a function of the wavenumber q (Schmittbuhl et al., 1995b;Zimmerman et al., 2004;[START_REF] Candela | Characterization of fault roughness at various scales: Implications of three-dimensional high resolution topography measurements[END_REF]. The power spectrum of linear profiles vs. the wavenumber show a linear trend in a log-log plot (Fig. 2.3) as follows:

C(q) = ζq -1-2H , (2.5)
where ζ represents the intercept of the power law line, the so-called 'prefactor' of the power spectrum [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF], determining the roughness amplitude h rms at a given scale, i.e., root mean square (RMS) average. The Hurst exponent H is then directly linked to the power-law slope.

To generate a 2D self-affine surface, a Gaussian random field is first created. The spatial correlation is introduced through the Fourier transform of the random field. Next, by using inverse Fourier transform, the distribution of the surface can be finally obtained [START_REF] Candela | Characterization of fault roughness at various scales: Implications of three-dimensional high resolution topography measurements[END_REF]. The roughness amplitude h rms is introduced in order to adjust the roughness in the z-direction by normalizing the heights of the discrete points on the 2D grid to obtain a prescribed RMS of the whole height.

As an example, a synthetic self-affine surface is generated with an isotropic roughness exponent H ∥ = H ⊥ = 0.8 at the 512 m scale and h rms = 0.1 m. The scaling of the power spectral density of our synthetic fault surface (blue dots in Fig. 2.3) is consistent with that of a collection of major faults from field observations [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. The height probability follows a Gaussian distribution with the prescribed h rms (Fig.

Step 2: generation of a self-affine aperture distribution

The fracture aperture is defined as the distance between the two facing surfaces, perpendicular to the nominal fracture plane [START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF][START_REF] Méheust | Scale effects related to flow in rough fractures[END_REF]Neuville et al., 2010b;[START_REF] Marchand | A stochastic study of flow anisotropy and channelling in open rough fractures[END_REF]. If the fracture is modeled by two parallel plates, then the aperture is constant. However, for real rock fractures with rough walls, the aperture varies in space, resulting in an aperture distribution h(x, y) [START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF][START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF]. To obtain h(x, y), we first reproduce the generated self-affine surface from the previous section and flip it vertically (along z). We then place the two surfaces such that they face each other with a separation of d m . Because the two facing surfaces are unmated and self-affine, the resulting aperture field h(x, y) is also self-affine and shares the same Hurst exponent (Neuville et al., 2010b). Consequently, the aperture distribution h(x, y) can be alternatively generated as a single self-affine object z(x, y) plus the average aperture d m . Accordingly, the upper and lower boundaries of the aperture field can be written as:

s 1 (x, y) = + z(x, y) + d m 2 , (2.6) s 2 (x, y) = - z(x, y) + d m 2 , (2.7)
where the Fourier power spectra of s 1 (x, y) and s 2 (x, y) should correspond to Fig. 2.3. Fig. 2.5(a) shows a sketch (i.e., vertical exaggeration) of a 2D cross-section of the self-affine aperture model h(x, y) from two unmated symmetrical surfaces s 1 (x, y) and s 2 (x, y) when d m is larger than the maximum height. In this case, it is a fully open fracture where the two surfaces are completely separated.

Step 3: fracture closure

When a normal displacement is imposed on the fracture, i.e., to close the fracture along its normal direction, there are regions where the two opposing faces of the fracture wall virtually overlap each other (Fig. 2.5b). These regions are called contact areas and commonly exist in real rock fractures [START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF]. In this case, the mean aperture d m computed over the whole surface assuming a zero aperture at the contact areas will be smaller than the mean aperture of the open areas occupied by the fluid. In this study, since we consider both fractures with and without contact, we hereinafter refer to d m as the mean aperture of the open area in both cases. For the aperture of overlapping asperities, we set it to zero assuming that the solid part is perfectly plastic. The plastic limit is defined as the strength of the material over which the material will be eroded This is effectively different from an elasto-plastic model. If the stress concentration on asperities is very high, this approximation is relevant. This contact model is based on the interpenetration approach, which is simple but fast and effective [START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF][START_REF] Méheust | Scale effects related to flow in rough fractures[END_REF][START_REF] Pei | Finite element modeling of elasto-plastic contact between rough surfaces[END_REF][START_REF] Watanabe | Determination of aperture structure and fluid flow in a rock fracture by high-resolution numerical modeling on the basis of a flow-through experiment under confining pressure[END_REF][START_REF] Liu | Normal-stress dependence of fracture hydraulic properties including two-phase flow properties[END_REF][START_REF] Kling | Numerical simulations and validation of contact mechanics in a granodiorite fracture[END_REF].

Step 4: finite element meshing

The three-dimensional (3D) fracture aperture distribution h(x, y) is a spatial function of the x and y coordinates along the mean plane of the fracture. However, when it is applied to the flow simulation, the requirement of fine isotropic meshes in the 3D volume may lead to difficulties, especially when the scale discrepancy between the aperture (along z) and the fracture size (along x and y) is large.

For instance, if we assume that the fracture fluctuations are on the order of h rms and the length is L, then the aspect ratio of the fracture is h rms /L. If the surface is self-affine, then h rms ∼ L H . Consequently, the aspect ratio can scale as L H-1 . If H < 1 and L becomes very large, then the aspect ratio goes to zero [START_REF] Schmittbuhl | Fracture morphology and viscous transport[END_REF]. In this case, the 3D isotropic mesh requires very fine elements due to the small scale in the z-direction, leading to a large number of elements for large fractures. However, this fine meshing is unnecessary along the x-axis and y-axis since aperture variations along x and y are on a larger length scale. Accordingly, anisotropic elements for a 3D volume are required. Therefore, we use prism elements with a mean aspect ratio of 1:25 for the volume representation. To reduce the numerical cost for large-scale reservoir simulation, we can approximate the volumetric representation of the 3D fracture by a 2D surface representation, as shown in Fig. 2.6. A validation of this 3D to 2D approximation is discussed later in the paper. For a 2D surface representation, we make use of quadrilateral mesh elements. In both cases, the meshing stage is done by relying on the open source software Gmsh [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in preand post-processing facilities[END_REF].

From the Navier-Stokes equation to the Reynolds equation

The fluid flow of an incompressible Newtonian fluid is governed by the Navier-Stokes equation, which expresses momentum and mass conservation over the fracture void space as follows [START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF][START_REF] Guyon | Physical hydrodynamics[END_REF]:

∂u ∂t + (u • ∇)u = - 1 ρ ∇p + µ ρ ∇ 2 u, (2.8) ∇ • u = 0, (2.9)
where u is the fluid velocity vector, p denotes the fluid pressure, and ρ and µ are the fluid density and dynamic viscosity, respectively.

The inertial forces represented by the nonlinear term (u • ∇)u in Eq. (2.8) render these equations difficult to solve. This nonlinear inertial term can be neglected if the viscous forces dominate over the inertial forces, that is, for small Reynolds numbers, i.e., Re ≪ 1 [START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF]. The Reynolds number is defined as the ratio of the inertial terms to the viscous terms:

Re = ρU l 2 z µl h , (2.10)
where U is the magnitude of the horizontal average velocity and l z and l h are the scales of vertical and horizontal velocity variations, respectively.

Considering a steady fluid flow in the low Reynolds number regime, the N-S equation can be linearized to [START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF]:

∇p -µ∇ 2 u = 0.
(2.11)

For a rough fracture, assuming that the local aperture h(x, y) variations satisfy the lubrication approximation (i.e., ||∇h ≪ 1||), Eq. (2.11) can be recast in the form of the Reynolds equation as [START_REF] Zimmerman | Fluid flow in rock fractures: From the navier-stokes equations to the cubic law[END_REF]:

∇ • h(x, y) 3 12µ ∇p = 0.
(2.12)

The Reynolds equation has been applied to solve fluid flow through rough-walled fractures by many authors [START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF][START_REF] Renshaw | On the relationship between mechanical and hydraulic apertures in roughwalled fractures[END_REF][START_REF] Brush | Fluid flow in synthetic rough-walled fractures: Navier-stokes, stokes, and local cubic law simulations[END_REF][START_REF] Marchand | A stochastic study of flow anisotropy and channelling in open rough fractures[END_REF], with the local cubic law for computing the flux [START_REF] Oron | Flow in rock fractures: The local cubic law assumption reexamined[END_REF][START_REF] Klimczak | Cubic law with aperture-length correlation: implications for network scale fluid flow[END_REF].

The pressure diffusion equation

Assuming the validity of the lubrication approximation, a surface-represented fracture (Fig. 2.6b) can be parameterized in terms of an effective aperture. Within the fracture plane, the fluid pressure is governed by the following equation [START_REF] Cacace | Flexible parallel implicit modelling of coupled thermalhydraulic-mechanical processes in fractured rocks[END_REF]:

h K f ∂p ∂t + ∇ • (hq) = 0, (2.13)
where K f is the fluid bulk modulus. The fluid flow is computed by Darcy's law in the local coordinate system as follows:

q = - k f rac µ • ∇p, (2.14)
where k f rac denotes the fracture permeability tensor. Assuming the local cubic law is valid in the laminar flow within the fracture plane, the isotropic permeability in the local coordinate system can be identified as:

k f rac = h 2 12 1, (2.15)
where 1 is the unit tensor.

Finite element modeling using the MOOSE/GOLEM framework

Steady flow through a rough fracture -a validation of the volume representation vs. surface representation

The incompressible Navier-Stokes equations are applied to solve for the dynamic pressure and the velocity field inside the fracture opening. These equations correspond exactly to the Navier-Stokes equations in 'Laplace' form in the INS Module of MOOSE [START_REF] Peterson | Overview of the incompressible navier-stokes simulation capabilities in the moose framework[END_REF]. A stabilized Petrov-Galerkin finite element method is used to solve these equations with appropriate initial and boundary conditions.

For a self-affine fracture aperture model with Hurst exponent H, the Reynolds lubrication approximation is [START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF]:

l c ≫ ∆h(L) L H 1 1-H , (2.16)
where l c denotes the lower bound of self-affine scaling and ∆h(L) is the maximum distance between the rough walls at a scale L.

In this study, we consider a uniform mesh size of 1 m. When using l c = 1 m and H = 0.8, for fracture length L = 64 m, the maximum ∆h(L) is ≈ 0.15 m, and we have l c ≫ 4.5 × 10 -12 m. For L = 512 m, ∆h(L) ≈ 0.6 m, and we have l c ≫ 1.13 × 10 -12 m. Based on Eq. (2.16), we can conclude that our model satisfies the Reynolds lubrication approximation.

When the fluid flow reaches the steady state, the fluid pressure governing equation of the surface-represented fracture is equivalent to the Reynolds equation (Eq. 2.12). To test the degree of validity of this approximation, for our transient flow model, we take the following steps. First, a self-affine fracture surface is synthetically generated at an intermediate scale of 64 m with isotropic Hurst exponent H = 0.8 and roughness amplitude h rms = 0.025 m. A set of apertures is then created with the ratio of the mean aperture to the roughness amplitude d m /h rms in the range of 1 to 15. Among them, the fracture remains fully open when the ratio is greater than or equal to 3, whereas fracture-wall contact occurs when d m /h rms < 3. The ratio d m /h rms = 1 corresponds to ∼ 20% of the contact area. Next, we build the volume-represented and surface-represented fractures for each aperture distribution, along with one of their finite element meshes, which are shown in Fig. 2.6a and Fig. 2.6b, respectively. Finally, fluid flow simulations are conducted for both the surface and volume representations of the fracture. A pressure difference ∆p = 10 -8 Pa between the two ends of the model along either the x-axis or the y-axis is set as the boundary conditions (i.e., 10 -8 Pa and 0 Pa pressure for the inlet and outlet boundaries, respectively; see Fig. 2.7).

To compare the results of the two types of simulations, we introduce the following definition of the hydraulic aperture d h :

d h = 12 V µ |∆p| 1/3
.

(2.17)

The total volume flux V along the flow of the rough fracture is obtained from the mean of the local flux q(x, y). Considering that the pressure drop is along the x-direction, q(x, y) is then the product of the velocity profile v x and the local aperture h(x, y).

For the volume-represented fracture case, the velocity has a z component. We then average the velocity over z such that v x = v x (x, y, z). Since the local cubic law is assumed to be valid, the velocity locally follows a parabolic shape along the z-direction, as in the parallel plate model. The estimation of v x (x, y, z) is given in Appendix A. For the surface-represented fracture case, we have v x = v x (x, y), which can be directly obtained from the velocity profile.

Transient flow in a parallel plate configuration -a reference model

In this section, we check the transient flow in a parallel plate model with an analytical solution. Consider an instantaneous increase in the pressure difference along the x-axis from 0 to L. The initial pore pressure p 0 is given at t = 0. During the whole diffusive process, a fixed pressure p 1 is applied at x = 0, whereas at the end of the plate (x = L), the pressure is maintained constant as per its initial value (p 0 ):

p = p 0 at t = 0 f or 0 ≤ x ≤ L p = p 1 at x = 0 f or t > 0 p = p 0 at x = L f or t > 0 .
(2.18)

Based on Eq. (2.18), the 1D solution of the linear diffusion equation Eq. (2.1) is given by [START_REF] Turcotte | Geodynamics[END_REF][START_REF] Carlsaw | Conduction of heat in solids[END_REF]:

p(x, t) = p 1 + (p 0 -p 1 ) x L + 2 π ∞ n=1 1 n exp - Dn 2 π 2 t L 2 sin nπx L . (2.19)
A characteristic time t * is defined as:

t * = L 2 /D.
(2.20)

When t ≫ t * , the pressure reaches an equilibrium linear profile, and the solution becomes:

p = p 1 + (p 0 -p 1 ) x L .
(2.21)

A simulation is performed on the parallel plate fracture with dimensions of L × L (512 m × 512 m) and an aperture of 0.1 m. The inlet pressure p 1 and outlet pressure p 0 are set to 10 -5 Pa and 0 Pa along the x-axis, respectively. This small pressure gradient is chosen to separate the effect of fracture deformation on diffusivity. Therefore, a direct relation of the effects of fracture roughness and stepwise fracture closure on diffusivity can be obtained. Regarding the fluid properties, the fluid modulus and the fluid dynamic viscosity are 2.2 GPa and 0.001 Pa•s, respectively. In this case, fracture apertures are represented as a surface and modeled by single-node interfacial elements [START_REF] Jacquey | Modelling coupled fluid flow and heat transfer in fractured reservoirs: description of a 3d benchmark numerical case[END_REF].

On the basis of the settings in the last paragraph, we obtain the fracture's hydraulic diffusivity D f = 1.833 × 10 9 m 2 /s (Eq. 2.13) and the corresponding characteristic time t * = 1.43 × 10 -4 s (Eq. 2.20). We then substitute D f = 1.833 × 10 9 m 2 /s and L = 512 m into Eq. (2.19) to obtain the pressure distribution p a (x, t). On the other hand, the numerically computed pressure distribution along the parallel plate fracture is p ∥ (x, y, t), which is then averaged along the y-axis as < p ∥ (x, t) > y . We then plot both p a (x, t) and < p ∥ (x, t) > y at t 1 = 5 × 10 -6 s, t 2 = 1.5 × 10 -5 s and t 3 = 4.5 × 10 -5 s in Fig. 2.9. The consistency between p a (x, t) and < p ∥ (x, t) > y illustrates that the numerical result p a (x, t) can be seen as the analytical solution for the parallel plate model averaged over the transverse direction, providing a reference for the rough fracture study.

Figure 2.9: Comparison between the simulation results for the parallel plate model (red squares) and analytical solutions (solid lines) of pressure diffusion at t 1 = 5 × 10 -6 s, t 2 = 1.5 × 10 -5 s and t 3 = 4.5 × 10 -5 s

Results

Based on the results described in the previous paragraph, we hereafter consider only a 2D surface representation of the rough fracture to quantify the effective hydraulic diffusivity (and its variations) under transient flow conditions.

Spatiotemporal evolution of the pressure in a large-scale fracture

The self-affine aperture distributions are generated with L = 512 m, isotropic H = 0.8 and h rms = 0.1 m. When the mean aperture d m varies, the fracture remains fully or partially open. As an example, Fig. 2.10 illustrates the aperture map for d m = 0.3 m, i.e., d m /h rms = 3, corresponding to a totally open fracture but close to contact of the asperities on the opposite fracture walls. Given this size and the mean aperture of the fracture, the characteristic time of the equivalent parallel plate model is t * = 1.5 × 10 -5 s (Eq. 2.20).

For transient fluid flow simulations, all parameters used, initial conditions and boundary con- ditions are consistent with Sect. 2.3.3.2. The time step is set small enough to capture the transient pressure diffusion process (e.g., 1 × 10 -7 s in the case of Fig. 2.10). To investigate the anisotropy, simulations are performed with an imposed macroscopic pressure gradient along either the x-or along the y-direction.

The pressure distributions and the fluid flux distributions of the rough fracture with d m = 0.3 m at t 1 = 5 × 10 -7 s, t 2 = 2 × 10 -6 s and t 3 = 5 × 10 -6 s are presented in Fig. 2.11 and Fig. 2.12, respectively. For comparison, the pressure and flux distribution for a parallel plate model with the same d m aperture are also shown (Fig. 2.11 and 2.12a, b, c). As expected, the roughness causes deviations of the pressure diffusion from the parallel plate model. Visually, the pressure diffuses faster in the y-direction (Fig. 2.11g, h, i) than in the x-direction (Fig. 2.11d, e, f), where the pressure front is propagating significantly slower than that in the parallel plate configuration with the same mean aperture (Fig. 2.11a, b, c). As fluid flows along the fracture, it follows preferential pathways (i.e., channels) due to the impact of aperture variations (Fig. 2.12). Moreover, compared to the xdirection (Fig. 2.12d, e, f), a stronger channeling effect is observed in the y-direction (Fig. 2

.12g, h, i). (a) t 1 = 5 × 10 -7 s (b) t 2 = 2 × 10 -6 s (c) t 3 = 5 × 10 -6 s (d) t 1 = 5 × 10 -7 s (e) t 2 = 2 × 10 -6 s (f) t 3 = 5 × 10 -6 s (g) t 1 = 5 × 10 -7 s (h) t 2 = 2 × 10 -6 s (i) t 3 = 5 × 10 -6 s

Effective hydraulic diffusivity

To quantify the pressure diffusion along the rough fracture, an effective hydraulic diffusivity D e is obtained by fitting the pressure solution in time and space for the rough aperture by the parallel plate solution. The approach is similar to the assessment of the hydraulic aperture that is defined by fitting the effective hydraulic flux of a rough fracture by a parallel plate model. The procedure is as follows.

The numerical pressure distribution p R (x, y, t) for the rough aperture is first averaged along the yaxis as < p R (x, t) > y . We then optimize the hydraulic diffusivity of a parallel plate model with a pressure distribution < p ∥ (x, t, D f ) > y to match < p R (x, t) > y in the least square error sense.

Noting the consistency between < p ∥ (x, t, D f ) > y and the analytical solution < p a (x, t, D f ) >, we have:

< p R (x, t, D e ) > y =< p ∥ (x, t, D f ) > y = p a (x, t, D f ). (2.22)
Details of the procedure are given in Appendix B.

Fig. 2.13 and Fig. 2.14 illustrate how the effective diffusivity D e is different from the hydraulic diffusivity of the parallel plate model D m , which has the same mean aperture d m as the rough fracture. The comparison between the parallel plate model and the rough fracture is shown in terms of the pressure distribution at t 1 = 5×10 -7 s, t 2 = 2×10 -6 s and t 3 = 5×10 -6 s and along the x-and y-axes. When they have the same mean aperture (Fig. 2.13a, Fig. 2.14a), the pressure diffusion of the rough fracture significantly deviates from that in the parallel plate model. In contrast, when they have the same hydraulic diffusivity (Fig. 2.13b, Fig. 2.14b), the pressure diffusion curves match well with some slight discrepancies at some positions, which demonstrates that the effective hydraulic diffusivity reflects the rate of the pressure diffusion along the rough fracture as a whole.

Effect of fracture closure

In this section, we consider normal closure in the direction perpendicular to the mean plane of the fracture, which leads to fracture surfaces contacting each other. Closure is obtained by imposing a normal displacement stepwise along the whole open fracture surface. Under the perfect plastic assumption, owing to the self-affine property of the aperture, an increase in normal displacement leads to a decrease in the mean aperture, which follows a linear trend with the increase in the contact area when the closure becomes significant, as shown in Fig. 2 is relatively large (relatively small contact), there is a linear behavior between the two quantities:

d m = 0.863d h + 0.053.
Interestingly, the slope is not one, showing that the mean aperture d m is decreasing slower than the hydraulic aperture d h . Additionally, there exists a residual mean aperture at zero hydraulic aperture, showing that immobile fluid is trapped at the percolation threshold.

When approaching the percolation threshold, the decrease rate of d h is faster than that of d m . This behavior is attributed to the strong increase in tortuosity and channeling of the flow as the contact area increases [START_REF] Nolte | The fractal geometry of flow paths in natural fractures in rock and the approach to percolation[END_REF][START_REF] Unger | Numerical study of the hydromechanical behavior of two rough fracture surfaces in contact[END_REF][START_REF] Sahimi | Flow and transport in porous media and fractured rock: from classical methods to modern approaches[END_REF]. For the simulations, we set the pressure drop either along the x-axis or the y-axis. We increase the normal displacement step by step until the fracture aperture reaches the percolation threshold (i.e., loss of the hydraulic connection from the inlet to the outlet and zero fluid velocities) in the two directions. We calculate the hydraulic aperture d h and the effective hydraulic diffusivity D e for each stage.

The effective diffusivity D e is plotted as a function of the hydraulic aperture d h in log-log space (Fig. 2.16). As a check, we plot the hydraulic diffusivity calculated by Eq. (2.3), which is suitable for the parallel plate model (i.e., D m , dashed line in Fig. 2.16). It shows that the effective hydraulic diffusivity D e is close to the prediction from Eq. (2.3) only when the mean aperture is large (i.e., at a relatively low contact area). As d m decreases, the decrease in effective diffusivity is either faster or slower than that in the diffusivity of the equivalent parallel plate model according to the orientation of the imposed pressure drop.

In the x-direction, the effective diffusivity drops by 7 orders of magnitude from ∼ 10 9 m 2 /s without contact to ∼ 10 2 m 2 /s when it is close to the percolation threshold (∼ 66%). By contrast, we obtain a 4 orders of magnitude reduction in the y-direction, but with a lower percolation threshold (∼ 31%). Fig. 2.17 shows the pressure distribution as a function of contact area in the x-direction. An increase in contact area causes an increase in tortuosity and channeling of the flow field in the aperture distribution (Fig. 2.18). Consequently, a longer time is required to reach the steady state, resulting in a decrease in the hydraulic diffusivity. When the contact area exceeds 50%, there are still large diffusivities since large channels exist in the aperture distribution. When the contact area approaches the percolation threshold, these channels are drastically reduced, and the diffusivity shows more significant changes (only a small change in contact area can lead to a large decrease in diffusivity).

Discussion

Anisotropy of hydraulic diffusivity

When the mean aperture of our fracture is large enough such that the largest asperities are about to touch each other (d m = 0.3 m), the effective hydraulic diffusivity is D e = 0.423D m and D e = 1.691D m with the imposed pressure drop in the x-direction (Fig. 2.13) and in the y-direction (Fig. 2.14), respectively, which shows anisotropic behavior (Fig. 2.16). At very large apertures ). Interestingly, when closing the fracture, the anisotropy, defined here as the ratio of the effective diffusivities for the pressure drop along the x-or y-direction, is maximum when d m /h rms ≈ 3, which is approximately when the (few) asperity contacts start to develop.

Closing the fracture further reduces the diffusivity in both directions. The decrease in diffusivity is more accentuated the closer the system is to the flow percolation threshold (Fig. 2.16). However, this decrease is different for the two pressure drop directions. At d m = 0.002 m, the difference between the y-and x-directions reaches two orders of magnitude (Fig. 2.16). In other words, the anisotropy increases when closing the fracture. Anisotropy in the hydraulic behavior depends on the geometrical heterogeneity as well as on the self-affinity of the fracture surfaces/apertures [START_REF] Méheust | Flow enhancement of a rough fracture[END_REF]. The roughness exponent (self-affinity) introduces spatial correlations to the roughness amplitude and therefore to the aperture distribution. Such long-range correlations (up to the fracture scale) of self-affine apertures induce strong channeling of the flow [START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF][START_REF] Neuville | Hydraulic transmissivity and heat exchange efficiency of open fractures: a model based on lowpass filtered apertures[END_REF].

Spatially correlated fractures tend to have only a few dominant flow paths compared to uncorrelated fractures (Pyrak- [START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF]. Although the aperture variation in x and y is statistically isotropic, the resulting aperture distribution is heterogeneous (Fig. 2.10). This leads to different flow channels along the x-and y-directions and is therefore responsible for the anisotropy of the fluid flow. As the fracture closes, the channeling effect becomes more prominent. Accordingly, the anisotropy becomes more noticeable (Fig. 2.16). The anisotropic flow behavior has been verified by lab experiments [START_REF] Méheust | Flow enhancement of a rough fracture[END_REF] and numerical studies [START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF][START_REF] Marchand | A stochastic study of flow anisotropy and channelling in open rough fractures[END_REF], both targeting self-affine surfaces with an isotropic Hurst exponent

H = 0.8.
In our study, we observed that the decrease in the hydraulic diffusivity is enhanced along the y-direction and inhibited along the x-direction. This is specific to the chosen surface, i.e., choice of the seed used to generate a random number in the generator of the self-affine surface [START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF][START_REF] Candela | Characterization of fault roughness at various scales: Implications of three-dimensional high resolution topography measurements[END_REF]. In Fig. 2.20, the behavior for two other choices of the seed while keeping the Hurst exponent H = 0.8 and the RMS h rms = 0.1 m are shown. This illustrates the variability of the behavior within the same general trend: beginning of departure from the parallel plate model for d m /h rms ≈ 3 and a strong drop in the diffusivity when approaching the percolation threshold. However, the specific sensitivity to the pressure drop orientation is different for the different orthogonal directions. The anisotropy of the effective diffusivity has an important influence on the resulting pore pressure diffusion. Some authors found that isotropic diffusivity poorly describes pressure diffu-sion (seismicity migration) compared to anisotropic diffusivities. For instance, [START_REF] Noir | Fluid flow triggered migration of events in the 1989 dobi earthquake sequence of central afar[END_REF] estimated an isotropic diffusivity of 1.2 × 10 4 m 2 /s, whereas the anisotropic hydraulic diffusivities were D xx ∼ 3 × 10 4 m 2 /s, D yy ∼ 3 × 10 3 m 2 /s and D zz ∼ 3 × 10 3 m 2 /s during the 1989 Dobi earthquake sequence, showing that the fastest seismic migration was along the x-direction. Similarly, [START_REF] Antonioli | Fluid flow and seismicity pattern: Evidence from the 1997 umbria-marche (central italy) seismic sequence[END_REF] obtained an ∼90 m 2 /s isotropic diffusivity, while the maximum value of the anisotropic diffusivity was 275 m 2 /s from the 1997 Umbria-Marche seismic sequence. The maximum diffusivity direction coincides with the strike of the active faults. The author concluded that the large diffusivity is associated with high permeability rough fractures within the damage zone of the active fault system, which essentially supports this study. In turn, if the anisotropic fracture diffusivity can be predicted properly, then it is possible to determine the orientation of the preferential earthquake migration direction.

Comparison with hydraulic measurements and implications

Our results under large closure (hydraulic diffusivities are on the order of 10 2 m 2 /s -10 4 m 2 /s) are consistent with the values derived from the analysis of some earthquake sequences [START_REF] Noir | Fluid flow triggered migration of events in the 1989 dobi earthquake sequence of central afar[END_REF][START_REF] Antonioli | Fluid flow and seismicity pattern: Evidence from the 1997 umbria-marche (central italy) seismic sequence[END_REF][START_REF] Malagnini | Control of pore fluid pressure diffusion on fault failure mode: Insights from the 2009 l'aquila seismic sequence[END_REF][START_REF] Dempsey | Response of induced seismicity to injection rate reduction: Models of delay, decay, quiescence, recovery, and oklahoma[END_REF]. These earthquakes were assumed to be triggered by the diffusion of pore pressure perturbations in a fractured medium, and the seismicity migration was then evidenced to be compatible with pore pressure relaxation. The hydraulic diffusivity estimated by [START_REF] Noir | Fluid flow triggered migration of events in the 1989 dobi earthquake sequence of central afar[END_REF] for the 1989 Dobi earthquake sequence of Central Afar ranges between 10 3 and 10 4 m 2 /s, which corresponds to a characteristic width (i.e., effective aperture) of 1 mm -3 cm. The consistency with our results indicates that our model might be used to predict potential earthquake migration, particularly when a single fault path dominates the fluid flow. Compared to diffusivities estimated from direct hydraulic tests, the values obtained from our simulations are somewhat large. The discrepancy could be attributed to several aspects. First, there is an issue regarding the representative elementary volume (REV) of the measurement. For example, in the laboratory, the tested target is typically an intact rock sample, where fluid flow is restricted by interconnected pores, resulting in small diffusivities generally ranging between 10 -7 m 2 /s and 10 -2 m 2 /s (Song and Renner, 2006,0;[START_REF] Rempe | Influence of effective stress and pore fluid pressure on fault strength and slip localization in carbonate slip zones[END_REF][START_REF] Kranz | Hydraulic diffusivity measurements on laboratory rock samples using an oscillating pore pressure method[END_REF][START_REF] Wibberley | Hydraulic diffusivity of fault gouge zones and implications for thermal pressurization during seismic slip[END_REF]. In contrast, hydraulic tests in the field are inherently dominated by discrete fracture conduits. This results in an order of magnitude for the diffusivity that spans from 10 -1 m 2 /s to 10 1 m 2 /s [START_REF] Renner | Periodic pumping tests[END_REF][START_REF] Cheng | Exploratory use of periodic pumping tests for hydraulic characterization of faults[END_REF][START_REF] Maineult | Self-potential signals induced by periodic pumping tests[END_REF][START_REF] Talwani | In situ measurements of hydraulic properties of a shear zone in northwestern south carolina[END_REF][START_REF] Becker | Cross-hole periodic hydraulic testing of inter-well connectivity[END_REF][START_REF] Sayler | Understanding the geometry of connected fracture flow with multiperiod oscillatory hydraulic tests[END_REF]. In some cases, the observed hydraulic diffusivity might be much lower [START_REF] Doan | In situ measurement of the hydraulic diffusivity of the active chelungpu fault, taiwan[END_REF][START_REF] Xue | Continuous permeability measurements record healing inside the wenchuan earthquake fault zone[END_REF] or higher [START_REF] Becker | Cross-hole periodic hydraulic testing of inter-well connectivity[END_REF][START_REF] Guiltinan | Measuring well hydraulic connectivity in fractured bedrock using periodic slug tests[END_REF][START_REF] Sayler | Understanding the geometry of connected fracture flow with multiperiod oscillatory hydraulic tests[END_REF] depending on the site geology. Moreover, different test methods may also provide different diffusivity values, e.g., lower hydraulic diffusivities were observed in constant rate tests than in periodic tests [START_REF] Guiltinan | Measuring well hydraulic connectivity in fractured bedrock using periodic slug tests[END_REF]. At geothermal sites, the hydraulic diffusivity obtained by fitting seismic events also commonly varies between 10 -1 m 2 /s and 10 1 m 2 /s [START_REF] Shapiro | Estimating the crust permeability from fluid-injectioninduced seismic emission at the ktb site[END_REF][START_REF] Shapiro | Fluid-induced seismicity: Pressure diffusion and hydraulic fracturing[END_REF]. However, it is worth noting that these values represent the averaged hydraulic diffusivity of the whole tested fractured rock system. They are the combination of the matrix diffusivity and the fracture diffusivity [START_REF] Ortiz | Two-dimensional numerical investigations on the termination of bilinear flow in fractures[END_REF][START_REF] Sayler | Understanding the geometry of connected fracture flow with multiperiod oscillatory hydraulic tests[END_REF]. Hence, it is not surprising that our single fault model renders higher values of diffusivity. These results are also supported by previous studies (e.g., [START_REF] Sayler | Understanding the geometry of connected fracture flow with multiperiod oscillatory hydraulic tests[END_REF]), which evidenced that flow between an interval with large diffusivities (up to 10 3 m 2 /s) might be dominated by a constrained planar fracture.

For real-world case applications, other factors may also alter the hydraulic diffusivity, such as flow exchange between the fracture and matrix, mineral sealing and the temperature [START_REF] Wibberley | Hydraulic diffusivity of fault gouge zones and implications for thermal pressurization during seismic slip[END_REF]. Furthermore, hydraulic diffusivity has also been correlated with resolved fault movement. For instance, Guglielmi et al. (2015b) reported a wide range of diffusivities (10 -9 m 2 /s and 10 3 m 2 /s) during injection-induced fault reactivation experiments. These observations require further studies. In this work, we focused on understanding the impact of the fracture geometry on diffusivity, which is a fundamental topic. As such, one implication of our results is to provide a reference for complex numerical models. As an example, on the basis of the linear diffusion equation, [START_REF] Haagenson | Seismic diffusivity and the influence of heterogeneity on injection-induced seismicity[END_REF] used 2.2 × 10 2 -3.3 × 10 3 m 2 /s for the hydraulic diffusivity of each single fracture (compatible with our results) as an input in their 3D discrete fracture network and matrix (DFNM) numerical model and obtained an effective hydraulic diffusivity of 0.29 m 2 /s for the whole system (a common value in the field). We infer that the result might be improved if considering roughness and anisotropy (e.g., varied aperture distribution) for the input single fracture diffusivity.

Conclusion

We studied the effect of fracture surface roughness and fracture closure on pressure diffusion by numerically simulating transient fluid flow. The effect was evaluated quantitatively in terms of the effective hydraulic diffusivity D e . We considered the self-affinity property for the fracture surfaces as well as the fracture aperture. The implemented fracture geometry was based on synthetically generated surfaces/apertures following field observations. We performed transient pressure diffusion modeling in surface-represented rough fractures for different stages of fracture closure and ob-served that the roughness could significantly affect the effective hydraulic diffusivity of the fracture. At large openings, the rough fracture exhibits hydraulic behavior similar to the parallel plate model. As the fracture is gradually closed, the effective hydraulic diffusivity increasingly deviates from the parallel plate model and shows anisotropic behavior by enhancing or reducing the diffusivity according to the orientation of the pressure drop. Furthermore, when it approaches the percolation threshold, the increase in the fracture contact area and tortuous flow channels strongly decreases the effective hydraulic diffusivity by 7 and 4 orders of magnitude in the x-and y-directions, respectively. However, owing to the self-affinity property, a large residual opening (large diffusivity) exists even with a small hydraulic aperture. Although the method is based on a simple linear diffusion equation, our results show good consistency with some previously obtained field observations. Therefore, this study could have important implications for understanding the measurement of hydraulic properties as well as the associated fluid-induced seismicity pattern. The influence of the rock matrix and elastic fracture closure (for the volume representation) will be considered in future studies.

Appendix A

For 3D fracture flow, the fluid velocity distribution has the z component. To obtain the average velocity, the local velocity profile v x (z) in the z-direction is first assumed to follow a parabolic equation (similar to the velocity profile of the parallel plate model):

v x (x, y, z) = az 2 + bz + c.
(2.23)

The average velocity is defined as the average of the integral of the velocity function along the z-direction:

v x (x, y, z) = 1 h h 2 -h 2 v x (x, y, z) dz.
(2.24)

Appendix B

The procedure of searching for the effective hydraulic diffusivity of a rough fracture is as follows:

1) First, from the simulation results of a rough fracture flow, we obtain a time-dependent pressure distribution p R (x, y, t), which is averaged along the y-axis as < p R (x, t) > y . Fig. (2.21) shows the pressure distribution of the parallel plate model (d m = 0.3 m) at x = L/2 (averaged along the y-axis) plotted as a function of normalized time t/t * . The pressure is stabilized at t s ≈ 0.5t * (i.e., the fluid velocity is zero). Accordingly, we use t s as the upper bound of time to estimate the effective hydraulic diffusivity (i.e., from t = 0 to t = t s ). The time resolution matches the spatial resolution, i.e., we have 512 × 512 points on the fracture plane, so the number of time step is also 512. This is tested by a finer meshing resolution (1024 × 1024) that the results have tiny variations (less than an error of 0.1%), although the error is also very small with a more sparse meshing size 256 × 256 (< 1%). 2) The pressure distribution p ∥ (x, t, D f ) is calculated using different hydraulic diffusivities (hereinafter referred to as the test diffusivities) from the analytical expression Eq. (2.19). These test diffusivities are given in units of D m , which is the diffusivity for a parallel plate model with the same mean aperture d m as the rough fracture. For example, when we expect that the pressure diffusion of the rough fracture is reduced compared to the equivalent parallel plate model, we can build the set of test diffusivities in the range [0.1 : 0.9]D m with a step ∆D = 0.1D m . Alternatively, if it is expected to be enhanced, then the test diffusivity set can be [1 : 0.1 : 2]D m .

3) The differences between < p R (x, t) > y and p ∥ (x, t, D f ) in the least square sense is computed and the differences summed, followed by calculating the least square error for each test diffusivity.

The diffusivity with the minimum error is assigned as the diffusivity for the rough fracture. The best diffusivity is referred to as D e1 .

4) A set of diffusivities near D e1 with a smaller step such as ∆D = 0.01D m is rebuilt, and step 3 is repeated to obtain a more accurate diffusivity D e2 . The step is reduced to ∆D = 0.001D m and D e3 obtained. In theory, the smaller the step of the test diffusivity used, the higher the accuracy of the diffusivity. In this study, we calculate the test diffusivity up to three decimal places. The final best D e3 is regarded as the effective hydraulic diffusivity D e of the rough fracture. Note that for the obtained best diffusivity at each step, if it is at the boundary of the set of the test diffusivities, we rebuild the set by including this diffusivity inside it and repeat step 3 until an optimal hydraulic diffusivity inside this set is found.

5) The effect of different time resolutions is tested. When the resolution is halved to 256 time steps, the results remain almost unaffected. Even when 100 time steps are used, the results only show an error of less than 1%. Therefore, the results obtained at the time resolution we use are considered robust. L'écoulement de Navier-Stokes est résolu dans l'ouverture rugueuse tridimensionnelle de la fracture, et l'écoulement de Darcy est résolu dans la matrice poroélastique environnante. Nous étudions l'évolution de la rigidité mécanique et de la perméabilité aux fluides du système fracture-roche pendant la fermeture de la fracture en considérant le rendement de l'aspérité et le dépôt du matériau de remplissage de la fracture dans l'espace ouvert de la fracture rugueuse. L'évolution du volume de la faille, de la rigidité normale de la fracture et de la perméabilité est suivie jusqu'à ce que les seuils de percolation des fluides soient dépassés dans deux directions orthogonales du gradient de pression imposé. Enfin, nous proposons une loi basée sur la physique pour l'évolution de la rigidité et de la perméabilité en fonction du volume de la faille. Il est démontré que pendant la fermeture, la rigidité augmente exponentiellement lorsque le volume de la faille diminue. Une forte anisotropie de la perméabilité de la fracture est également mise en évidence lorsque les seuils de percolation sont atteints. Nous discutons de l'influence de ces résultats sur l'utilisation de la stimulation chimique dans les réservoirs fracturés.

Abstract

Natural or artificial fluid flow in deep fractured reservoirs, such as Enhanced Geothermal Systems (EGS), is primarily controlled by open fractures and faults, and is considered a key element for hydraulic performance. Flow along these fractures is strongly affected by channeling between fracture asperities and by deposits sealing the open fracture space due to mineral precipitation. Fracture asperities and fracture sealing also impact the mechanical behavior of fractures, especially their mechanical stiffness. Here, we study both the permeability and the stiffness of a rough fracture at the field scale during its closure. We base our approach on a well-established self-affine geometrical model for fracture roughness. We develop a finite element model based on the MOOSE/GOLEM framework and conduct numerical flow experiments in a 256 × 256 × 256 m 3 granite reservoir hosting a single, partially sealed fracture under variable normal loading conditions. Navier-Stokes flow is solved in the embedded 3-dimensional rough aperture, and Darcy flow is solved in the surrounding poroelastic matrix. We study the evolution of the mechanical stiffness and fluid permeability of the fracture-rock system during fracture closure by considering the asperity yield and the depositing of fracture-filling material in the open space of the rough fracture. The evolution of the fault volume, fracture's normal stiffness and permeability are monitored until fluid percolation thresholds are exceeded in two orthogonal directions of the imposed pressure gradient. Finally, we propose a physically based law for the stiffness and permeability evolution as a function of the fault volume. It is demonstrated that during closure, stiffness increases exponentially as the fault volume decreases. A strong anisotropy of the fracture permeability is also evidenced when reaching percolation thresholds. We discuss how the results influence the use of chemical stimulation in fractured reservoirs.

Introduction

Fractures play a major role as preferential fluid pathways for numerous georesource applications, such as deep geothermal reservoirs [START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF][START_REF] Cox | Coupling between deformation, fluid pressures, and fluid flow in oreproducing hydrothermal systems at depth in the crust[END_REF][START_REF] Gudmundsson | Rock fractures in geological processes[END_REF][START_REF] Huenges | Geothermal energy systems: exploration, development, and utilization[END_REF]. Indeed, they provide major conduits at all scales through which fluids can flow [START_REF] Caine | Fault zone architecture and permeability structure[END_REF][START_REF] Berkowitz | Characterizing flow and transport in fractured geological media: A review[END_REF]. They are also critical for the mechanical behavior of fractured rock masses owing to their ubiquity and the fact that they are typically mechanically weak objects (Zimmerman and Main, 2004). Subsequently, the hydromechanical behavior of rock fractures is of central importance for determining the transport and deformation properties of faulted reservoirs. This behavior can be studied on the scale of a fractured rock mass that includes fracture networks with many fractures but also at the scale of a single fracture [START_REF] Adler | Fractures and fracture networks[END_REF][START_REF] De Dreuzy | Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks (dfn)[END_REF][START_REF] Cornet | Elements of crustal geomechanics[END_REF][START_REF] Scholz | The mechanics of earthquakes and faulting[END_REF][START_REF] Kluge | Hydraulic-mechanical properties of microfaults in granitic rock using the punch-through shear test[END_REF][START_REF] Deng | Hydraulic diffusivity of a partially open rough fracture[END_REF]. The behavior of a single fracture must be addressed first to understand that of the whole fractured mass.

As demonstrated by numerous studies, the fundamental properties of individual fractures strongly rely on the geometry of the fracture void space and the contact area between the facing fracture walls [START_REF] Barton | The shear strength of rock joints in theory and practice[END_REF][START_REF] Jaeger | Effect of fluid pressure heterogeneity on injection-induced fracture activation[END_REF][START_REF] Gudmundsson | Rock fractures in geological processes[END_REF][START_REF] Scholz | The mechanics of earthquakes and faulting[END_REF]Kluge et al., 2021a). Fracture roughness, which controls the fracture's void space geometry, has been extensively analyzed in the past by using various techniques to measure the fracture surface topography [START_REF] Tse | Estimating joint roughness coefficients[END_REF][START_REF] Brown | Broad bandwidth study of the topography of natural rock surfaces[END_REF][START_REF] Power | Euclidean and fractal models for the description of rock surface roughness[END_REF]Schmittbuhl et al., 1995a;[START_REF] Bouchaud | Scaling properties of cracks[END_REF]. Among the techniques that characterize the fracture and fault roughness from topography measurements (Schmittbuhl et al., 1995b;[START_REF] Candela | Characterization of fault roughness at various scales: Implications of three-dimensional high resolution topography measurements[END_REF], spectral analysis is more generally used to demonstrate the power law behavior of the height-height correlation function [START_REF] Jacobs | Quantitative characterization of surface topography using spectral analysis[END_REF]. It is shown that the fracture surfaces of different rock materials, with different fracture modes and within different geological contexts share a self-affine scaling invariance that extends up to a multi-kilometric fault scale with no roll-off wavelength at large scales (Schmittbuhl et al., 1993a;[START_REF] Fardin | The scale dependence of rock joint surface roughness[END_REF][START_REF] Sagy | Evolution of fault-surface roughness with slip[END_REF][START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. In some cases, the self-affine geometrical model of the fault surface exhibits anisotropy related to the main slip direction with anisotropy in the Hurst or roughness exponent [START_REF] Renard | High resolution 3d laser scanner measurements of a strike-slip fault quantify its morphological anisotropy at all scales[END_REF][START_REF] Brodsky | Faults smooth gradually as a function of slip[END_REF][START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF].

Single-phase fluid flow in the open void space of rough fractures has been extensively studied, particularly in the framework of the lubrication approximation, where the roughness is considered as a perturbation to a parallel plate model [START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF]. The cubic law is then applied to relate the mechanical aperture of the fracture, i.e., the average of the geometrical aperture, to the fracture permeability. In addition, based on its hydraulic properties, the hydraulic aperture of a fracture is defined. Numerous models and experiments have been developed to compare the mechanical and hydraulic apertures of fractures for different flow regimes [START_REF] Mourzenko | Permeability of a single fracture; validity of the reynolds equation[END_REF][START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF][START_REF] Méheust | Flow enhancement of a rough fracture[END_REF][START_REF] Matsuki | Size effect on aperture and permeability of a fracture as estimated in large synthetic fractures[END_REF][START_REF] Méheust | Geometrical heterogeneities and permeability anisotropy of rough fractures[END_REF][START_REF] Ishibashi | Beyond-laboratory-scale prediction for channeling flows through subsurface rock fractures with heterogeneous aperture distributions revealed by laboratory evaluation[END_REF]. These works demonstrate the impact of the fracture roughness on the fracture permeability with regard to the enhancement or inhibition of the flow owing to strong channeling of the flow for a specific orientation of the pressure drop. Moreover, the largest wavelengths in the aperture field are shown to dominate the influence of the roughness (Neuville et al., 2010a,0;Wang et al., 2016). In the case of many independent fractures, fracture roughness is shown, on average, to inhibit flow when a hydraulic aperture is smaller than the mechanical aperture as the fracture surfaces come into contact. However, significant departures from this average behavior exist since the behavior is dominated by the largest wavelengths, which are not stationary in space and are accordingly specific to each surface.

As a fracture is mechanically loaded or chemically sealed, its contact area is expected to increase, and a percolation threshold for fluid flow might be reached [START_REF] Sahimi | Flow phenomena in rocks: from continuum models to fractals, percolation, cellular automata, and simulated annealing[END_REF]Schmittbuhl et al., 1993b;[START_REF] Dapp | Self-affine elastic contacts: percolation and leakage[END_REF]. The contact mechanics of self-affine surfaces have been explored by numerous studies [START_REF] Borri-Brunetto | Scaling phenomena due to fractal contact in concrete and rock fractures[END_REF][START_REF] Batrouni | Elastic response of rough surfaces in partial contact[END_REF][START_REF] Persson | On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion[END_REF][START_REF] Hyun | Finite-element analysis of contact between elastic self-affine surfaces[END_REF][START_REF] Pei | Finite element modeling of elasto-plastic contact between rough surfaces[END_REF][START_REF] Persson | Contact mechanics for randomly rough surfaces[END_REF][START_REF] Carbone | Asperity contact theories: Do they predict linearity between contact area and load[END_REF][START_REF] Vallet | Sampling effect on contact and transport properties between fractal surfaces[END_REF][START_REF] Almqvist | Interfacial separation between elastic solids with randomly rough surfaces: comparison between theory and numerical techniques[END_REF][START_REF] Pastewka | Finitesize scaling in the interfacial stiffness of rough elastic contacts[END_REF]. They show specific properties such as 1) a linear relationship between the contact area and normal loading; 2) an exponential evolution of the normal stress with respect to the normal displacement; and 3) a linear increase in the normal fracture stiffness with normal stress. These studies have also shown that the normal stress field along asperities exhibits self-affine properties [START_REF] Hansen | Normal stress distribution of rough surfaces in contact[END_REF][START_REF] Schmittbuhl | Are stress distributions along faults the signature of asperity squeeze[END_REF].

Natural sealing of fractures derives from many processes and is dependent on the fluid chemistry, temperature and pressure [START_REF] Parry | Fault-fluid compositions from fluid-inclusion observations and solubilities of fracture-sealing minerals[END_REF]. Following the pioneering work of [START_REF] Durney | Incremental strains measured by syntectonic crystal growths[END_REF], the following three fracture filling modes are defined: the syntaxial mode for inward growth, antitaxial mode for outward growth, and stretching mode for an unspecific growth direction. The growth of seals in deforming rocks can arise from a single event, but it is often the result of multiple crack sealing events following fracture growth [START_REF] Bons | A review of the formation of tectonic veins and their microstructures[END_REF]. Then, the fibrous habits of crystalline materials induced by pressure solution in the rock matrix are shown [START_REF] Ramsay | The crack-seal mechanism of rock deformation[END_REF][START_REF] Laubach | Practical approaches to identifying sealed and open fractures[END_REF]. Dissolution and precipitation along fracture walls owing to secondary mineral de-position or artificial reactive fluid injection might also significantly modify the permeability and stiffness of fractures [START_REF] Laubach | Practical approaches to identifying sealed and open fractures[END_REF][START_REF] Liu | A fully-coupled hydrological-mechanical-chemical model for fracture sealing and preferential opening[END_REF][START_REF] Elkhoury | Dissolution and deformation in fractured carbonates caused by flow of co2-rich brine under reservoir conditions[END_REF][START_REF] Noiriel | 3d analysis of geometry and flow changes in a limestone fracture during dissolution[END_REF][START_REF] Lang | Hydraulic sealing due to pressure solution contact zone growth in siliciclastic rock fractures[END_REF]. Changes in stiffness from sealing bridge development along fractures can have large-scale implications by inducing, for instance, seismic velocity variations [START_REF] Sayers | The effect of mineralization on the ratio of normal to tangential compliance of fractures[END_REF].

Our study aims to develop a numerical model to address a specific aspect of the mechanical closure of a fracture related to sealing processes, i.e., the coupling between the normal stiffness and the permeability of the fracture, and we propose that this model be applied at the field scale using the multiscale knowledge of fracture geometry. For the description of the sealing process, we consider syntaxial mineral growth, i.e., from the walls of the fracture towards its center in multiple separate steps. We therefore assume nucleation of new crystals to be negligible compared to the extension of existing crystals. Minerals are supposed to be typically elongated in the growth direction perpendicular to the fracture walls, forming homogeneous layers in which the thickness is controlled by precipitation kinetics and the fracture opening rate [START_REF] Ramsay | The crack-seal mechanism of rock deformation[END_REF][START_REF] Bons | A review of the formation of tectonic veins and their microstructures[END_REF]. Successive layers then fill the open space of the fracture and modify the permeability and stiffness of the fracture. Their effects on the fracture stiffness result from the competition between the progressive evolution of the contact area owing to the mechanical closure of the fracture filling with competent cement and the generation of new bridges between the fracture walls in the asperities neighbourhood. One specific aspect of our approach is the consideration of the multiscale realistic geometry of the fracture, i.e., a self-affine geometry [START_REF] Renard | High resolution 3d laser scanner measurements of a strike-slip fault quantify its morphological anisotropy at all scales[END_REF][START_REF] Brodsky | Faults smooth gradually as a function of slip[END_REF][START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF], and the analysis of the impact of the fracture closure at a field scale of 256 m [START_REF] Blöcher | Hydraulic history and current state of the deep geothermal reservoir groß schönebeck[END_REF]. We simplify the details of the sealing mineral deposition to a progressive filling of the self-affine aperture by a set of discrete homogeneous layers. We analyze the impact of the healing rate (i.e., the amount of healing material in the open fracture) on the flow channeling related to the fracture roughness and on the stiffness of the fracture when sealing is reaching the percolation threshold for fluid flow. Finally, we discuss the hydraulic and geomechanical implications of the sealing process for deep EGS geothermal reservoirs where chemical stimulation is considered to modify the fracture healing material [START_REF] Nami | Chemical stimulation operations for reservoir development of the deep crystalline hdr/egs system at soultz-sous-forêts (france)[END_REF][START_REF] Portier | Chemical stimulation techniques for geothermal wells: experiments on the three-well egs system at soultz-sous-forêts, france[END_REF][START_REF] Huenges | Geothermal energy systems: exploration, development, and utilization[END_REF].

Model description

A partially open fracture embedded in a porous rock under normal load

We model the fluid circulation with a large-scale partially open fracture embedded in a 3D poroelastic block (Fig. 3.1). The dimensions of the entire block are 256 × 256 × 256 m 3 , which represent a characteristic scale of a typical fractured/faulted reservoir, e.g., the Soultz-sous-Forêts deep geothermal reservoir [START_REF] Vallier | Large scale hydro-thermal circulation in the deep geothermal reservoir of soultz-sous-forêts (france)[END_REF]. The fracture is partially open due to the presence of fracture asperities that define open domains where fluid can flow. We assume drained conditions, i.e., fluid is able to move along the fracture due to an imposed pressure difference at the boundaries. Fluid circulation is modeled by imposing a pressure gradient (constant pressures at the inlet and outlet boundaries) along either the x-direction or y-direction of the block to investigate a possible anisotropy in the fluid flow behavior. The contact areas between asperities are the domains where stress is transferred with no fluid flow, i.e., closed domains. We assume that the fracture is mechanically loaded by a normal stress σ n applied along the top of the block. Because of the partial contact along the fracture, the applied stress on the whole block boundary is amplified along the contact zones following Amonton's law, as shown in Fig. 3.1a [START_REF] Scholz | The mechanics of earthquakes and faulting[END_REF]. Fig. 3.1b shows an example of the fluid flow both in the pores of the matrix and in the open space of the fracture for a case when the pressure gradient is imposed along the x-direction.

Fracture aperture

Faults and fractures are typically considered two-dimensional objects. This is because their thicknesses are often considered negligible with respect to their extension. They are then approximated as planar surfaces. However, even as 2D objects, they show multiscale corrugations that can be described by a self-affine topography over a very large range of scales (up to 9 decades of length scales) [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. A 2D cross-section of the fault will then remain statistically invariant under any of the following scaling transformation: δx → λδx, δz → λ H δz whatever the value of the coefficient λ, where δx is the coordinate along the section, δz is the topography amplitude and H is the Hurst or roughness exponent (Schmittbuhl et al., 1995a;[START_REF] Meakin | Fractals, scaling and growth far from equilibrium[END_REF]. Moreover, [START_REF] Renard | High resolution 3d laser scanner measurements of a strike-slip fault quantify its morphological anisotropy at all scales[END_REF] and [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF] showed that faults/fractures exhibit a certain anisotropy in their roughness exponent in response to the fault slip history, resulting in a roughness exponent on the order of H ∥ = 0.6 along the slip direction and H ⊥ = 0.8 in the slip perpendicular direction.

Here, we consider a 256 m × 256 m fault that reproduces the geometry of natural faults over the maximum range of scales in terms of scaling of the power spectral density (PSD) and the root mean square (RMS) of the fracture geometry (see supplementary material, Fig. 3.18) [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. The spectral properties of the synthetic fault geometry are consistent with those of a natural fault. We assume that the time scale of fluid flow is much shorter than that of the long-term shear slip along the fault responsible for the anisotropy in the fault geometry. The fault geometry is then considered fixed during the flow simulation. Based on the fault topography, we build the fault aperture with the following assumption regarding the fault slip history: the cumulative fault shear offset is supposed to be large compared to the size of the fault area we are considering [START_REF] Plouraboué | Aperture of rough cracks[END_REF]. In this case, correlations between the two facing surfaces are lost, and the fault aperture is composed of the facing of two independent fault surfaces and subsequently has the same scaling properties of each fracture surface.

Accordingly, to reproduce field observations, we generate the fault aperture as an anisotropic self-affine surface with H ∥ = 0.6 along the historical slip direction (y-axis) of the fault, H ⊥ = 0.8 along the perpendicular direction (x-axis) and an RMS amplitude of 0.09 m for the fault topography at the scale of 256 m (see Fig. 3.2a) [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. Fig. 3.2b shows the histogram of the aperture distribution in the entire fault domain of 256 m × 256 m with no normal deformation of the asperities or no filling material in the openings of the fault. The measured distribution is well approximated by a Gaussian distribution. It is of interest to note that such a fault model predicts a typical magnitude for fault asperities on the order of 0.09 m for a fault with an extension of 256 m. The aspect ratio of such a 3D fault geometry is then on the order of 0.09 m/256 m ≈ 3 • 10 -4 . We also see that the maximal amplitude of the asperities is on the order of 3h rms ≈ 0.27 m. In the (unrealistic) case of a fully rigid matrix, i.e., no deformation of the largest asperities, this would correspond to the typical aperture of the fracture. As discussed below, the aperture will actually decrease when considering elasto-plastic deformation of the asperities under normal mechanical loading of the fault. The aperture is also expected to decrease when healing processes develop over time.

To explore the variability of the fracture geometry, the following six independent synthetic faults are generated: F1, F2, F3, F4, F5 and F6. They all share the same self-affine scaling properties, i.e., the same roughness exponents and roughness amplitude at the same scale. The only difference among them is that different seeds are selected for generating the initial white noise signal, which is used to construct the self-affine fault surface [START_REF] Candela | Characterization of fault roughness at various scales: Implications of three-dimensional high resolution topography measurements[END_REF]. As such, the six faults have different surface or aperture geometries but still share the same statistical properties regarding their distribution. The aperture map of F1 prior to any mechanical loading or sealing is shown in Fig. 3.2a. The maps for the other faults are given in Fig. 3.19.

A kinematic model for reconstructing the fracture closure

In this study, we assess the impact of sealing on the mechanical and hydraulic properties of a fracture by modeling its closure following a specific conceptual route. We reconstruct the closure of the fracture assuming that it is the result of two different mechanisms (Fig. 3.3). First, we consider that the fracture is only mechanically closed due to an increase in normal displacement applied along the top of block u b from normal loading σ n (see Fig. 3.4). No sealing process is considered during this phase. The matrix and the fracture are shown as strings in series for which the bulk displacement u b is the sum of the displacement of the matrix u m and that of the fracture u f . In other words, the fracture displacement u f is equivalent to the difference between the vertical displacement of the fractured rock u b and the intact rock matrix u m , as follows:

u f = u b -u m (3.1)
Note that in this paper, the direction along negative z is considered positive. Therefore, the fracture closure u f is positive with respect to the zero displacement of the first contact.

To simplify the contact rheology, we assume that the deformation of asperities follows a rigidplastic behavior (i.e., rigid solid with a low yield stress). The normal loading then creates contact areas where stresses strongly concentrate and rapidly overpass the elastic limit, leading to a perfectly plastic behavior [START_REF] Pei | Finite element modeling of elasto-plastic contact between rough surfaces[END_REF]. This is similar to an 'overlap' model where overlapping asperities are erased with volume losses (Fig. 3.3c) [START_REF] Méheust | Scale effects related to flow in rough fractures[END_REF][START_REF] Kim | Fracture flow simulation using a finite-difference lattice boltzmann method[END_REF][START_REF] Walsh | Numerical modeling of stress-permeability coupling in rough fractures[END_REF]. This approach is a proxy for a full elastoplastic model. This minimizes contact areas since contact regions develop mainly along the highest asperities in compacted clustered regions and the deformation is local without volume conservation. In contrast, in a fully elastic model, contact areas would spread out over multiple regions because of non-local deformation and the additional requirement of volume preservation. Our choice of an 'overlap' model stems from this model being supported by laboratory observations where deformation in contact/overlap areas produces gouge particles that are flowing out of the contact areas (Kluge et al., 2021b). The initial mechanical load leads to rigid plastic deformation and accumulation of elastic strain, and a stress release would lead to an elastic relaxation as indicated by the dashed red lines on the x-axis (Kluge et al., 2021b); in this case, plastic deformation will remain. We consider mechanical closure to be part of the initial history of the fracture, and it occurs before sealing begins to be effective. We assume that this phase leads to an overall contact area on the order of 20%, as observed in fracture closure experiments for a normal load on the order of several tens of MPa (Kluge et al., 2021b). During this mechanical phase of closure, the normal stress along the contact areas increases as the fracture displacement rises, leading to a more significant closure of the fracture, as shown in Fig. 3.5.

The second mechanism for fracture closure is supposed to be due to a sealing process that occurs in the later stage, upon termination of the mechanical phase (i.e., typically when the contact area > 20%) at constant imposed normal displacement and normal stress (Fig. 3.3c and Fig. 3.5). Fracture sealing is considered to be the sole result of mineral precipitation of successive layers, as commonly observed in fractured reservoirs (e.g., [START_REF] Griffiths | Geothermal implications for fracture-filling hydrothermal precipitation[END_REF]). Homogeneous layers of secondary minerals are observed to be deposited with a rather constant thickness, following a syntaxial growth of the filling material, both in the field [START_REF] Durney | Incremental strains measured by syntectonic crystal growths[END_REF][START_REF] Davis | Energy of injectioninduced seismicity predicted from in-situ experiments[END_REF] and in laboratory hydrothermal experiments [START_REF] Okamoto | Textures of syntaxial quartz veins synthesized by hydrothermal experiments[END_REF]. Specifically, during syntaxial growth, filling materials are progressively added from the fracture walls to the middle (inward) as opposed to during antitaxial growth, where the filling materials are added from the middle to the walls (outward) [START_REF] Davis | Energy of injectioninduced seismicity predicted from in-situ experiments[END_REF]. Sealing would increase the stiffness of the system, as indicated by the extended dashed line in Fig. 3.5. Therefore, the system could sustain more stress without additional plastic deformation. This will also lead to a transformation into a more elastic regime.

It is then of interest to note that this second part of the closure modeling can also be viewed as an 'overlap' model since asperities are also erased due to a different physical mechanism, namely, during chemical sealing, the contact areas are expanding because of the outer growth of the asperity boundaries, while during mechanical closure, plastic deformation induces the inner movement of the asperity boundaries. Subsequently, we expect that in both cases, the stiffness and permeability of the fracture will be impacted in a similar manner and, thus, that they can be described by a common generic law. The whole fracture closure process can then be viewed as a progressive, continuous translation of one surface towards the other from the first contact points when the closure is zero to the full contact when the closure is at its maximum. This corresponds to a kinetic description of the aperture geometry evolution caused by steps of imposed displacements of one facing fracture penetrating the other. During the mechanical phase, the closure is equal to the relative imposed displacement u f from its value at the first contact. During the sealing phase, the closure steps correspond to the thickness of the deposited layer. By gradually closing a fracture using this kinematic model, the contact area y = A c /A is shown to follow an error function law, i.e., y = (1 + erf (x/ √ 2))/2, as a function of the normalized closure x = (u f -h)/h rms , where u f is the closure, h is the mechanical aperture of the initial fracture and h rms its RMS. This law describes the cumulative distribution function (CDF) of a normal distribution; see the orange line in Fig. 3.6. The pink area shows the mechanical closure (below 20% of the contact area), while the light blue area shows the sealing closure (above 20% of contact area).

Numerical method

Analyzing the stress, pressure and velocity fields in the fracture-matrix system requires an adequate physical model of the processes involved and their coupling. Here, the Multiphysics Object Oriented Simulation Environment (MOOSE) in combination with GOLEM, a MOOSE-based application for modeling coupled thermal-hydraulic-mechanical (THM) processes in fractured and faulted geothermal reservoirs, is used [START_REF] Jacquey | Modelling coupled fluid flow and heat transfer in fractured reservoirs: description of a 3d benchmark numerical case[END_REF][START_REF] Peterson | Overview of the incompressible navier-stokes simulation capabilities in the moose framework[END_REF]. The rock matrix is considered a poroelastic medium, formulated in a Lagrangian framework, as implemented in [START_REF] Cacace | Flexible parallel implicit modelling of coupled thermalhydraulic-mechanical processes in fractured rocks[END_REF]. For the solid phase in the matrix, the mass conservation reads as follows:

(1 -n) ρ s D s ρ s Dt - D s n Dt + (1 -n)∇ • v s = 0 (3.2)
where n is the porosity, ρ s is the density of the solid, D s /Dt is the Lagrangian derivative with respect to the moving solid and v s is the solid velocity. The moment balance equation leads to the following for solving the displacement vector in terms of the effective stress σ ′ :

∇ • σ ′ -αp f 1 + ρ b g = 0 (3.3)
where α, g, 1 is the Biot coefficient, gravitational acceleration and the rank-two identity tensor.

ρ b = nρ f + (1 -n)ρ s
is the bulk density of the fluid-solid mixture.

For the fluid phase in the matrix, we assume a Darcy's flow, as follows:

q D = n(v f -v s ) = - k m µ f • (∇p f -ρ f g) (3.4)
where q D is the volumetric flow rate per unit surface area, k m is the permeability tensor of the rock matrix, and µ f is the dynamic fluid viscosity. The fluid mass conservation then renders the following:

n ρ f D f ρ f Dt + D s n Dt + n∇ • v s + ∇ • q D = 0 (3.5)
where ρ f is the density of the fluid and D f /Dt is the Lagrangian derivative with respect to the moving fluid.

For the fluid phase in the open space of the fracture, the incompressible Navier-Stokes (INS) equation is implemented as follows [START_REF] Peterson | Overview of the incompressible navier-stokes simulation capabilities in the moose framework[END_REF]:

ρ f D f v f Dt + ∇p f -µ f ∇ 2 v f -ρ f g = 0 (3.6)
The system of equations is spatially discretized by the finite-element method, while temporal discretization relies on a finite difference implicit Euler scheme [START_REF] Cacace | Flexible parallel implicit modelling of coupled thermalhydraulic-mechanical processes in fractured rocks[END_REF]. The Navier-Stokes module provides the solution for the pressure and the velocity field inside the fracture only [START_REF] Blöcher | Permeability of matrix-fracture systems under mechanical loading-constraints from laboratory experiments and 3-d numerical modelling[END_REF]. For the rock matrix, the primary variables are the pore pressure and the solid displacements. The fluid velocity in the porous matrix and the strain and stress fields are derived from these primary quantities. To derive a coupling between the two physical systems, i.e., the open fracture and the rock domain, we define four exchange functions along the interfaces between the rock and the fracture (Fig. 3.7). Since the fluid pressure inside the fracture occurs almost instantaneously when compared to the time scale of pore pressure diffusion in the rock, we use this pressure as a boundary condition to solve for the pore pressure in the rock matrix system (Fig. 3.7a). Furthermore, this pressure acts as an additional body force on the rock matrix. Therefore, we use the pressure as a stress boundary condition in the rock matrix system (Fig. 3.7a). Based on this pore pressure, we calculate the fluid velocity in the rock domain following Darcy's law. The fluid velocity is then assigned as a slipping boundary for the fracture flow inside the fault (Fig. 3.7b). To impose the mechanical coupling and the corresponding deformation of the fault, the displacement values obtained in the rock matrix are assigned as first-order boundary conditions for the fault surfaces (Fig. 3.7b). Since the incompressible Navier-Stokes flow does not consider the solid displacement as a primary variable, we use the displacement vector field to geometrically alter the fault mesh at each time step, resulting in a displaced mesh configuration that is used to update the final solution for the dynamic pressure and fluid velocity field [START_REF] Blöcher | Permeability of matrix-fracture systems under mechanical loading-constraints from laboratory experiments and 3-d numerical modelling[END_REF].

To derive a solution for the fluid flow in the fracture-rock system and measure its permeability, we apply the following initial and boundary conditions (Fig. 3.7). The initial condition for the rock matrix is a zero displacement field and zero pore pressure. Furthermore, the fluid pressure and the velocity inside the fracture are also set to zero. In addition to the boundary conditions applied at the fault-matrix interface, we impose a flow either in the x-or y-direction by setting a small inlet overpressure of 10 Pa either along the x-or y-direction, respectively. This is to prevent the fluid velocity from being too high (see Fig. 3.12). At the outlet, we keep the pressure constant at zero during the total transient simulation. For the mechanical part, we set a zero displacement in the z-direction for the bottom and no displacement in the x-and y-directions for the lateral borders, as well as for the inlet and outlet. At the top of the block, we apply normal loading to assess the mechanical fracture closure. All material proprieties used in the modeling are listed in Tab. 3.1. These values are representative of the well-known EGS project Soultz-sous-Forêts [START_REF] Vallier | Large scale hydro-thermal circulation in the deep geothermal reservoir of soultz-sous-forêts (france)[END_REF].

Mesh generation

The open fracture space is used as a 3D domain for finite element mesh generation. The domain is based on a regular grid along the mean plane of the fracture with a grid spacing of 1 m and dimensions of 256 m by 256 m (Fig. 3.8a). The aperture ranges from a minimum of 0 m (i.e., contact areas) to a maximum of ∼ 0.5 m (i.e., at zero closure).

The magnitude of the fracture opening is on the order of 0.1 m (Fig. 3.2c), a characteristic value for fractures at the studied scale (∼ 256 m) [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF]. The challenge for the model is to account for the extreme aspect ratio of the fracture volume, which can be as high as 1:1000, making its meshing by isotropic voxels or elements particularly difficult for a similar resolution in the zdirection compared to the x-and y-directions. To represent a fracture opening with such an aspect ratio, we instead implement prism elements within an anisotropic mesh (Fig. 3.8b). Normal to the fracture mean plane (i.e., along the z-direction), we use a discretization of the fracture open space by 8 layers of prism elements, with an aspect ratio of 1:25. The minimum element sizes normal and parallel to the fracture are 0.0125 m and 1 m, respectively. For meshing the solid part, we use an isotropic tetrahedral mesh with an element size of ∼ 0.1 m close to the fracture, while the element size coarsens gradually with increasing distance from the fracture (Fig. 3.8a). Meshing is performed by Gmsh [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in preand post-processing facilities[END_REF], an open source finite element mesh generator. For the open fracture space, we generate between 71,152 and 262,080 prism elements according to the closure degree of the fracture. For the solid part, we generate between 1,563,350 and 2,620,330 tetrahedral elements.

Based on the mesh, we use the flood-fill algorithm [START_REF] Nosal | Flood-fill algorithms used for passive acoustic detection and tracking[END_REF] to determine the flow percolation threshold of a fault, i.e., the state where the fracture flow loses connection between the inflow and outflow borders. If no flow path is found, the fault is marked as having reached/exceeded the percolation threshold.

Fracture stiffness evolution during closure

The normal elastic response of a fracture is often quantified by the fracture's normal stiffness κ [START_REF] Goodman | A model for the mechanics of jointed rock[END_REF][START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF][START_REF] Morris | Numerical investigation of alternative fracture stiffness measures and their respective scaling behaviours[END_REF][START_REF] Jaeger | Effect of fluid pressure heterogeneity on injection-induced fracture activation[END_REF][START_REF] Lang | Hydraulic sealing due to pressure solution contact zone growth in siliciclastic rock fractures[END_REF]:

κ = dσ n du f (3.7)
where σ n is the applied normal stress and u f is the associated displacement of the fracture, considered a deformable body (Fig. 3.4). Eq. (3.7) depicts the rigidity of the fracture as stress changes due to the closure of the fracture and the induced changes in the contact areas along the fracture.

Stress-linear stiffness approach

Fracture stiffness measurements have been obtained using a variety of methods, such as laboratory experiments [START_REF] Sun | Rock joint compliance tests for compression and shear loads[END_REF][START_REF] Schrauf | Laboratory studies of gas flow through a single natural fracture[END_REF][START_REF] Bandis | Fundamentals of rock joint deformation[END_REF][START_REF] Cook | Natural joints in rock: mechanical, hydraulic and seismic behaviour and properties under normal stress[END_REF][START_REF] Zhao | Hydro-thermo-mechanical properties of joints in the carnmenellis granite[END_REF][START_REF] Kluge | Hydraulic-mechanical properties of microfaults in granitic rock using the punch-through shear test[END_REF], field tests [START_REF] Jung | Hydraulic in situ investigations of an artificial fracture in the falkenberg granite[END_REF][START_REF] Pratt | Elastic and transport properties of an in situ jointed granite[END_REF][START_REF] Rutqvist | Determination of fracture storativity in hard rocks using high-pressure injection testing[END_REF], and numerical simulations (Pyrak- [START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF][START_REF] Wang | Development of an empirical model relating permeability and specific stiffness for rough fractures from numerical deformation experiments[END_REF][START_REF] Morris | Numerical investigation of alternative fracture stiffness measures and their respective scaling behaviours[END_REF][START_REF] Lavrov | Fracture permeability under normal stress: a fully computational approach[END_REF]. Most studies measure a nonlinear increase in normal stress as a function of ongoing closure due to a linear increase in the stiffness from the applied stress, at least during the initial load of the fracture [START_REF] Cook | Natural joints in rock: mechanical, hydraulic and seismic behaviour and properties under normal stress[END_REF]. Such a process has also been shown to be irreversible when cyclic loading is applied [START_REF] Kluge | Hydraulic-mechanical properties of microfaults in granitic rock using the punch-through shear test[END_REF]. In fact, [START_REF] Berthoud | Shear stiffness of a solid-solid multicontact interface[END_REF] directly showed this linear behavior of the fracture stiffness, as follows:

κ = σ n u 0 (3.8)
where u 0 is a characteristic length related to the roughness of the surface. As proposed by [START_REF] Cook | Natural joints in rock: mechanical, hydraulic and seismic behaviour and properties under normal stress[END_REF], this linear property of the stiffness (Eq. 3.8) together with its definition (Eq. 3.7) lead to a simple differential equation [START_REF] Pastewka | Finitesize scaling in the interfacial stiffness of rough elastic contacts[END_REF], as follows:

dσ n du f = σ n u 0 (3.9)
which has a solution in the form of a semilogarithmic behavior of the stress, as classically observed in the following equation [START_REF] Zangerl | Normal stiffness of fractures in granitic rock: A compilation of laboratory and in-situ experiments[END_REF]:

ln (σ n ) = χ • u f + ln (σ ref n ) (3.10)
if one introduces the characteristic stiffness χ = 1/u 0 and σ ref n as the reference normal stress, which refers to the effective normal stress level at the beginning of the test when u b = 0, and is generally a small quantity [START_REF] Goodman | A model for the mechanics of jointed rock[END_REF]. The characteristic stiffness χ denotes the slope of the logstress against the closure curve and defines the rate of change of normal stiffness and normal stress [START_REF] Evans | The effects of fracture normal compliance on the long term circulation behavior of a hot dry rock reservoir: a parameter study using the new fully-coupled code'fracture[END_REF].

From Eq. (3.8) and Eq. (3.10), the behavior of the stiffness κ(u f ) can be predicted as follows:

ln (κ) = χ • u f + ln (χ • σ ref n ) (3.11)
In our model, to compute the fracture stiffness for a given closure u f , we apply a small stress perturbation (∆σ n = 1 MPa) along the top of the model and compute the corresponding volume change of the fracture ∆V f . Dividing the volume change ∆V f by the apparent fracture area A provides an estimate of the induced closure change ∆u f = ∆V f /A. The fracture stiffness κ at each closure step u f can then be estimated as κ = ∆σ n /∆u f . Eq. (3.11) is then used to fit the evolution of the stiffness during closure and to obtain the following missing parameters: χ and σ ref n . Fig. 3.9a shows the evolution of the computed stiffness during closure for fault F1 (see Fig. 3.20 for the other fault configurations). It demonstrates that Eq. (3.11) is a good framework to describe the stiffness evolution since the general behavior is well captured and validates the approximation made (in particular, the stress-linear stiffness approach -Eq. (3.8), see Fig. 3.9b). The stiffness characteristic χ for fault F1 is estimated as χ = 23.2 m -1 (= 0.0232 mm -1 ), which is a small value compared to the values (10 -70 mm -1 ) measured at the laboratory scale [START_REF] Zangerl | Normal stiffness of fractures in granitic rock: A compilation of laboratory and in-situ experiments[END_REF]. For the other five faults, we obtain similar stiffness characteristic values of 18.7 m -1 , 23.2 m -1 , 23.0 m -1 , 24.6 m -1 and 27.6 m -1 (see supplementary material,Fig. 3.20). Previous studies (e.g., [START_REF] Pastewka | Finitesize scaling in the interfacial stiffness of rough elastic contacts[END_REF]) have shown that the reciprocal of the stiffness characteristic 1/χ (the characteristic length u 0 ) is on the order of the RMS roughness of the surfaces for nonadhesive solid contact surfaces. In our study, we obtain stiffness characteristics in the range of 18 -27 m -1 for the 6 faults, indicating a characteristic length in the range of 0.04 -0.05 m, which is consistent with the RMS of the fault surface (i.e., 0.045 m).

A yield-contact area approach for closure

Eq. ( 3.9) provides a physical link between the characteristic stiffness χ and the geometry of the fracture. Once the stiffness and the characteristic stiffness are known, we obtain the normal stress by means of Eq. (3.8). Fig. 3.10 shows the normal stress evolution during the mechanical closure stage using the stress-linear stiffness approach for fault F1 (see Fig. 3.21 for other faults). For comparison, we estimate the evolution of the absolute normal stress σ n using an independent approach based on the assumption that contact areas act as yield regions. Because the ratio of the applied normal stress and the normal stress along the fault asperities (σ n /σ asperities n ) is equal to the ratio of the apparent area to the actual contact area (A/A c ), the stress concentration along the asperities (i.e., along contact areas) can be accessed. For example, with a 20% contact area, the applied stress is magnified by a factor of five at the asperities. Subsequently, with an applied stress of a few tens Figure 3.10: Comparison of different approaches to obtain the normal stress against fracture closure during the mechanical closure phase for fault F1: the red curve indicates the stress-linear method in the Bandis model; the blue curve indicates the elastic-limit method; and the green curve indicates the elasto-plastic (EP) model using the solution from the boundary element method as in [START_REF] Röttger | contact. engineering-create, analyze and publish digital surface twins from topography measurements across many scales[END_REF].

of MPa, the asperities might reach their elastic limits and enter the plastic domain. Following the 'overlap' model approximation, we assume that the stress on all asperities is homogeneous when reaching the yield stress of the material σ p :

σ asperities n ≈ σ p (3.12)
Therefore, the applied normal stress can be approximated as follows:

σ n ≈ σ p • A c A (3.13)
Considering the relationship between the contact area and fracture closure shown in Fig. 3.6 and assuming that all asperities are at a stress level close to their compressive strength, the applied normal stress can be written as a function of the fracture closure as follows:

σ n ≈ σ p • (1 + erf ((u f -h)/( √ 2h rms ))/2 (3.14)
The behavior of the applied normal stress with closure is then compared to that from the stresslinear approximation (Eq. 3.10). Fig. 3.10 presents the corresponding curves for both approaches. It is found that the best effective yield stress value is 45.5 MPa for fault F1 (similar values are obtained for other faults, as shown in the supplementary material, Fig. 3.21). This effective yield stress is relatively low compared with the common elastic limit of rock (typically lower than the uniaxial compression strength (U CS)). However, it is consistent with the order of U CS if the scale effects of rock strength are considered [START_REF] Thuro | Scale effects in rock strength properties. part 1: Unconfined compressive test and brazilian test[END_REF][START_REF] Zhang | Study of scale effect on intact rock strength using particle flow modeling[END_REF]. Indeed, according to the relationship between U CS and specimen size U CS/U CS c50 = (50/d) 0.18 proposed by [START_REF] Hoek | Underground excavations in rock[END_REF], where U CS c50 is the U CS of a 50-mm diameter specimen and d the diameter, we can estimate U CS ≈ 43 MPa for our model (d is approximated by the block dimension) if we assume U CS c50 is on the order of 200 MPa. Here, U CS and d are in units of MPa and mm, respectively.

The good agreement between the two models illustrates that the normal stiffness during mechanical closure might be simply predicted from a single parameter, i.e., the yield stress σ p , and the evolution of the contact area during closure (Fig. 3.6). In the supplementary material (Fig. 3.22a and Fig. 3.22b), we also check the linear behavior of the stiffness-stress curve obtained from the linear-stress approach and the yield-contact area approach, showing χ = 23.2 m -1 and χ = 23.1 m -1 , respectively. To complete this check, we compute the stiffness-stress curve considering elastoplastic rheology based on a boundary element method (BEM) [START_REF] Röttger | contact. engineering-create, analyze and publish digital surface twins from topography measurements across many scales[END_REF]. We use a yield stress σ p = 46.7 MPa, which provides a very similar evolution of the normal stress as a function of the fracture closure (Fig. 3.10) and a linear relationship between stiffness and stress (see supplementary material, Fig. 3.22c). This, in turn, validates our stress-linear approximation.

Normal fracture stiffness and fault volume

The overall evolution of the normal stress against fracture closure is shown in Fig. 3.5. The normal stress increases with fracture closure during mechanical closure (pink zone), while it remains unchanged during sealing (cyan zone). The numerically computed normal stiffness is plotted as the slope at each closure u f (red lines in Fig. 3.5).

Since fracture closure u f is not easily and directly accessible in the field, we propose in Fig. 3.11 to describe the closure state of the fault in a different way. Specifically, we consider the open space volume of the fault V f , which might be accessible through an injection test of the fault. Fig. 3.6 shows that the fault closure varies between [-3h rms 3h rms ], where h rms is the RMS of the fracture surface topography under free boundary conditions, i.e., before any normal loading. In Fig. 3.23, we show that the fault volume is well approximated by V f ≈ A • (6h rms -u f ), where A is the apparent area of the fracture. Therefore, the fault closure can be written as follows:

u f ≈ 6h rms -V f /A (3.15)
Substituting Eq. ( 3.15) into Eq. ( 3.11) and using χ = 1/h rms yields the following: The evolution of the normal stiffness as a function of the open fault volume for the 6 synthetically generated faults, plotted in a log-lin space in Fig. 3.11, shows that the predicted results (pink curve, Eq. 3.16) are a good proxy of the numerical results (dots). This additionally verified the validity of the linear-stress approach to estimate the normal stress.

κ ≈ κ 0 exp (-V f /V c f ) ( 3

Fracture permeability

Permeability measurement

To compute the fracture permeability, we assume laminar flow conditions and apply Darcy's law to the entire fracture-matrix system [START_REF] Blöcher | Permeability of matrix-fracture systems under mechanical loading-constraints from laboratory experiments and 3-d numerical modelling[END_REF]. Since we impose the same pressure boundary conditions for the rock matrix and the fracture, the pressure at the inlet and that at the outlet of the block are constant. Furthermore, there is almost no change in the total area of the block side A side during the entire simulation. The permeability k of the overall system can then be determined after the steady state as follows:

k = - µ f L V A side p in f -p out f (3.17)
where L is the length of the block and V is the total fluid flux of the block obtained as the mean of the inlet flux Vin and outlet flux Vout , i.e., V = ( Vin + Vout )/2, since they are almost identical (Eq. 3.6). Vin and Vout can be computed as the sum of the inlet flux and the outlet flux of the fracture and the rock matrix, respectively:

Vin = v in f A in f + v in m A in m and Vout = v out f A out f + v out m A out m
, where v f and v m are the fluid velocities of the fracture and matrix, and A f and A m are the areas of the fracture and matrix. The subscripts in and out indicate the inlet and outlet of the corresponding object (fracture or matrix), respectively.

The permeability is assessed by applying a small pressure gradient of 10 Pa/256 m to the block and computing the resulting fluid flux across the fault. Due to the low permeability of the matrix, most of the flow takes place along the open space of the fracture between the contact areas. The distribution of the velocity field in the fault-rock system also indicates that the permeability is dominated by the flow along the fault. Due to the surface roughness and the spatially correlated aperture distribution, a strong channeling effect is observed (with a maximum velocity magnitude of 0.5 m/s at the given pressure gradient of 10 Pa/256 m) along the flow path of the fault (Fig. 3.12). In contrast, the maximum velocity magnitude in the matrix is approximately 1 × 10 -13 m/s. Because of the low flow velocities in the matrix, it takes almost 10 years for the pressure and velocity field in the rock domain to reach a steady state.

In addition, we observe a significant permeability anisotropy (Fig. 3.12). The permeability of fault F1 in the x-direction (perpendicular to the fault slip direction) is smaller than that in the y-direction (perpendicular to the fault slip direction) during closure, showing flow anisotropic be-havior.

Permeability evolution

For either the mechanical or mineral sealing process, we impose a closure of ∼ 0.03 m per step, except when approaching the percolation threshold, extra steps are considered between these steps. In the supplementary material, Fig. 3.24 illustrates some selected steps of the progressive closure (denoted as c0 -c11) as maps of the open fault space for F1, and Fig. 3.25 shows the related permeability in the x-and y-directions against the fault volume at each step of closure for all 6 faults. In terms of orders of magnitude, the fluxes along the x-direction and y-direction for case c0 (Fig. 3.24a) are 4.4 and 8.7 m 3 /s, respectively. The computed permeability is in the range of 8 × 10 -7 m 2 . Due to progressive closure, the fault reaches the flow percolation threshold. For fault F1, the flow percolation threshold in the x-direction is between stages c5 (Fig. 3.24f) and c6 (Fig. 3.24g). The transition of the fracture flow reduces the permeability of the overall system to 1.63 × 10 -15 m 2 , approximately 8 times the matrix permeability 2 × 10 -16 m 2 . By contrast, the permeability in the y-direction is only reduced from 3.77 × 10 -7 m 2 (c0) to 2.3 × 10 -7 m 2 (c6). The percolation threshold in the y-direction is between stages c8 (Fig. 3.24i) and c9 (Fig. 3.24j), which reduces the permeability from 3.16×10 -8 m 2 to 5.26 × 10 -16 m 2 . The transition behavior of the permeability when reaching the flow percolation threshold is also observed for all other faults (Fig. 3.25b-f), with a similar magnitude reduction from ∼ 10 -8 m 2 to ∼ 10 -15 m 2 . In Fig. 3.13, we compile all results in one figure by plotting the permeability against the percolated volume (i.e., the difference between the fault volume above and at the percolation threshold). An overall quadratic evolution of the permeability above the percolation threshold is observed, as follows:

k ≈ (B/A 2 ) • (V f -V perc f ) 2 (3.18)
which is an expected behavior from the linear relationship between the fault volume and the closure (Eq. 3.15) assuming a proxy as a parallel plate model (i.e., a cubic law) where the permeability is related to the square of the fracture aperture, as follows:

h = V f /A ≈ (6h rms -u f ). 1 0 0 1 0 0 0 1 0 0 0 0 1 E -9 1 E -8 1 E -7 1 E -6 1 E -5 1 E -4 P e
ᛓR ၓយ 0 . 9 4 3 7 ᰁ५ E X P _ 2 ( U s e r ) យӛ A * x ^2 নࣞ ၖਾᝐᕎ A 2 . 8 9 2 6 1 E -1 4 ± 5 . 5 3 3 7 E -1 6 R e d u c e d C h i -S q r 3 . 3 9 2 0 9 E -1 R ၓយ(C O D )
0 . 9 2 5 9 2 ᛓR ၓយ 0 . 9 2 5 9 2 )) for the 6 studied faults, which is fitted by k

≈ (B/A 2 ) • (V f -V perc f ) 2
, where A is the total block area and B = 1.24 × 10 -4 is a geometrical prefactor.

Stiffness-permeability relationship

The permeability prediction is of key interest for subsurface reservoir engineering, such as enhanced geothermal systems, particularly when fluid circulation is dominated by a single fault. One way to predict the permeability is by fracture stiffness, which can be obtained by geophysical probing, e.g., seismic wave measurements (Pyrak- Nolte et al., 1990,9). This requires a quantitative relationship between permeability and stiffness. Several previous studies have been devoted to linking permeability and stiffness for lab-scale fractures under normal stress (Pyrak- [START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF][START_REF] Petrovitch | Scaling of fluid flow versus fracture stiffness[END_REF][START_REF] Wang | Development of an empirical model relating permeability and specific stiffness for rough fractures from numerical deformation experiments[END_REF][START_REF] Li | On the relationship between normal stiffness and permeability of rock fractures[END_REF]. Here, based on the expressions established above for normal stiffness vs. fault volume (Eq. 3.16) and permeability vs. fault volume (Eq. 3.18), the stiffness and permeability can be related as follows:

κ = κ 0 exp(-(A k/B + V perc f )/V c f ) (3.19)
where V perc f at the percolation threshold might vary for different faults or in different directions. When we use a mean value of 3000 m 3 (see Fig. 3.15), it represents the numerical data well; see the pink curve in a log-lin plot (Fig. 3.14) for stiffness against permeability.

Due to different roughness scaling properties in the two directions (i.e., H x =0.8 -perpendicular to slip and H y =0.6 -along slip direction), faults are generally more permeable along the ydirection. This is shown in Fig. 3.14, as the data set in the y-direction (squares) lies essentially on top of the x-direction (dots).

Permeability anisotropy at the percolation threshold

Close to the fluid flow percolation threshold, permeability decreases significantly and gradually approaches the matrix permeability (this transition is illustrated in the supplementary material for each fracture, see Fig. 3.25). For the large closure case (small volume), the permeability anisotropy is expected to be at a maximum, as the percolation threshold might be reached in one direction (at

V Pmax f
) but not the other (at V P min f ), emphasizing the importance of permeability contrast.

1 0 0 0 1 0 0 0 0 1 E -0 2 1 E + 0 0 1 E + 0 2 1 E + 0 4 1 E + 0 6 1 E + 0 8 V f P m i n k y / k x F a u l t v o l u m e ( m 3 )
F 1 -M e c a h n i c a l F 2 -M e c h a n i c a l F 3 -M e c h a n i c a l F 1 -S e a l i n g F 2 -S e a l i n g F 3 -S e a l i n g We plot in a log-scale in Fig. 3.15 the ratio of the permeability in the y-direction to the permeability in the x-direction as a function of the fault volume. When faults are below the percolation threshold in both directions because of a large closure (i.e., small fault volume below V P min f ) or above because of a low closure (i.e., large fault volume above V Pmax f ), the ratio k y /k x is close to one. When the permeability shows anisotropy, k y /k x jumps approximately 7 orders of magnitude to ∼ 10 7 . This order is comparable to the ratio of the fracture permeability to the matrix permeability. We also observe that the emergence of the percolation threshold in the x-direction (at V Pmax f ) is not far from the transition from mechanical to sealing closure (transition from red to blue color in Fig. 3.15). This shows that the preset transition from mechanical to chemical closure is in fact related to a strong change of the hydraulic fault properties. In addition, we show that fault volumes with large anisotropy (volumes between percolation thresholds V Pmax f and V P min f ) might cover different ranges, as illustrated by the behaviors of faults F1, F2 and F3, where the domain of high anisotropy is either reduced (F1), intermediate (F3) or large (F2). Note that F4 and F5 are similar to F3, and F6 is similar to F2; therefore, they are not shown here for better readability.

V f P m a x F 1 F 2 F 3

Discussion

Comparison with laboratory and in-situ fracture stiffness and stiffness characteristic measurements

Previous studies have shown that fracture normal stiffness is related to a variety of factors, such as the elastic properties of the host rock, the distribution of the contact area, and the number of asperities (Pyrak- [START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF][START_REF] Hopkins | Fracture stiffness and aperture as a function of applied stress and contact geometry[END_REF][START_REF] Scholz | The mechanics of earthquakes and faulting[END_REF]. It is also essential to incorporate the fracture size when dealing with fracture normal stiffness. This is because rock samples measured in the lab are commonly of small size (∼ mm to ∼ cm), and therefore, there is a large difference in the fracture properties between the lab and the in-situ scale (∼ m to ∼ km). According to the power law scaling of the roughness amplitude for self-affine fracture surfaces [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF], longer fractures typically have wider apertures and larger volumes and thereby become more compliant, as also revealed by previous studies (Pyrak [START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF][START_REF] Hobday | Field measurements of normal and shear fracture compliance[END_REF]. Fig. 3.16 shows an experimental dataset of fracture stiffness (blue dots) derived from a wide range of fracture length scales [START_REF] Hobday | Field measurements of normal and shear fracture compliance[END_REF]. For comparison, we present here the stiffness data from our numerical simulations at the 256 m scale with contact regions of ∼ 20% (red squares) and ∼ 60% (red stars), representing two degrees of fracture closure, i.e., without any sealing at the end of the mechanical closure and with significant sealing, respectively. To illustrate the scale effect, we also conduct numerical modeling on several smaller fractures (i.e., 25.6 m and 0.256 m) generated synthetically following the same power spectral density and RMS scaling as for the 256 m fractures (Fig. 3.3). Although variations in the rock elastic properties are not considered, the overall numerical results of Fig. 3.16 match well with the experimental data and show good consistency with the stiffness-dimension scaling relationship (dashed line) proposed by [START_REF] Worthington | The scaling of fracture compliance[END_REF].

The stiffness characteristic has been found to be rather constant for numerous fractures [START_REF] Zangerl | Normal stiffness of fractures in granitic rock: A compilation of laboratory and in-situ experiments[END_REF][START_REF] Evans | The effects of fracture normal compliance on the long term circulation behavior of a hot dry rock reservoir: a parameter study using the new fully-coupled code'fracture[END_REF], from well-mated fractures in laboratory tests to reactivated hydrofractures in in-situ borehole tests. The values of the stiffness characteristic cover a broad range from 3 mm -1 to 720 mm -1 depending on the fracture area (blue and pink dots in Fig. 3.17). In this study, we obtain the stiffness characteristic ∼ 100 mm -1 for an area of 6.55×10 -2 m 2 (0.256 m scale), ∼ 0.2 mm -1 for an area of 6.55×10 2 m 2 (25.6 m scale) and ∼ 0.02 mm -1 for an area of 6.55×10 4 m 2 (256 m scale), as shown by the red squares in Fig. 3.17. Although [START_REF] Zangerl | Normal stiffness of fractures in granitic rock: A compilation of laboratory and in-situ experiments[END_REF]

concluded 1 0 -2 1 0 -1 1 0 0 1 0 1 1 0 2 1 0 3 1 0 1 1 0 2 1 0 3 1 0 4 1 0 5 1 0 6 1 0 7 1 0 8 ( H o b d a y & W o r t h i n g t o n , 2 0 1 2 ) ( W o r t h i n g t o n
, 2 0 0 7 ) N u m e r i c a l r e s u l t s ( 2 0 % o f c o n t a c t a r e a ) N u m e r i c a l r e s u l t s ( 6 0 % o f c o n t a c t a r e a )

N o r m a l s t i f f n e s s ( M P a / m ) F r a c t u r e s c a l e ( m ) Figure 3.16: Scaling of fracture normal stiffness as a function of fracture length; modified from [START_REF] Hobday | Field measurements of normal and shear fracture compliance[END_REF] and [START_REF] Worthington | The scaling of fracture compliance[END_REF]. that there is only weak evidence of systematic scale effects for fracture cross-sectional areas larger than 100 cm 2 , we obtain obvious scale effects and small variations in the stiffness characteristics of all studied faults at a given scale. This in turn suggests that under certain stress conditions, the prediction of normal stiffness is possible for such faults with self-affinity.

Implications for chemical stimulation

In this study, we divide fracture closures into mechanical and sealing processes and simulate the coupled fluid flow and mechanical response during closure. Compared to previous studies that dealt with permeability and stiffness separately [START_REF] Petrovitch | Scaling of fluid flow versus fracture stiffness[END_REF][START_REF] Wang | Development of an empirical model relating permeability and specific stiffness for rough fractures from numerical deformation experiments[END_REF][START_REF] Li | On the relationship between normal stiffness and permeability of rock fractures[END_REF], our model allows access to both permeability and stiffness simultaneously. Moreover, we consider fluid exchange between the matrix and the fracture, the latter being a more realistic situation in practice compared to other models that generally ignore the matrix effect [START_REF] Méheust | Flow enhancement of a rough fracture[END_REF][START_REF] Pyrak-Nolte | Single fractures under normal stress: The relation between fracture specific stiffness and fluid flow[END_REF][START_REF] Wang | Development of an empirical model relating permeability and specific stiffness for rough fractures from numerical deformation experiments[END_REF]. Although this effect seems to be negligible here when the matrix permeability is relatively low, the topic might be of great interest in future study.

The fracture sealing due to mineral precipitation is assumed to be time-independent and the precipitation growth rate is therefore constant everywhere. However, due to the channeling effect, the fluid velocity along the fracture will be different. This would affect the location of mineral pre-

1 0 1 1 0 2 1 0 3 1 0 4 1 0 5 1 0 6 1 0 7 1 0 8 1 0 9 1 0 -3 1 0 -2 1 0 -1 1 0 0 1 0 1 1 0 2 1 0 3
L a b o r a t o r y t e s t s f r o m l i t e r a t u r e I n -s i t u t e s t s f r o m l i t e r a t u r e R e s u l t s f r o m t h i s s t u d y S t i f f n e s s c h a r a c t e r s t i c ( m m -1 ) F r a c t u r e a r e a ( c m 2 ) Figure 3.17: Fracture stiffness characteristic as a function of fracture area; modified from [START_REF] Zangerl | Normal stiffness of fractures in granitic rock: A compilation of laboratory and in-situ experiments[END_REF]. cipitation and change the permeability [START_REF] Ortoleva | Geochemical self-organization i; reactiontransport feedbacks and modeling approach[END_REF][START_REF] Sausse | Hydromechanical properties and alteration of natural fracture surfaces in the soultz granite (bas-rhin, france)[END_REF][START_REF] Wetzel | Hydraulic and mechanical impacts of pore space alterations within a sandstone quantified by a flow velocity-dependent precipitation approach[END_REF]. Time-dependent sealing may result in completely sealed fractures or partially sealed fractures at the same scale [START_REF] Griffiths | Geothermal implications for fracture-filling hydrothermal precipitation[END_REF], i.e., the mineral growth might be nonhomogeneous in the fracture space. Moreover, previous studies on the alteration of fracture permeability by coupled mechanical and chemical processes have shown that contact regions can dissolve faster than the open fracture space under high flow rate conditions (e.g., [START_REF] Ameli | Fracture permeability alteration due to chemical and mechanical processes: a coupled high-resolution model[END_REF]). The nonhomogeneous precipitation/dissolution would modify the percolation threshold of the fracture flow in both directions compared to the homogeneous precipitation/dissolution. This would require the development of a coupled H-M-C model. Despite these assumptions and given the low flow rate assumed, the main conclusions derived in this study remain valid.

In addition, we only consider mineral precipitation inside the fracture. Previous studies have shown that it might also occur in the pores [START_REF] Griffiths | Geothermal implications for fracture-filling hydrothermal precipitation[END_REF]. This may modify the permeability of the rock matrix and therefore the fluid exchange between the fracture and the matrix. However, the bulk effect would primarily lower the effective permeability of the matrix, so our conclusions, which are based on the assumption that the fracture flow dominates over the bulk matrix flow, will remain unaltered. In this case, the pore filling may have a limited influence on the whole permeability.

Other factors, such as temperature, may directly affect the mineral growth rate during the sealing process [START_REF] Griffiths | Geothermal implications for fracture-filling hydrothermal precipitation[END_REF]. Accounting for these additional effects requires a fully coupled T-H-M-C model, which will be considered in our future studies.

Finally, we can identify several implications in terms of permeability anisotropy during closure for chemical stimulation of faults. First, if the fault volume is larger than the percolation volume V pmax f , that is, the fluid channels are connected in both directions, a chemical treatment might provide only a limited improvement in permeability. Second, if the fault volume is lower than the percolation volume V p min f or if it is between the percolation volume V pmax f and V p min f , that is, the fluid channels are fully sealed or only connected in one direction, a chemical treatment may significantly enhance the permeability in at least one direction. In the latter case, one should also consider how thick the minerals are that need to be dissolved to improve the permeability. It is neither economically nor technically available if too much chemical solvent is required. For example, according to Fig. 3.15, if V = 400 m 3 , we only need to dissolve a small layer of minerals (∼ 0.01 m) for F2, but it requires three times or more for F1 and F6 to improve their permeability noticeably. In comparison, if V = 4000 m 3 , the chemical treatment might be available for all three faults (the dissolved thickness < 0.01 m). The fault volume, which is directly related to the degree of sealing, might be used as a key parameter in the chemical treatment. Moreover, it can be determined once the normal stiffness is known because of their quantitative function relationship (Fig. 3.11). The normal stiffness can in turn be derived relatively easily by field measurements. Therefore, through stiffness and permeability tests, it is possible to estimate the degree of sealing and evaluate a chemical treatment.

However, one should also note that if a critical stiffness is exceeded during the chemical treatment, i.e., the lowest stiffness to maintain fault stability, faults are likely to be closed. Here, the critical state is defined as the transition from the mechanical to the sealing regime at 20% of the contact area (see Fig. 3.15). When close to the transition zone, chemical stimulation should be avoided since fault collapse may in turn reduce the permeability. This critical state requires further evaluation for faults in a given environment. In summary, the degree of sealing, anisotropic flow behavior and critical stiffness are the key factors in the chemical stimulation of reservoir connectivity dominated by a single fault.

Conclusion

By means of a fully 3-dimensional finite element approach (MOOSE/GOLEM), we simulated the hydro-mechanical behavior of a partly sealed rough fault embedded in a granite reservoir. Here, faults with an aspect ratio of 1:1000 were handled. The progressive closure of the fracture was modeled in two successive phases. First, a mechanical closure following a specific rigid-plastic model that we bench-marked by comparing the following three approaches: one based on a stress-linear stiffness approximation as broadly empirically observed, one based on a yield contact area assessment and one on an independent boundary element method to reproduce an elastoplastic closure. Second, a sealing of the fault was simulated by adding a stepwise coating to the fault surfaces, which progressively reduces the fault aperture. We defined the switch from one phase to the other at approximately 20% of the contact area based on experimental observations. It appears from our computation that this limit also corresponds to a specific closure of the fracture (approximately h rms of the fracture roughness) and to the flow percolation threshold along the fault slip perpendicular direction. At the percolation threshold, there is a dramatic reduction in permeability by ∼ 8 orders of magnitude with a strong anisotropy. We also showed that the evolution of the fracture stiffness during closure follows a generic exponential law with the fracture volume, including a characteristic volume that is related to features of the fracture geometry and stiffness characteristic classically measured during mechanical tests. Moreover, a stiffness-permeability relationship has been obtained in terms of fault volume.

The success of a chemical treatment to increase the permeability of a partly sealed fault strongly depends on the degree of sealing. If flow paths are fully sealed at least in one direction or anisotropic flow behaviour is observed (a directional percolation threshold exists), a chemical treatment may have a major impact but also depending on how close it is to the percolation threshold. Otherwise, if flow paths are open in all directions, a chemical treatment will have a minor impact. To quantify the degree of sealing, the normal stiffness can be used. The stiffness of the fault is a function of fault volume and, henceforth, a function of the degree of sealing. If the normal stiffness can be derived by field measurements, the degree of sealing can be estimated, and the success of a chemical treatment can be evaluated. In summary, chemical treatments for fault-rock systems should consider multiple factors, such as the directional percolation threshold, anisotropic permeability and appropriate fracture stiffness. Meanwhile, the collapsing of the fault by chemical stimulation should be avoided when approaching the mechanical-sealing transition state. Finally, we made some assumptions in this study for mineral sealing and ignored some effects, such as time and temperature, which could be improved in the future. ing elasto-plastic contact problems on rough surfaces and for fruitful discussions with Prof. Lars Pastewka and Mr. Antoine Sanner. This work was performed under the framework of the Laboratory of Excellence LABEX ANR-11-LABX-0050-G-EAU-THERMIE-PROFONDE and the Interdisciplinary Thematic Institute GeoT, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm, supported by IdEx Unistra (ANR 10 IDEX 0002) and by SFRI STRAT'US project (ANR 20 SFRI 0012) under the framework of the French Investments for the Future Program.

Complementary materials

1 0 -4 1 0 -3 1 0 -2 1 0 -1 1 0 0 1 0 1 1 0 -1 0 1 0 -5 1 0 0 1 0 5 1 0 1 0 1 0 1 5 X -d i r e c t i o n Y -d i r e c t i o n P ⊥ ( j ) = 3 . 5 ×1 0 -4 j -2 . 6 P // ( j ) = 3 . 5 ×1 0 -4 j -2 . 2
P o w e r s p e c t r a l d e n s i t y P ( j ) ( m 3 ) W a v e n u m b e r j ( m -1 ) andY). The dashed lines indicate the domain of measured surfaces from the field and laboratories [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF] (top line is C(q) = 2 × 10 -2 q -3 , medium line is C(q) = 2 × 10 -4 q -3 and bottom line is C(q) = 2 × 10 -6 q -3 ); (b) RMS h rms evolution with scale L along X and Y for the same synthetic fault.

1 0 0 1 0 1 1 0 2 1 0 3 1 0 -3 1 0 -2 1 0 -1 1 0 0 R M S h
1 E -1 1 E -1 1 E -1 1 E -1 1 E - 1 E - 1 E - (

Chapter 4

Hydro-mechanical responses during hydraulic fracturing

Résumé

Nous développons un modèle hydro-mécanique entièrement couplé pour simuler le glissement de faille dû à l'injection de fluide. Nous considérons l'interaction entre une fracture hydraulique et les failles préexistantes ainsi que l'échange de fluide entre la fracture/la faille et la matrice poreuse. Afin de considérer un mécanisme de diffusion de la pression, nous fixons une perméabilité relativement élevée autour du chemin stimulé. Notre étude paramétrique montre que plusieurs facteurs affectent l'activation de la faille et son comportement de glissement, tels que les propriétés de la faille, les propriétés de friction et le scénario d'injection. Nous observons que la diffusion de la pression interstitielle induit un changement de contrainte poroélastique, capable de produire des événements de fermeture avec un décalage temporel et spatial. Ce mécanisme affecte également le comportement de glissement pendant l'injection, en particulier lorsque la perméabilité environnante est élevée (par exemple, jusqu'à 1e-13 m 2 ), et fournit un nouvel éclairage pour comprendre l'apparition d'événements sismiques plus forts après la fermeture par rapport à la phase d'injection. En outre, nous montrons que de petites perturbations peuvent déclencher un glissement de faille sismique important, ce qui souligne le rôle clé de l'état de contrainte initial de la faille. Les résultats ont des implications profondes pour l'ingénierie liée à l'injection de fluides en profondeur ainsi que pour les stratégies d'injection cyclique douce visant à atténuer le risque de grands séismes.

Introduction

Hydraulic stimulation has been applied in many underground engineering projects, such as coalbed methane mining [START_REF] Mcdaniel | Hydraulic fracturing techniques used for stimulation of coalbed methane wells[END_REF][START_REF] Colmenares | Hydraulic fracturing and wellbore completion of coalbed methane wells in the powder river basin, wyoming: implications for water and gas production[END_REF], oil/gas recovery [START_REF] Law | Gas in tight reservoirs-an emerging major source of energy[END_REF][START_REF] Barati | A review of fracturing fluid systems used for hydraulic fracturing of oil and gas wells[END_REF], shale gas production [START_REF] Johri | The evolution of stimulated reservoir volume during hydraulic stimulation of shale gas formations[END_REF][START_REF] Lei | Fault reactivation and earthquakes with magnitudes of up to mw4. 7 induced by shale-gas hydraulic fracturing in sichuan basin, china[END_REF], waste water disposal [START_REF] Keranen | Sharp increase in central oklahoma seismicity since 2008 induced by massive wastewater injection[END_REF][START_REF] Brudzinski | Seismicity induced by hydraulic fracturing and wastewater disposal in the appalachian basin, usa: A review[END_REF] as well as deep geothermal energy extraction [START_REF] Fehler | Stress control of seismicity patterns observed during hydraulic fracturing experiments at the fenton hill hot dry rock geothermal energy site, new mexico[END_REF][START_REF] Murphy | Hydraulics and well testing of engineered geothermal reservoirs[END_REF][START_REF] Evans | Microseismicity and permeability enhancement of hydrogeologic structures during massive fluid injections into granite at 3 km depth at the soultz hdr site[END_REF][START_REF] Gaucher | Induced seismicity in geothermal reservoirs: physical processes and key parameters[END_REF][START_REF] Schill | Hydraulic performance history at the soultz egs reservoirs from stimulation and long-term circulation tests[END_REF]. It has been proven to be effective to improve fracture connectivity, increase reservoir transmissivity and enhance fluid circulation. However, due to high pressure fluid injection, such operation is often accompanied by seismic activities. Sometimes they can be felt by people on the ground and destroy ground facilities, therefore cause strong social disputes [START_REF] Majer | Induced seismicity associated with enhanced geothermal systems[END_REF][START_REF] Häring | Characterisation of the basel 1 enhanced geothermal system[END_REF][START_REF] Langenbruch | Physics-based forecasting of man-made earthquake hazards in oklahoma and kansas[END_REF][START_REF] Lee | Managing injection-induced seismic risks[END_REF][START_REF] Lei | The december 2018 ml 5.7 and january 2019 ml 5.3 earthquakes in south sichuan basin induced by shale gas hydraulic fracturing[END_REF]. Furthermore, those man-made earthquakes not only occur during the injection phase, but also after the operation termination, termed as 'post injection induced seismicity', which can be observed both in the laboratory and on the field [START_REF] Stanchits | Fracturing of porous rock induced by fluid injection[END_REF]Hsieh and Bredehoeft, 1981;[START_REF] Evans | Microseismicity and permeability enhancement of hydrogeologic structures during massive fluid injections into granite at 3 km depth at the soultz hdr site[END_REF][START_REF] Albaric | Monitoring of induced seismicity during the first geothermal reservoir stimulation at paralana, australia[END_REF]. The post injection induced seismicity sometimes exhibits discontinuities both in time and space. For instance, the hydraulic stimulation performed in June 2013 in the Rittershoffen geothermal site induced seismicity during injection but a second swarm of events occurred over 100 meters away from the first earthquake sequence 4 days after shut-in [START_REF] Lengliné | Seismicity related to the hydraulic stimulation of grt1, rittershoffen, france[END_REF]. Moreover, these additional seismic events may have large magnitude, some of them are even stronger than the seismicity during the injection [START_REF] Häring | Characterisation of the basel 1 enhanced geothermal system[END_REF][START_REF] Baisch | A numerical model for fluid injection induced seismicity at soultz-sous-forêts[END_REF][START_REF] Albaric | Monitoring of induced seismicity during the first geothermal reservoir stimulation at paralana, australia[END_REF][START_REF] Mcclure | Generation of large postinjection-induced seismic events by backflow from dead-end faults and fractures[END_REF][START_REF] Baujard | Hydrothermal characterization of wells grt-1 and grt-2 in rittershoffen, france: Implications on the understanding of natural flow systems in the rhine graben[END_REF]. Those observations increase the uncertainty and complexity of the hydraulic operation response, and pose a challenge for earthquake risk prediction and control.

Over the past few decades, a large number of researches have been devoted to studying fluidinduced seismicity mechanisms and its hazard control, across scales ranging from ∼ cm to ∼ km. Those studies are based on either a statistical [START_REF] Shapiro | Large-scale in situ permeability tensor of rocks from induced microseismicity[END_REF][START_REF] Parotidis | Back front of seismicity induced after termination of borehole fluid injection[END_REF][START_REF] Barth | Probability of earthquake occurrence and magnitude estimation in the post shut-in phase of geothermal projects[END_REF] or physical description of the induced dynamics. Physics based investigation comprises laboratory experiments [START_REF] Stanchits | Fracturing of porous rock induced by fluid injection[END_REF][START_REF] Passelègue | Fault reactivation by fluid injection: controls from stress state and injection rate[END_REF][START_REF] Ji | Unloading-induced rock fracture activation and maximum seismic moment prediction[END_REF][START_REF] Wang | Injection-induced seismic moment release and laboratory fault slip: Implications for fluid-induced seismicity[END_REF], field tests [START_REF] Cornet | Analysis of induced seismicity for stress field determination and pore pressure mapping[END_REF][START_REF] Blöcher | Evaluating micro-seismic events triggered by reservoir operations at the geothermal site of groß schönebeck (germany)[END_REF][START_REF] Davis | Energy of injectioninduced seismicity predicted from in-situ experiments[END_REF], as well as numerical modeling [START_REF] Baisch | A numerical model for fluid injection induced seismicity at soultz-sous-forêts[END_REF][START_REF] Yoon | Discrete bonded particle modelling of fault activation near a nuclear waste repository site and comparison to static rupture earthquake scaling laws[END_REF][START_REF] Maurer | Role of fluid injection on earthquake size in dynamic rupture simulations on rough faults[END_REF]. The main triggering mechanisms for injection-induced seismicity, can be attributed to pore pressure diffusion [START_REF] Shapiro | Large-scale in situ permeability tensor of rocks from induced microseismicity[END_REF][START_REF] Brown | Small earthquakes matter in injection-induced seismicity[END_REF], poroelastic stress response [START_REF] Mcclure | Investigation of injection-induced seismicity using a coupled fluid flow and rate/state friction model[END_REF][START_REF] Chang | Injection-induced seismicity on basement faults including poroelastic stressing[END_REF][START_REF] Jacquey | Far field poroelastic response of geothermal reservoirs to hydraulic stimulation treatment: Theory and application at the groß schönebeck geothermal research facility[END_REF], aseismic slip (Guglielmi et al., 2015a;[START_REF] Lengliné | Seismicity related to the hydraulic stimulation of grt1, rittershoffen, france[END_REF], as well as earthquake interactions [START_REF] Yeo | Causal mechanism of injection-induced earthquakes through the m w 5.5 pohang earthquake case study[END_REF]. Yet other mechanisms are also reported to be responsible for post-injection induced seismicity, such as the superposed stress variations [START_REF] De Simone | Superposition approach to understand triggering mechanisms of post-injection induced seismicity[END_REF] and fracture normal closure [START_REF] Ucar | Postinjection normal closure of fractures as a mechanism for induced seismicity[END_REF].

Although those studies have formed the basic framework of fluid injection induced seismicity, none of them can deal with all cases due to the complexity of the underground structure and the coupling process. Some questions, e.g., the role of related control parameters, have been partially answered and still need to be further investigated. Thus, it is of primary importance to explore the multi-physical processes occurring in deep reservoirs from the injection stage to the post injection stage, in order to minimize the risk of earthquakes and achieve successful hydraulic stimulation.

The goal of hydraulic stimulation is to stimulate permeability enhancement in the reservoir. As such it requires a good conceptual model for the stimulation. Currently, there are four main stimulation mechanisms are considered: (i) pure opening (tensile) mechanism, aiming to create new fractures (i.e., hydro-fracturing) [START_REF] Zoback | Laboratory hydraulic fracturing experiments in intact and pre-fractured rock[END_REF][START_REF] Ren | Simultaneous hydraulic fracturing of ultra-low permeability sandstone reservoirs in china: Mechanism and its field test[END_REF]; (ii) pure shear mechanism, which assumes that the stimulation mainly acts on the shear slip of the pre-existing fractures [START_REF] Xie | Initiation and propagation of fracture shearing during hydraulic stimulation in enhanced geothermal system[END_REF][START_REF] Ye | Injection-induced shear slip and permeability enhancement in granite fractures[END_REF]; (iii) primary fracturing with shear stimulation leakoff, under which new fractures are continuously created and propagated from the wellbore with fluid leaks off into natural fractures [START_REF] Pearson | The relationship between microseismicity and high pore pressures during hydraulic stimulation experiments in low permeability granitic rocks[END_REF][START_REF] Wang | Hydraulic fracturing with leakoff in a pressuresensitive dual porosity medium[END_REF]; (iv) the mixed mechanism, i.e., flow pathways are connected with both newly created and pre-existing fractures [START_REF] Norbeck | An embedded fracture modeling framework for simulation of hydraulic fracturing and shear stimulation[END_REF][START_REF] Lei | Modelling fluid injection-induced fracture activation, damage growth, seismicity occurrence and connectivity change in naturally fractured rocks[END_REF]. Different communities tend to rely on a different stimulation concept. For example, the opening mode is commonly the preferred stimulation in oil/gas industry, while hydro-shearing is widely used in the context of Enhanced Geothermal System (EGS). Numerous field observations are also supportive of the mixed stimulation [START_REF] Albaric | Monitoring of induced seismicity during the first geothermal reservoir stimulation at paralana, australia[END_REF][START_REF] Norbeck | Field observations at the fenton hill enhanced geothermal system test site support mixed-mechanism stimulation[END_REF][START_REF] Krietsch | Changing flow paths caused by simultaneous shearing and fracturing observed during hydraulic stimulation[END_REF]. Moreover, a recent study suggests an important role of the stress transfer from hydraulic fracture opening in induced seismicity distribution [START_REF] Kettlety | Stress transfer from opening hydraulic fractures controls the distribution of induced seismicity[END_REF].

In the context of the mixed stimulation strategy, there are two processes that could trigger seismicity, that is, the brittle failure of intact rocks by hydro-fracturing (HF) and the activation of pre-existing fractures. It has been evidenced that HF-induced seismicity in terms of magnitude is negligible compared to the fault shearing induced seismicity [START_REF] Zoback | Reservoir geomechanics[END_REF][START_REF] Lei | Modelling fluid injection-induced fracture activation, damage growth, seismicity occurrence and connectivity change in naturally fractured rocks[END_REF]. Hence, many studies only focus on the shear slip of natural fractures in fluid injection-induced seismicity [START_REF] Mcclure | Investigation of injection-induced seismicity using a coupled fluid flow and rate/state friction model[END_REF][START_REF] Rutqvist | Modeling of fault reactivation and induced seismicity during hydraulic fracturing of shale-gas reservoirs[END_REF][START_REF] Yoon | Discrete bonded particle modelling of fault activation near a nuclear waste repository site and comparison to static rupture earthquake scaling laws[END_REF].

In this study, we aim at investigating fault slip behavior during active fluid injection and after shut-in. We develop a fully hydro-mechanical coupling model for a fractured porous rock while considering the interaction between existing faults and the hydraulic fractures. In doing so, we consider only a single path for the stimulated fracture thereby being able to better control the model behavior. In section 4.2, basic governing equations are described. Section 4.3 introduces the model setup for the simulations, which is followed by section 4.4 regarding results with parametric studies on the matrix permeability, the fault friction coefficient, fault orientation, as well as different injection scenarios. Further discussion and main conclusions are given in section 4.5 and section 4.6, respectively. A validation of hydraulic fracture propagation is also presented in the end.

A coupled hydro-mechanical cohesive zone model

Generally, when involving hydraulic fracturing, there are four physical processes to be considered, i.e., porous rock mass deformation, pore fluid flow, fracture fluid flow and fracture propagation [START_REF] Carrier | Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model[END_REF]. For problems which involve interactions with pre-existing faults, the fault slip should also be defined [START_REF] Chen | Finite-element simulation of a hydraulic fracture interacting with a natural fracture[END_REF][START_REF] Li | Modeling hydraulic fracture propagation using cohesive zone model equipped with frictional contact capability[END_REF].

Porous media deformation

Under isothermal conditions, rock mass is considered as an isotropic, poroelastic material. Assuming small strains, the poro-elastic constitutive relation is expressed as [START_REF] Coussy | Poromechanics[END_REF]:

σ ij -σ ij,0 = 2Gε ij + K - 2 3 G ϵδ ij -b(p -p 0 )δ ij (4.1)
where σ ij and ε ij are the total stress and strain; ϵ the volumetric strain; p the pore pressure; G and K are the dry elastic shear and bulk moduli, and b the Biot's coefficient. The subscript '0' represents the initial state of each variable.

Eq. ( 4.1) in terms of effective stress can be simplified to:

σ ′ ij -σ ′ 0 ij = 2Gε ij + K - 2 3 G ϵδ ij (4.2)
Following Biot's theory [START_REF] Biot | General theory of three-dimensional consolidation[END_REF], the effective stress tensor is defined for fully saturated media as:

σ ′ ij = σ ij -bpδ ij (4.3)

Pore fluid flow

It is assumed that fluid flows through an interconnected pore network following Darcy's law:

q = - k η ∇p (4.4)
where q is the Darcy's velocity, η is the pore fluid dynamic viscosity, and k is the permeability.

The continuity flow equation for the fluid within the pores is given by:

1 M b ∂p ∂t + b ∂ϵ ∂t + ∇ • q = 0 (4.5)
where M b is Biot's modulus defined by:

1 M b = ϕ 0 K f + b -ϕ 0 K s (4.6)
where K f , K s denotes the pore fluid bulk modulus and the porous medium solid grain bulk modulus, and ϕ 0 is the initial porosity.

Combining Eq. ( 4.4) with Eq. ( 4.5), the pore fluid diffusion equation is obtained as follows:

1 M b ∂p ∂t + b ∂ϵ ∂t = k η ∇ 2 p (4.7)

Fluid flow inside the fracture

Assuming Newtonian and incompressible fluid, the flow inside a fracture can be split into a tangential flow within the fracture and a normal flow across the fracture walls, as shown in Fig. 4.1. The tangential flow is governed by Reynold's lubrication theory defined by the continuity equation [START_REF] Detournay | Propagation regimes of fluid-driven fractures in impermeable rocks[END_REF]:

∂w ∂t + ∂q f ∂x + v = 0 (4.8)
where w is the fracture aperture, q f is the longitudinal fluid flow rate, and v denotes the normal flow rate leaking into the porous medium through the fracture surfaces. The latter parameter can be calculated from the difference between the fluid pressure inside the fracture p f and the pore fluid pressure:

v = c (p f -p) (4.9)
where c is the leak-off coefficient. Eq. ( 4.9) imposes a pressure continuity between the fracture and the rock and defines a pressure-flow relationship between the cohesive element's middle nodes and its adjacent surface nodes. In homogeneous porous reservoirs, leakoff only occurs normal from the fracture into the pore system. This kind of leakoff can be described by Carter's leakoff model considering a filter-cake zone, an invaded zone, and a reservoir compaction zone [START_REF] Howard | Optimum fluid characteristics for fracture extension[END_REF][START_REF] Liu | Leakoff characteristics and an equivalent leakoff coefficient in fractured tight gas reservoirs[END_REF]. In such a model the leakoff coefficient as shown in Eq. ( 4.9) is pressure independent. In the case of additional natural fractures, the physics of leakoff will alter due to pressure dependent flow behavior of the fracture. In that case the leakoff coefficient becomes pressure dependent [START_REF] Liu | Leakoff characteristics and an equivalent leakoff coefficient in fractured tight gas reservoirs[END_REF]. In our study we only consider the first scenario, the pressure diffusion from a fracture into a porous matrix and henceforth only consider a pressure independent leakoff coefficient (Eq. 4.9). In addition, as leak-off coefficient and permeability increase, the fracture length decreases [START_REF] Yao | Linear elastic and cohesive fracture analysis to model hydraulic fracture in brittle and ductile rocks[END_REF]. When the permeability and the leak-off coefficient are sufficient small, the effect on the fracture length will be negligible. In this work, considering the range of permeability (≤ 1e-13 m 2 ), we set the leak-off coefficient 1e-11 m 3 /kPa•s, e.g., [START_REF] Yao | Linear elastic and cohesive fracture analysis to model hydraulic fracture in brittle and ductile rocks[END_REF].

Assuming the fracture walls are parallel and smooth, the longitudinal fluid flow rate q f is related to the pressure gradient inside the fracture as (the Poiseuille law):

q f = - w 3 12η ∂p f ∂x (4.10)
where q f is the flow flux inside the fracture.

Substituting Eq. ( 4.9) and Eq. ( 4.10) into Eq. ( 4.8) yields the governing equation of the fluid flow inside the fracture:

∂w ∂t + v = ∂ ∂x w 3 12η ∂p f ∂x (4.11)

Fracture initiation and propagation

For hydraulic fracturing, stress-pressure cohesive elements coupled with displacement and pore pressure degrees of freedom (DOFs) are simultaneously used. They are embedded in solid continuum elements (also with pressure DOF) to define a predefined crack path, such that the fracture growth is constrained to this path. Cohesive elements have been widely used to study fracture problems in rock-like materials [START_REF] Zhuang | A comparative study on unfilled and filled crack propagation for rock-like brittle material[END_REF][START_REF] Saadat | A cohesive discrete element based approach to characterizing the shear behavior of cohesive soil and clay-infilled rock joints[END_REF]. During the fracture propagation, a small fracture process zone (FPZ) is formed to capture stress concentration around the crack tip. In this study, we assume the damage initiation and evolution follow a linear tractionseparation law (Fig. 4.2). The relation between the traction (strength) T and the separation δ is given by [START_REF] Yao | Linear elastic and cohesive fracture analysis to model hydraulic fracture in brittle and ductile rocks[END_REF]:

T = (1 -D)K 0 δ (4.12)
where D(0 ≤ D ≤ 1) is the damage variable and K 0 is a penalty stiffness which can avoid total compliance of the whole model before damage initiation [START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF]. The total separation (i.e., fracture propagation) is dependent on the fracture energy G c and the damage initiation is determined by the maximum nominal stress criterion:

max < t n > t 0 n , t s t 0 s , t t t 0 t = 1 (4.13)
where t 0 n , t 0 s , t 0 t is the peak value of the nominal stress when the deformation is either purely normal to the interface or purely in the first or the second shear direction, respectively. The Macaulay bracket symbol "<>" signifies that a pure compressive stress does not initiate damage.

Fault slip

The Coulomb friction law is used to model the shear slip behavior of fault. Assuming the cohesion strength of the fault surface is negligible, i.e., for a cohesion-less fault, the shear strength τ c is changing with the normal stress σ n as follows:

τ c = µ • (σ n -p) (4.14)
where µ is the coefficient of friction. Typically, µ varies between 0.6 and 1 for intact rock, 0.3 and 0.6 or even lower for pre-existing faults [START_REF] Zoback | Reservoir geomechanics[END_REF][START_REF] Ellsworth | Injection-induced earthquakes[END_REF][START_REF] Gaucher | Induced seismicity in geothermal reservoirs: physical processes and key parameters[END_REF]. In this work, µ = 0.2, 0.4 and 0.6 are selected for parametric study.

As soon as the shear stress exceeds the fault shear strength, the fault becomes unstable and slips along the fault plane. There are several causative processes responsible for the shear slip. For example, a reduction in the effective normal stress, or a decrease in the friction coefficient.

The induced shear slip is often accompanied by the release of stored strain energy, hence induces seismicity.

Model setup

We consider hydraulic stimulation based on the assumption of mixed mechanism (i.e., both fracture opening and shearing). As shown in Fig. 4.3(a), we define two pre-existing faults F1 and F2 (blue solid line) in a 2D model. Those two faults are connected with potential hydraulic fracturing paths (blue dashed line), which are assumed to be parallel to the maximum principle stress σ 1 [START_REF] Zoback | Reservoir geomechanics[END_REF] and perpendicular to the minimum principle stress σ 3 . Each individual segment has a length of 10 m and both faults are oriented at a degree of α, where α is the angle between the maximum principal stress σ 1 and the fault plane. All boundaries are fixed with zero displacement and initial pore pressure of the model. The injection point is located at the bottom of the model (red dot) and fluid is injected at a flow rate Q. To investigate the after shut-in effect, we intentionally set the stop position between F1 and F2 (black dot). It will not exceed the half length of the segment (5 m) from F1 in all the simulations. To precisely describe the pore pressure diffusion, we set a relatively permeable zone near the faults and hydraulic paths (black dashed lines) for which we assumed a constant width (1 m). It is modeled by an increased permeability (considered as an effective value averaging local micro-structural heterogeneities) than the intact rock matrix, where the fracture density is high. In addition, we also consider an exchange of fluid mass from the fractures to the permeable zone, in order to explore the impact of its permeability (linked with pressure diffusion) on the fault slip.

At the intersection of hydraulic fracturing path and fault, we assume that the fluid is directly diverted to the fault. This assumption is supported by previous studies which indicated that for a ratio of the initial maximum to the minimum stress larger than 1, and for a friction coefficient smaller than 0.65 and fracture orientation lower than 45°, a hydraulic fracture will be diverted into the pre-existing fault [START_REF] Gu | Hydraulic fracture crossing natural fracture at nonorthogonal angles: a criterion and its validation[END_REF].

In addition, the porous medium is assumed to be fully saturated. The simulations are performed in Abaqus package. All parameters used are listed in Tab. 4.1.

Results and analysis

Fig. 4.4 shows the fluid pressure distribution at different time, and Fig. 4.5 presents the fluid pressure evolution at the injection point and the center of F1, F2. The fault F1 is activated at 38.1 s while F2 is triggered at 520.1 s (i.e., 360.1 s after shut-in). This can be identified by the damage variable (Fig. 4.3b).

At the beginning, due to high-pressure fluid stimulation, the injection pressure rises to approximately 67 MPa from an initial value 23.7 MPa in a very short time (≈ 2.7 s). Due to this drastic pressure change, the induced tensile stress is sufficient to overcome the rock tensile strength. As the hydraulic fracture propagates, the fluid pressure tends to drop to a stable level. Meanwhile, the fluid pressure at the center of F1 and F2 increases, since the high-pressure fluid inside the hydraulic fracture leaks into the surrounding porous rock mass and then diffuses all around. At 38.1 s (point A in Fig. 4.5), the injected water diverts to F1, which lowers the injected pressure and forces the fluid pressure of F1 to increase. This is accompanied by a shear slip along F1 as its shear strength is overcome by the high pressure (see Fig. 4.8a). Meanwhile, the pressure in the center of F2 witnesses a slight jump due to the poro-elastic effect. Immediately afterwards, the pressure within F1 also decreases due to the generation of a new hydraulic fracture from the end of F1. This pressure drop in regions near fracture intersections is a common feature which has been already observed by other authors, e.g., [START_REF] Piris | Fluid pressure drops during stimulation of segmented faults in deep geothermal reservoirs[END_REF]. From this time onward, the fluid pressures at injection point and F1 almost remain at the same level with a tiny difference owing to the fluid flowing into the matrix. After 160 s of injection (point B in Fig. 4.5), the pumping stops and the fluid pressure starts to decrease. Pore pressure diffusion away from the pressurized hydraulic fracture can indeed be clearly observed (Fig. 4.6a). After approximately 360.1 s of termination of the injection, the pressure perturbation and the corresponding stress change promote instability along the fault. As the fault slips, a compressional zone and a dilation zone emerge around the fault tip (Fig. 4.6). High pressure and stress in the compressional zone and low pressure and stress in the dilation zone are observed, respectively. This leads to only a minor increase in pressure at the F2 center (point C in Fig. 4.5), but to a rather significant increase at Node 2 within the compressional zone (point C in Fig. 4.7). See also Fig. 4.6 for their relative locations. This poro-elastic stressing mechanism is similar to the one described by injection into a single fracture embedded in a poro-elastic medium [START_REF] Lei | Modelling fluid injection-induced fracture activation, damage growth, seismicity occurrence and connectivity change in naturally fractured rocks[END_REF], and the combined induced fault slip mechanism (i.e., interaction between delayed pressure diffusion and poro-elasticity) during shut-in is in agreement with the results from previous investigation, e.g., [START_REF] Chang | Seismicity rate surge on faults after shut-in: Poroelastic response to fluid injection[END_REF]. In fact, the sudden elevated fluid pressure, either owing to direct injection in F1 or pressure diffusion/poro-elasticity for F2, counteracts the normal stress and induces dynamic slip of the fault, which correlates to potential micro-seismic events. The evolution of slip and slip rate of F1 and F2 over time is presented in Fig. 4.8. The characterization of displacement jump and slip rate rise is similar to an unloading experiment of a sawcut fracture [START_REF] Ji | Unloading-induced rock fracture activation and maximum seismic moment prediction[END_REF]. It can also be seen that in most of time the slip rates of F1 and F2 are close to zero, indicating an aseismic behavior. The seismic and aseismic slip is distinguished with a slip velocity threshold 5 mm/s [START_REF] Mcclure | Investigation of injection-induced seismicity using a coupled fluid flow and rate/state friction model[END_REF][START_REF] Cappa | On the relationship between fault permeability increases, induced stress perturbation, and the growth of aseismic slip during fluid injection[END_REF]. In the following, the seismic slip and the associated average slip rate are calculated for each case.

Permeability

Pore pressure diffusion in a fractured porous medium is of central importance in fluid injection in faulted reservoirs [START_REF] Chang | Injection-induced seismicity on basement faults including poroelastic stressing[END_REF][START_REF] Shapiro | Large-scale in situ permeability tensor of rocks from induced microseismicity[END_REF]. The core controlling parameter is the hydraulic diffusivity, which is related to the permeability [START_REF] Deng | Hydraulic diffusivity of a partially open rough fracture[END_REF][START_REF] Shapiro | Large-scale in situ permeability tensor of rocks from induced microseismicity[END_REF]. In order to investigate the impact of the hydraulic pressure diffusion on the fault slip, we carry out a systematic analysis on the role of the permeable zone permeability by varying its value within three steps, that is, k f 1 = 1e-15 m 2 , k f 2 = 5e-15 m 2 and k f 3 = 1e-14 m 2 . We perform several simulations at each permeability by altering the injection duration, in order to capture the activation It is worth noting that the permeability surrounding the fracture path is higher than the matrix permeability (1e-16 m 2 ), hence allowing faster pore pressure diffusion in this zone. The time at which the hydraulic fracture intersects F1 scales with the permeability, since higher permeability values promote faster diffusion and therefore a slower fracture propagation [START_REF] Carrier | Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model[END_REF]. Meanwhile, due to fluid leak-off from the hydraulic fracture into the matrix, the pressure front propagates faster than the hydraulic fracturing front (Fig. 4.9). Here, the pressure front is defined as the pressure contour above the initial pore pressure. We noticed that when the matrix permeability is relatively low (case A2 and A8), the difference in the timing between the two propagating fronts is small (Fig. 4.9a and Fig. 4.9b). When the permeability is high (case A12), the pressure front travels significantly ahead of the fracturing front (Fig. 4.9c). If the leak-off is ignored (e.g., the KGD problem described in the Appendix A), the pressure front and the fracturing front coincide Therefore, the fault F1 will be overpressurized before a hydraulic connection with the injection source is reached. This over-pressure combined with poro-elastic stressing results in a reduction of the effective normal stress and induces slip along F1 because the induced shear stress exceeds its shear strength. The time lag of the over-pressurized front and fracturing front for cases A2 and A8 is very small (< 0.1 s) such that the sliding looks more continuous. However, for case A12 with a higher surrounding permeability, this time lag becomes obvious (∼ 12 s), and F1 undergoes two distinct slips due to over-pressurization and hydraulic shearing, respectively. During the first slip, the total shear slip and slip rate are smaller than those obtained for case A2 and A8, since the pressure perturbation is lower. During the second slip, there is mainly aseismic slip and this is because most of the accumulated strain energy has been released during the first slip. As a result, the slip and slip rate of F1 exhibit an obvious downward trend as permeability increases (Fig. 4.10a).

For each tested value for the permeability, we observed that F2 is activated by the pressure diffusion and poro-elasticity though the activation time differs for each case. A higher permeability requires less fluid injection, i.e., 150 -160 s for k f 1 = 1e-15 m 2 , 120 -130 s for k f 2 = 5e-15 m 2 , and only 100 -110 s for k f 3 = 1e-14 m 2 . When the injection time is identical, the higher the permeability, the shorter the time lag for F2 reactivation. For example, with an injection duration of 160 s, F2 has a quiescence of 360.1 s and 34.7 s after shut-in for case A2 and case A8, whereas it is even activated during the injection for case A12. The permeability has a limited influence on the magnitude of the

Friction coefficient

The effect of the friction coefficient on fault slip behavior due to direct injection can be seen from Fig. 4.11(a), corresponding to cases A2, B1 and B2 in Tab. 4.3. Triggering by direct fluid injection, faults with a higher friction coefficient tend to display a higher level of slip and higher slip rates. As stated by [START_REF] Ngo | Modeling of fault slip during hydraulic stimulation in a naturally fractured medium[END_REF], this is because more accumulated energy are required to activate faults with higher friction coefficient, and accordingly more energy is released during slip. In addition, unlike faults with a smaller friction coefficient, which are more likely to be activated by pressure diffusion and poro-elastic stressing before hydraulic shearing, the faults with higher friction coefficient tend to slip only during hydraulic shearing and therefore more violently.

As for the shut-in stage, because F2 cannot be activated with permeability k f = 1e-15 m 2 in case B1, we keep the friction coefficient of F1 constant as µ 1 = 0.2 and set case A12 as the base case (k f = 1e-14 m 2 ). When the friction coefficient of F2 is 0.4, F2 activation can be observed either by increasing injection duration to 180 s (case B4) or permeability to 1e-13 m 2 (case B5). For µ 2 = 0.6, F2 cannot be triggered unless these two parameters are simultaneously increased (case B7). These illustrate that a higher friction coefficient hinders fault reactivation. However, once F2 slips, higher slip and slip rate are observed when the friction coefficient is larger (Fig. 4.11(b)), similar to

Fault orientation

On the basis of case A2, we conduct additional simulations where we varied the fault orientation, by changing the fault angle α (tested value are 45°, 20°and 10°, see Tab. 4.4).

During the fluid injection stage, the hydraulic connection of the injection well-bore and the fault F1 occurs at 40.68 s, 39.48 s, 38.12 s, and 35.83 s when the fault angle increases from 10°t o 45°, suggesting that a smaller fault dip tends to hinder hydraulic fracture propagation. We infer In a second stage, we keep the injection rate as constant (0.001 m 3 /s), but we change the injection protocol. We consider a step-wise injection schedule (a period of constant injection followed by a pause) by maintaining the same total volume of fluid injected (Fig. 4.13b).

In case E4, the stimulation is divided evenly into 4 injections and 3 pauses (i.e., a cycling stimulation), each injection/pulse has a duration of 40 s. Since F1 reactivation occurs at 38.12 s < 40 s in case A2, this setting will not change the response of F1 during injection. However, the relaxation of stimulation during the injection favours pressure relaxation, which prevents the reactivation of F2. In case E5, the first injection and suspending time in case E4 are changed to 30 s and 50 s, respectively. In this case, F1 is activated during the first pause (47.83 s) due to pressure diffusion and poro-elastic stressing, while mainly hydraulic shearing occurs during the second injection (80.27 s). Compared to case E4, not only is F2 reactivation after shut-in prevented in case E5, but the injection-induced slip and slip rate of F1 are significantly reduced owing to the double slip effect, i.e., 1.36 mm and 50.4 mm/s for case E4 vs. 0.75/0.18 mm and 26.97/25.95 mm/s for case E5. It is equivalent to artificially extend the interval between the pressure diffusion and hydraulic shearing, producing a longer and more stable energy release.

Discussion

Comparison with earthquake fault scaling relations

It is well known that unstable fault slip is accompanied by a cluster of induced seismic events [START_REF] Zoback | Reservoir geomechanics[END_REF][START_REF] Ellsworth | Injection-induced earthquakes[END_REF]. To quantify their strength, [START_REF] Aki | Generation and propagation of g waves from the niigata earthquake of june 16, 1964.: Part 2. estimation of earthquake moment, released energy, and stress-strain drop from the g wave spectrum[END_REF] first defined the scalar seismic moment M 0 as:

M 0 = GdA (4.15)
where G is the shear modulus, d and A denote the fault shear displacement and fault rupture area, respectively.

The moment magnitude M w can be then computed by assuming a GR(Gutenberg-Richter) statistics as [START_REF] Hanks | A moment magnitude scale[END_REF]:

M w = 2 3 [log 10 (M 0 ) -9.1] (4.16)
where we consider a GR b-value equal to 1 in deriving Eq. ( 4.16).

We calculate the moment magnitude for each case based on the obtained seismic slip. Considering that the fault slip occurs during a very short time, we assume that the rupture covers the whole fault, therefore the rupture length is identical to the fault length, and we can approximate the sliding distance as the average seismic slip computed [START_REF] Lei | Modelling fluid injection-induced fracture activation, damage growth, seismicity occurrence and connectivity change in naturally fractured rocks[END_REF]. For calculating the rupture area, we assume that the rupture width is equal to the rupture length, which is supported by a large number of reported earthquakes [START_REF] Yoon | Discrete bonded particle modelling of fault activation near a nuclear waste repository site and comparison to static rupture earthquake scaling laws[END_REF].

We compare the results with two earthquake fault scaling relations. The first one describes the relationship between induced magnitude, fault size, shear slip and stress drop [START_REF] Zoback | Earthquake triggering and large-scale geologic storage of carbon dioxide[END_REF]. In our simulations we compute a shear slip between 0.22 and 1.57 mm, corresponding to Mw from -0.39 to 0.18. The related stress drop ranges from 4.29 MPa to 19.79 MPa, suggesting a good consistency with adopted scaling laws (Fig. 4.15). The second comparison relates the magnitude with the rupture area. Although with a smaller rupture size, our results are still in good agreement with numerous earthquakes (Fig. 4.16).

Comparison of induced seismicity during injection and after shut-in

Our results show that pre-existing faults can be reactivated and exhibit dynamic slip under certain conditions even when faults and hydraulic fractures are non-hydraulically connected, illustrating that seismic events can be induced during the post-injection period triggered by pore pressure diffusion and poro-elastic stress change mechanisms. This is in agreement with field observations as those reported in Soultz-sous-Forêt EGS site [START_REF] Baisch | A numerical model for fluid injection induced seismicity at soultz-sous-forêts[END_REF]. In addition, this post-injection reactivation of fault occurs after a period of quiescence, sometimes even exceeding twice of the injection time in our simulations (e.g., case A2 and case E2). This delayed phenomenon has been also observed in other EGS projects as in Rittershoffen [START_REF] Lengliné | Seismicity related to the hydraulic stimulation of grt1, rittershoffen, france[END_REF].

For most of the cases that display reactivation of F2 after shut-in we note that induced slip and rates computed at F1 during injection are larger than those obtained in F2 during the post injection stage. However, there are also some cases showing an inverse trend, i.e., F2 produces more shear [START_REF] Yoon | Discrete bonded particle modelling of fault activation near a nuclear waste repository site and comparison to static rupture earthquake scaling laws[END_REF] slip and higher slip rate than F1. We single out these cases and find that they generally share similar conditions, i.e., the diffusivity (permeability) in their nearby regions is extremely high (1e-14 m 2 or 1e-13 m 2 ), e.g., cases A10 -A12, B3 -B5. In these cases, F1 reactivation due to large pressure gradient occurs earlier than the direct hydraulic shearing, but this pressure disturbance is smaller compared to F2 after shut-in. Thus, the strain energy of F1 is released more evenly than F2. The results may provide an explanation for a common observation in the field, which is that the post injection might induce larger magnitude earthquakes [START_REF] Häring | Characterisation of the basel 1 enhanced geothermal system[END_REF][START_REF] Mukuhira | Characteristics of large-magnitude microseismic events recorded during and after stimulation of a geothermal reservoir at basel, switzerland[END_REF]. They also shed light on mitigating the risk of fluid induced earthquake. As shown in case E5, which displays double slip events in F1 during active injection, a lower surrounding permeability (1e-15 m 2 ) significantly reduces the induced seismic slip and its slip rates. By injecting the same amount of fluid, although cases E5 and E3 both lower the slip of F1 and evade the post-activation of F2, case E5 is more efficient because it consumes less time (280 s vs. 320 s). Therefore, on the basis of our analysis, we can conclude that a soft cycling stimulation strategy might be effective in controlling fluid injection induced seismicity, as suggested by [START_REF] Hofmann | Cyclic soft stimulation (css): a new fluid injection protocol and traffic light system to mitigate seismic risks of hydraulic stimulation treatments[END_REF][START_REF] Zang | How to reduce fluid-injection-induced seismicity[END_REF][START_REF] Zhuang | Cyclic hydraulic fracturing of pocheon granite cores and its impact on breakdown pressure, acoustic emission amplitudes and injectivity[END_REF].

Our results also suggest that the permeability beyond the stimulated zone is of great importance for the shut-in phase. After shut-in, high pressure in the stimulated region diffuses to the surrounding porous matrix, and coupled to induced poro-elastic stress changes can lead to activation of F2 and induced seismicity, thereby extending the seismic region. The work of [START_REF] Mcclure | Generation of large postinjection-induced seismic events by backflow from dead-end faults and fractures[END_REF] showed that the fluid pressure is redistributed and pressure front progresses further after shut-in, extending the pressurized area and induces seismicity. [START_REF] Mukuhira | Pore pressure behavior at the shutin phase and causality of large induced seismicity at basel, switzerland[END_REF] studied the behavior of pore pressure migration at shut-in phase in Basel and pointed out that large events tend to occur in regions with enhanced permeability. In our case, we set large permeability outside the stimulated region, causing pore pressure to migrate far field the seismic zone after termination.

Role of initial stress state

Our parametric study shows that there are multiple factors affecting fault slip behavior both during injection and after shut-in. For example, higher permeability reduces the seismic slip but more via double slip effect during the injection period and has little impact on the induced fault slip during the shut-in stage (Fig. 4.10), and the soft cycle injection lowers the seismic slip rate effectively (Fig. 4.14). In contrast, the fault orientation seems to have the greatest influence on the fault slip. As shown in Fig. 4.12, when the fault angle decreases, the seismic slip and slip rate drop significantly and the faults tend to slip aseismically regardless of during injection or after shut-in. This may be related to the initial stress state of the faults. To evaluate the initial potential for fault reactivation, we calculate the slip tendency ST, which is defined as the ratio between the absolute shear stress and effective normal stress [START_REF] Blöcher | Evaluating micro-seismic events triggered by reservoir operations at the geothermal site of groß schönebeck (germany)[END_REF]:

ST = |τ | |σ ′ n |
(4.17)

The initial shear stress and effective normal stress of fault are given by [START_REF] Jaeger | Effect of fluid pressure heterogeneity on injection-induced fracture activation[END_REF]:

τ = σ 1 -σ 3 2 sin(2α) (4.18) σ ′ n = σ ′ 1 + σ ′ 3 2 + σ ′ 1 -σ ′ 3 2 cos(2α) (4.19)
At the given initial stress state (Tab. 4.1), the slip tendency ST = 0.034, 0.069, 0.096, 0.108 for the fault angle α = 10°, 20°, 30°and 45°, suggesting an increasing critical condition for fault reactivation. Moreover, a larger initial slip tendency leads to higher seismic slip and slip rate under the same perturbation (see Fig. 4.12). To further validate this hypothesis, we conduct another simulation by only reducing the minimum principal stress σ 3 from 29 MPa to 24 MPa from case A2. Accordingly, the initial slip tendency increases from a nominal value of 0.096 to 0.155, approaching further the static friction coefficient (0.2). The results show that with only 8.05 s and 28.2 s fluid injection, F1 and F2 are both activated by coupled fluid pressure diffusion and poro-elastic stressing. This is consistent with the main conclusions derived in the study by [START_REF] Lee | Managing injection-induced seismic risks[END_REF], where they analyzed induced seismicity in the Pohang EGS project and found that the fault which is responsible for the earthquake, was critically stressed, susceptible to slip and very sensitive to even small perturbations. Furthermore, the slip rates of F1 and F2 also increase (mainly due to poro-elastic stressing since F1 slips aseismically when it connects the injection point at 16.87 s), indicating much stronger energy releases. This has important implications for controlling fluid injection-induced seismicity, in particular for post-injection induced seismicity targeting far-field faults. As suggested by [START_REF] Baisch | A numerical model for fluid injection induced seismicity at soultz-sous-forêts[END_REF], post-event magnitude is predominantly affected by the fault geometry and fluid pressure elevation brings stress conditions close to critical values. [START_REF] Mukuhira | Characteristics of large-magnitude microseismic events recorded during and after stimulation of a geothermal reservoir at basel, switzerland[END_REF] and [START_REF] Mukuhira | Pore pressure behavior at the shutin phase and causality of large induced seismicity at basel, switzerland[END_REF] also demonstrated that even a minor increase in pore pressure was sufficient to initiate shear slip and trigger large events at large differential stress. Indeed, many researches even stated that the initial stress controls the nature of the fluid induced seismicity [START_REF] Yoon | Discrete bonded particle modelling of fault activation near a nuclear waste repository site and comparison to static rupture earthquake scaling laws[END_REF][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] Passelègue | Initial effective stress controls the nature of earthquakes[END_REF]. This may also be one aspect to get an insight into an observation that a field area of little natural seismicity may have a lower propensity to produce felt or damaging earthquakes [START_REF] Evans | A survey of the induced seismic responses to fluid injection in geothermal and co2 reservoirs in europe[END_REF]. Therefore, to mitigate the unwanted seismicity, a good prior knowledge of the crust before treatment is required.

Conclusion

In this study, we conducted a coupled hydraulic-mechanical numerical model to investigate the induced fault slip behavior during injection and after shut-in. It shows the capability of generating seismic slip either owing to direct hydraulic shearing or fluid pressure diffusion coupled with poroelastic stressing when there is no connection between the hydraulic fracture and the pre-existing fault.

Through parametric study, we found that fault seismic slip can be attributed to the combined effect of multiple factors but with different emphases. Permeability along the fracture path is of great importance not only for the after shut-in stage, but also for the fracturing process. It allows fluid pressure diffusion induced poro-elastic stress change to activate the fault before the hydraulic shearing. This averages the strain energy release and lowers the seismicity magnitude during injection, providing an explanation for the occurrence of larger earthquakes after shut-in. It also highlights the effectiveness of soft cyclic injection since the fault strain energy can be partly released during each relaxation. Moreover, although the larger friction coefficient favors fault slip and the soft in-jection can alleviate the slip, the impact of fault orientation on the slip is more dominant. Finally, the combined effect of these factors can be largely attributed to the initial slip tendency and hence the initial stress state. We confirmed that critically-stressed faults are more likely to produce larger shear slip even with only small perturbations. Our results show good consistency with some field observations and previous studies, which have important implications for mitigating fluid injection induced seismicity. In the next step, the influence of fracture networks will be considered.

Acknowledgment

This work was supported by the China Scholarship Council (CSC) [grant number 201808510128].

Appendix. Validation of the hydraulic fracture propagation

For any complex numerical models involving hydraulic fracture simulation, [START_REF] Lecampion | Numerical methods for hydraulic fracture propagation: a review of recent trends[END_REF] emphasized the necessity of verification against the analytical solutions. Here the KGD (Kristianovic-Geerstma-de Klerk) model for plane strain hydraulic fracture is selected for validation [START_REF] Garagash | Plane-strain propagation of a fluid-driven fracture during injection and shutin: Asymptotics of large toughness[END_REF].

The scaling solutions for the hydraulic aperture w, the average net pressure p, and the fracture length l are given by [START_REF] Detournay | Propagation regimes of fluid-driven fractures in impermeable rocks[END_REF]: where x ′ = x/l(t) is the scaled coordinate, 0 ≤ x ′ ≤ 1. ε(t), L(t), P (t) denote respectively a small dimensionless parameter, a length scale if the same order as the fracture length l, and a dimensionless evolution parameter.

         w = ε(t) • L(t)Ω[x ′ , P (t)] p = ε(t) • E ′ Π[x ′ , P (t)] l(t) = L(t) • γ[P (t)]
The hydraulic fracture propagation is controlled by either a toughness-dominated regime (κregime) or a viscosity-regime (M -regime) [START_REF] Carrier | Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model[END_REF]. The dimensionless viscosity M and toughness κ are given by:

M = κ -4 = µ ′ E ′ 3 Q 0 K ′ 4 (4.21)
where µ ′ = 12µ;

E ′ = E 1 -v 2 ; K ′ = 4 2 π 1/2 K IC .
The toughness scaling is most appropriate for cases when the viscosity scaling M is small. The asymptotic solution F k (Ω k , Π k , γ k ) is given by:

F k = F k 0 + M F k 1 (4.22)
where F k 0 (Ω k 0 , Π k 0 , γ k 0 ) is the zero-viscosity solution given by: 4.24) where

               Ω k 0 = π 1/3 2 (1 -x ′2 ) 1/2 Π k 0 = π 1/3 8 γ k 0 = 2 π 2/3
F k 1 (Ω k 1 , Π k 1 , γ k 1 ) is:            Π k 1 = Π * k 1 + ∆Π k 1 Ω k 1 = 4Π * k 1 (1 -x ′2 ) 1/2 + 2 π 1 0 f (x ′ , x)∆Π k 1 (x) dx γ k 1 ≃ -2.72 (
∆Π k 1 (x) = 4 3π 2/3 ln(1 -x ′2 ) - 2x ′ arccos x ′ π 2/3 (1 -x ′2 ) 1/2 (4.25) f (x ′ , x) = - 4 π ln 1 -x ′2 + √ 1 -x 2 1 -x ′2 - √ 1 -x 2 (4.26)
For the toughness scaling,

         ε k = K ′4 E ′4 Q 0 t 1/3 L k = E ′ Q 0 t K ′ 2/3
(4.27) Fig. 4.17 shows the geometry and mesh used in the simulation. Here the 2D model has a dimension 45 m × 60 m, with a predefined hydraulic fracture path at the middle along the x-axis, simulated by cohesive elements with pore pressure degree of freedom. The well-bore is simplified Figure 4.17: Geometry and mesh of the KGD model to be a point and the point injection rate Q is assumed to be constant. Note that only half of the space is simulated due to geometric symmetry.

For cohesive elements, a coarse mesh density may reduce the accuracy [START_REF] Zielonka | Development and validation of fullycoupled hydraulic fracturing simulation capabilities[END_REF]. Here we use a finer mesh (the smallest mesh size around the hydraulic fracture path is 0.5 cm) to circumvent this effect [START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF][START_REF] Ngo | Modeling of fault slip during hydraulic stimulation in a naturally fractured medium[END_REF]. All material parameters used in the simulation are given in Tab. 4.6. According to Eq. ( 4.21), we obtain M = 0.0142, indicating the toughness-dominated regime.

In this simulation, we perform 20 s of hydraulic stimulation. The fracture length, fracture aperture and pressure at the injection point, and the corresponding relative errors are given in Fig. 4.18. Results show a good agreement between our numerical results and the analytical solutions except that the numerically obtained injection pressure has an increasing error as time increases. This is due to poro-elastic coupling effect [START_REF] Carrier | Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model[END_REF], which is not considered in the analytical model. 1) Steady-state flow solved by the Navier-stokes equation over a 3D self-affine rough fracture and transient flow solved by the pressure diffusion equation over the same fracture represented by a surface with aperture distribution as a parameter provide similar results. The effective hydraulic diffusivity of the rough fracture quantified by comparing it with the transient pressure distribution of the parallel plate model deviates from the results calculated using the hydraulic aperture. This effect is more pronounced with fracture closure due to the increase in roughness and tortuous channel flow. The hydraulic diffusivities show enhancement or inhibition in different flow directions and are comparable to field observed values when approaching the flow percolation threshold.

2) A 3D finite element model with a rough fracture embedded in the porous matrix shows the ability to determine the permeability and fracture normal stiffness simultaneously. Results obtained from 6 independent self-affine fractures demonstrate that the evolution of the permeability and stiffness during initial fracture mechanical closure (below 20% of the contact area) and subsequent minerals sealing (below the percolation thresholds) exhibits similar behaviors. Derived from a physical point of view, permeability and fault volume can be approximated by a quadratic equation, whereas stiffness and fault volume follow an exponential equation. These two relationships were validated by numerical simulations and the two quantities thus can be related by a generic law. In addition, the permeability at the percolation threshold shows anisotropy with variations in 7-8 orders of magnitude. This behavior is of great importance for chemical stimulation aimed at improving the permeability by an injection of fluids that contain chemical compounds (e.g., acids). Finally, assessing the effectiveness of chemical treatments by using the permeability-stiffness relationship would be a potential solution since the stiffness might be accessible (e.g., by geophysical probing).

3) A coupled hydro-mechanical numerical model using finite elements for porous rocks and cohesive elements for fractures was successfully developed to simulate the hydraulic stimulation of crystalline rocks . Both hydraulic fracturing and hydraulic shearing mechanisms were incorporated to reproduce the injection in the area around the wellbore, including several physical processes: porous rock deformation, pore fluid flow, fracture fluid flow, fracture initiation and propagation and fracture slip. Some important results were obtained from parametric study as follows. Pressure diffusion in the porous matrix induces poroelastic stress change that can reactivate pre-existing fracture prior to fluid inflow (i.e., before the intersection of hydraulic fracture and pre-existing fracture) but depending on the matrix permeability. This pre-reactivation reduces the fracture slip causing by the direct fluid inflow and thereby lowers the maximum magnitude of seismicity, which provides another insight to interpret the mechanism of soft cyclic stimulation. Fracture properties such as its orientation and friction coefficient impact the seismic slip behavior, but the initial fracture slip tendency seems has more profound influence because only small perturbations can reactivate the critically stressed fracture. These are consistent with other studies and observations.

Perspectives for future research

As stated in the first chapter of the thesis, this work attempts to improve the understanding of the hydraulic and mechanical properties of a single fracture by using numerical approaches. It presents some results under certain conditions, for example, the influence of roughness on the fracture flow was studied only at small pressure perturbations, which limits the application in real projects. Therefore, there are still a lot of questions remains unclear or unanswered. From our perspectives, some work with highly interest and potential are listed as follows.

1) Modeling fluid flow at larger pressure perturbations (e.g., fluid injection). In this work, fluid flow simulations on a single rough fracture were conducted at low pressure gradient to ensure laminar flow. However, at injection conditions, the fracture flow and deformation behaviors might differ. High flow rates are likely to change the normal stress and alter the fracture volume and morphology, leading to nonlinear fluid flow [START_REF] Murphy | Semianalytical solutions for fluid flow in rock joints with pressure-dependent openings[END_REF]. Although in the experiments of fluid injection into a single rough granodiorite fracture, [START_REF] Vogler | Experiments and simulations of fully hydro-mechanically coupled response of rough fractures exposed to highpressure fluid injection[END_REF] did not find significant impact of fluid injection pressures on fracture opening, the injection pressures deviate from numerical simulations predicted by the local cubic law.

2) Considering mineral precipitation and dissolution processes inside the fracture. In this work, to simplify the numerical procedure, minerals growth is assumed to be independent of time and position, and simulated by adding homogeneous layers. However, previous studies showed that mineral precipitation/dissolution processes are related to various factors, such as minerals type [START_REF] Li | Scale dependence of mineral dissolution rates within single pores and fractures[END_REF], flow rates [START_REF] Griffiths | Geothermal implications for fracture-filling hydrothermal precipitation[END_REF], temperature or boiling [START_REF] Dobson | Experimental and numerical simulation of dissolution and precipitation: implications for fracture sealing at yucca mountain, nevada[END_REF][START_REF] Polak | Permeability reduction of a natural fracture under net dissolution by hydrothermal fluids[END_REF], fluid saturation [START_REF] Detwiler | Permeability alteration due to mineral dissolution in partially saturated fractures[END_REF], pressure shadows, mineral availability, etc. In addition, fracture deformation are affected by coupled dissolution and mechanical stress, suggesting the importance of rock-water reactions [START_REF] Detwiler | Experimental observations of deformation caused by mineral dissolution in variable-aperture fractures[END_REF][START_REF] Ellis | Dissolutiondriven permeability reduction of a fractured carbonate caprock[END_REF]. Overall, the study of the chemical reactions within the fracture should not be separated from the hydro-mechanical model. It is essential to consider a coupled T-H-M-C model to better reproduce the entire process in subsurface.

3) Switching the research object from a single fracture to fracture networks. Although a single fracture is an essential element worthy of study, it does not fully reflect subsurface flow patterns influenced by many fractures, even if it serves as the main flow conduit (e.g., a large fault), since the surrounding zone (i.e., the damage zone) often contains thousands of fractures. In this regard, we have two proposals that are relevant to this thesis.

First, to consider the influence of fracture networks on the hydraulic stimulation around the injection wellbore. In this case, fracture distribution must affect the fracture propagation modes, the intersections between artificial fractures and natural fractures, fracture slip behavior, etc. Using such model to optimize the injection strategies under different stress conditions would be of interest.

Second, to measure the hydraulic diffusivity of a three dimensional fracture network by considering fracture roughness on each individual fracture. For example, De [START_REF] De Dreuzy | Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks (dfn)[END_REF] pointed out that the permeability reduction can reach a factor of 6 when using varied aperture instead of constant mean aperture in 3D discrete fracture networks (DFNs) where self-affine aperture distribution is spatially correlated and follows a truncated Gaussian law. [START_REF] Haagenson | Seismic diffusivity and the influence of heterogeneity on injection-induced seismicity[END_REF] corrected the effective hydraulic diffusivity of the heterogeneous rock formation to seismic diffusivity by using DFNs in the matrix (DFNM) under fluid injection conditions and successfully applied it to a case study of Soultz geothermal project. It is then promising to incorporate roughness into each fracture of such DFNM model and to see further effects.

To use the findings mentioned above for simulating the physical processes on a reservoir scale, a 2D representation of the fracture/fault is recommended. This was performed for the hydraulic response (Chapter 2) but is under development for the thermal, chemical and mechanical processes.

To have an optimized simulation the thermal, chemical and mechanical alteration must be represented by 2D fault surfaces embedded in a 3D porous matrix. This numerical approach is under development and it will be up-scaled to the real EGS cases, e.g., to investigate the induced fault slip caused by the fluid injection of the Geoven Deep Geothermal Project [START_REF] Schmittbuhl | Induced and triggered seismicity below the city of strasbourg, france from november 2019 to january 2021[END_REF], which was stopped due to an earthquake of magnitude ML3.6 induced below the city of Strasbourg on December 2020. Related work is in progress.

1. 7

 7 Sketch of cohesive zone model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Coupled deformation/pressure cohesive element with 6 nodes. (a) undeformed configuration; (b) deformed configuration. Red dots at corners denote displacement/pressure degrees of freedom, while blue dots at the middle-edge indicate only pressure degrees of freedom. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Intersection of a horizontal fracture and a vertical fracture. Four cohesive elements 1, 2, 3, 4 share the same node 10 at the intersection. . . . . . . . . . . . . . . . . . 2.1 Sketch of the whole simulation scheme: fracture surface/aperture generation; fracture closure; use of the Navier-Stokes equations for steady state flow through a volume-represented fracture and the pressure diffusion equation for transient flow along a surface-represented fracture; and determination of the effective hydraulic diffusivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XVIII 2.2 Statistical invariance illustrated by magnifying portions of the profiles from the Corona Heights fault surface in (a) parallel (H ∥ = 0.65 ± 0.04) and (b) perpen-

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.4 Histogram of the height distribution of our generated isotropic self-affine surface at the 512 m scale, with Hurst exponent H = 0.8 and roughness amplitude h rms = 0.1 m, consistent with the field observations of Candela et al. (2012) (see blue circles in Fig. 1). The largest height fluctuations are on the order of 3h rms ≈ ±0.3 m for a lateral extension of 512 m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.5 2D cross-section along the (x, z)-plane and (y, z)-plane (with a vertical magnification 500 times on the z-axis scale with respect to the x, y-axis) of the self-affine fracture aperture h(x, y) at two steps of the fracture closing. (a) Fully open fracture without contact. z = z 1 and z = z 2 represent the mean plane of the top surface s 1 (x, y) and the bottom surface s 2 (x, y), respectively. d m = z 1 -z 2 denotes the mean aperture. (b) Partial open fracture after the contact of some asperities. The shaded area indicates contact regions where the apertures are set to zero; d m here denotes the mean separation of the open area (light blue zones) . . . . . . . . . . . 35 2.6 Schematic diagram of the two representations of a rough fracture with L = 64 m. (a) Volume-represented fracture with 52,128 prism elements (8 layers along z); (b) surface-represented fracture with 4,096 quadrilateral elements. Both representations share the same mesh size of 1 m × 1 m along x and y . . . . . . . . . . . . . 36 2.7 Schematic of the initial conditions and boundary conditions applied in the x-and y-directions for each rough fracture. (a) Flow along the x-axis; (b) flow along the y-axis. The inlet and outlet boundaries are marked as red and blue, respectively . . 40 2.8 Comparison of the fluid flow between volume-and surface-represented rough fractures. The aperture evolution is given by the ratio of the hydraulic aperture d h to the mean aperture d m as a function of the h rms -normalized mean aperture d m /h rms . The horizontal dashed line denotes the parallel plate model in which d h /d m = 1, while the solid vertical line indicates the divide where the fracture fully opens (right part) and comes into contact (left part) . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 XIX 2.9 Comparison between the simulation results for the parallel plate model (red squares) and analytical solutions (solid lines) of pressure diffusion at t 1 = 5 × 10 -6 s, t 2 = 1.5 × 10 -5 s and t 3 = 4.5 × 10 -5 s . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Map of the aperture distribution with H = 0.8, h rms = 0.1 m and d m = 0.3 m at the 512 m scale. The fluid flows from left to right in the x-direction and from bottom to top in the y-direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Schematic diagram of the pressure field evolution in a rough fracture and the parallel plate model with the same size (512 m) and the same mean aperture (0.3 m) at t 1 = 5 × 10 -7 s, t 2 = 2 × 10 -6 s and t 3 = 5 × 10 -6 s. (a)-(c) Fluid flow along the x-axis for the parallel plate model; (d)-(f) fluid flow along the x-axis for the rough fracture; and (g)-(i) fluid flow along the y-axis for the rough fracture . . . . . . . . 2.12 Schematic diagram of the local flux evolution in a rough fracture and the parallel plate model with the same size (512 m) and the same mean aperture (0.3 m) at t 1 = 5 × 10 -7 s, t 2 = 2 × 10 -6 s and t 3 = 5 × 10 -6 s. (a)-(c) Fluid flow along the x-axis for the parallel plate model; (d)-(f) fluid flow along the x-axis for the rough fracture; and (g)-(i) fluid flow along the y-axis for the rough fracture . . . . . . . . 2.13 Pressure diffusion along the x-axis (averaged along y) with d m /h rms = 3 at t 1 = 5 × 10 -7 s, t 2 = 2 × 10 -6 s and t 3 = 5 × 10 -6 s. (a) Comparison between the rough fracture and the parallel plate model with the same d m = 0.3 m; (b) comparison between the rough fracture and the parallel plate model with the best fitting hydraulic diffusivity D e = 0.423D m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.14 Pressure diffusion along the y-axis (averaged along x) with d m /h rms = 3 at t 1 = 5 × 10 -7 s, t 2 = 2 × 10 -6 s and t 3 = 5 × 10 -6 s. (a) Comparison between the rough fracture and the parallel plate model with the same d m = 0.3 m; (b) comparison between the rough fracture and the parallel plate model with the best fitting hydraulic diffusivity D e = 1.691D m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 (left) Evolution of the contact area as a function of the normal closure (defined as the fracture's vertical displacement, which is positive along negative z direction) and mean aperture d m . (right) Evolution of the hydraulic aperture d h as a function of the mean aperture d m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.16 Log-log graph of the effective hydraulic diffusivity D e as a function of the hydraulic aperture d h . The dotted line corresponds to the parallel plate model (Eq. 2.19) . . . XX LIST OF FIGURES 2.17 Schematic diagram of the pressure diffusion evolution along the x-axis of the rough fracture with different contact areas. (a)-(c) t 1 = 2 × 10 -6 s, t 2 = 4 × 10 -6 s and t 3 = 1 × 10 -5 s with 5.43% contact; (d)-(f) t 1 = 2 × 10 -4 s, t 2 = 4 × 10 -4 s and t 3 = 0.001 s with 33.31% contact; and (g)-(i) t 1 = 0.002 s, t 2 = 0.004 s and t 3 = 0.01 s with 50.51% contact . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.18 Schematic diagram of the local flux evolution along the x-axis of the rough fracture with different contact areas. (a)-(c) t 1 = 2 × 10 -6 s, t 2 = 4 × 10 -6 s and t 3 = 1 × 10 -5 s with 5.43% contact; (d)-(f) t 1 = 2 × 10 -4 s, t 2 = 4 × 10 -4 s and t 3 = 0.001 s with 33.31% contact; and (g)-(i) t 1 = 0.002 s, t 2 = 0.004 s and t 3 = 0.01 s with 50.51% contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.19 Ratio of the effective hydraulic diffusivity and the parallel plate fracture diffusivity D e /D m as a function of the h rms -normalized mean aperture d m /h rms . The vertical black line corresponds to the first asperity contacts when closing the fracture at d m /h rms ≈ 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.20 Two examples of the evolution of the relative hydraulic diffusivity D e /D m for two other aperture fields with the same Hurst exponent H = 0.8 and the same RMS h rms = 0.1 m when changing the seed of the self-affine surface generator . . . . . . 2.21 Pore pressure averaged along the y-axis as a function of t/t * at x = L/2 in the case of the parallel plate model with d m = 0.3 m . . . . . . . . . . . . . . . . . . . . . 2.22 Least square error as a function of D e /D m . (a) D e1 = 0.4D m with the step of the test diffusivity ∆D = 0.1D m ; (b) D e2 = 0.42D m with ∆D = 0.01D m ; and (c) D e3 = 0.423D m with ∆D = 0.001D m . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Sketch of the model: a partially open fault (colored horizontal cross-section) is embedded in a 3D block (in dark gray) of size 256 × 256 × 256 m 3 loaded by a normal stress σ n . (a) The contact areas (in light gray) along the fracture support the applied normal stress with variable local amplification of the stress (reddish zones). The open part of the fracture where fluid circulates is shown in blue. (b) 3D view of the fluid circulation in the fracture (colored horizontal cross-section) and in the pore space of the matrix (white arrows) when a pressure drop is applied in the x-direction (the vertical thick, blue layers indicate the pressure field in the matrix). The contact areas are transparent here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Generation of the fault topology and aperture map: (a) aperture map with only a few contact points; (b) power density function (pdf) of the aperture distribution with a Gaussian fit (average aperture is h = 0.27 m; RMS is h rms = 0.09 m). . . . . . . . XXI 3.3 Cross-section through the block along the x-direction with a large vertical magnification to show the fracture closure process due to two different mechanisms. (a) Initial unloaded fracture aperture: fracture surfaces are just touching one another at a few contact points. (b) Phase 1 -mechanical closure (contact area < 20%): the contact area increases because of a normal loading (imposed normal displacement) applied on the top of the host rock; (c) asperities that overlap following a 'rigidplastic' behavior. (d) Phase 2 -sealing closure due to mineral deposits (contact area > 20%). Successive layers of sealing material are assumed to grow homogeneously with a constant thickness from the border to the center of the open space (syntaxial mode). The open space for fluid flow is shown in light blue, whereas the host rock is shown in light gray. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Definition of the fracture displacement u f . u f = u b -u m , where u b is the displacement of the fracture-matrix system, while u m is the displacement of the intact rock matrix under the same applied normal loading σ n . . . . . . . . . . . . . . . . . . . 3.5 Evolution of the normal stress over the fracture closure (imposed displacement). The slopes of the red lines at each step of the mechanical/sealing closure correspond to its normal stiffness. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Evolution of the relative contact area A c /A as a function of the normalized normal displacement imposed on the fractured block (u f -h)/h rms , where the mean aperture is h and the RMS is h rms . The pink domain corresponds to the mechanical closure of the fracture, assuming an overlapping model of contact (i.e., purely plastic rheology of the asperities). The light blue domain corresponds to the sealing closure due to homogeneous mineral deposits in the open fracture space. The evolution of the contact area can be fitted by the cumulative distribution function of a normal distribution y = (1 + erf (x/ √ 2))/2. The transition from a mechanical to a sealing closure is defined at 20% of contact, which corresponds to (u f -h)/h rms ∼ -1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Cross-section view of the initial conditions and boundary conditions for (a) the rock matrix and (b) the fracture. Transfer functions at the fracture-rock interface are also illustrated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 (a) Finite elements of the open fracture space (in green) and of the rock mass (in brown). Coarse-graining of the elements develops when moving away from the fracture. (b) Details of the meshing of the open space of the fracture by prism elements.

  3.12 Map view of the flow inside fault F1. (a) Pressure drop is applied along the xdirection. (b) Pressure drop is applied along the y-direction. The white zones indicate the contact regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Evolution of the permeability as a function of the percolation fault volume (the fault volume minus the fault volume at the percolation threshold (V f -V perc f

3. 14

 14 Fracture stiffness as a function of the square root of permeability. The pink line denotes the prediction from Eq. (3.19) with a representative fault volume of 3000 m 3 at the percolation threshold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Evolution of the permeability ratio (k y /k x ) as a function of the fault volume during the mechanical & sealing closure process for the three studied faults. . . . . . . . . 3.16 Scaling of fracture normal stiffness as a function of fracture length; modified from Hobday and Worthington (2012) and Worthington and Lubbe (2007). . . . . . . . . 3.17 Fracture stiffness characteristic as a function of fracture area; modified from Zangerl et al. (2008). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 19

 19 Initial aperture maps for all 6 faults F1, F2, F3, F4, F5 and F6. The transparent area represents initial contact points without any normal loading. . . . . . . . . . . . . 3.20 The stiffness as a function of closure to obtain the stiffness characteristic for all 6 faults. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XXIII 3.21 Comparison of different approaches to obtain the normal stress against fracture closure for all 6 faults. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.22 Linear behavior of the normal stiffness against the normal stress. (a) linear regression of stiffness-stress data of the stress-linear stiffness approach (χ = 23.2 m -1 ) for the rigid-plastic (RP) model; (b) linear regression of stiffness-stress data of the asperities yield approach (χ = 23.1 m -1 ) for the rigid-plastic (RP) model; (c) linear regression of stiffness-stress data of the elasto-plastic (EP) model. . . . . . . . . . 3.23 The scaling relationship between fault volume and fault closure showing that the fault volume V f can be approximated by A(6h rms -u f ). . . . . . . . . . . . . . . 3.24 Height steps of the progressive closure of the fracture opening for fault F1 due to 1) increase of applied normal stress related to (a) -(d): case c0 to case c4; 2) precipitation of fracture filling material as an homogeneous coating of the surfaces corresponds to (e) -(l): case c4 to case c11. . . . . . . . . . . . . . . . . . . . . . 3.25 Evolution of the permeability in the x-and y-direction as a function of the fault volume for the studied 6 faults. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Fluid flow model of the fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Linear elastic traction-separation cohesive model . . . . . . . . . . . . . . . . . . 4.3 (a) Geometry of the conceptual model. The red dot represents the injection point. The blue dashed line and blue solid line indicate hydraulic fracturing paths and preexisting faults, respectively. Each segment is 10 meters long. The two parallel black dashed lines denote the permeable fracture zone. (b) Schematic diagram of induced fault slip during and after injection. The point Q s represents the stop injection position. d is the quiescence space between Q s and the fault F2. The red and blue color denotes the damage variable D = 1 and D = 0, which indicates fully damage and zero damage respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Pore pressure distribution during injection and after termination. (a) onset of the injection; (b) F1 activation; (c) the second hydraulic fracture initiation; (d) stop injection; (e) F2 activation; (f) pore pressure diffusion in the porous rock matrix . . 4.5 Pore pressure evolution at the injection point and the center of fault F1 and F2 . . . 4.6 Evolution of (a) pressure and (b) stress of F2 during shut-in stage . . . . . . . . . . 4.7 Pore pressure evolution of Node 2 in the compressional zone . . . . . . . . . . . . 4.8 Evolution of (a) total slip and (b) slip rate of F1 and F2 over time . . . . . . . . . . 4.9 Hydraulic fracturing front and pressure front propagation. (a) case A2; (b) case A8; (c) case A12; (d) case with no leak-off . . . . . . . . . . . . . . . . . . . . . . . . Deep geothermal energy has been extensively exploited and utilized worldwide, with a typical depth of over 500 m. As a clean and renewable energy, it can be used as direct heat at lower temperatures (60 -150 °C, e.g., Lopez et al. (2010); Agemar et al. (2014); Limberger et al. (2018); Romanov and Leiss (

Figure 1 . 1 :

 11 Figure 1.1: Enhanced Geothermal System (EGS); source: https://www.egs-energy.com/projects/

Figure 1 . 4 :

 14 Figure 1.4: 3D view of the induced seismicity in Soultz-Sous-Forêts, modified from (Sausse et al., 2010); The circles indicate seismic events

Figure 1

 1 Figure 1.8: Coupled deformation/pressure cohesive element with 6 nodes. (a) undeformed configuration; (b) deformed configuration. Red dots at corners denote displacement/pressure degrees of freedom, while blue dots at the middle-edge indicate only pressure degrees of freedom.

Figure 2 . 1 :

 21 Figure 2.1: Sketch of the whole simulation scheme: fracture surface/aperture generation; fracture closure; use of the Navier-Stokes equations for steady state flow through a volume-represented fracture and the pressure diffusion equation for transient flow along a surface-represented fracture; and determination of the effective hydraulic diffusivity

Figure 2 . 2 :

 22 Figure 2.2: Statistical invariance illustrated by magnifying portions of the profiles from the Corona Heights fault surface in (a) parallel (H ∥ = 0.65 ± 0.04) and (b) perpendicular (H ⊥ = 0.83 ± 0.03) to the slip direction; derived from[START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF].

  2.4). The largest height fluctuations are on the order of 3h rms ≈ ±0.3 m for a lateral extension of 512 m.

Figure 2 . 4 :

 24 Figure 2.4: Histogram of the height distribution of our generated isotropic self-affine surface at the 512 m scale, with Hurst exponent H = 0.8 and roughness amplitude h rms = 0.1 m, consistent with the field observations of Candela et al. (2012) (see blue circles in Fig. 1). The largest height fluctuations are on the order of 3h rms ≈ ±0.3 m for a lateral extension of 512 m

  Figure 2.5: 2D cross-section along the (x, z)-plane and (y, z)-plane (with a vertical magnification 500 times on the z-axis scale with respect to the x, y-axis) of the self-affine fracture aperture h(x, y) at two steps of the fracture closing. (a) Fully open fracture without contact. z = z 1 and z = z 2 represent the mean plane of the top surface s 1 (x, y) and the bottom surface s 2 (x, y), respectively. d m = z 1 -z 2 denotes the mean aperture. (b) Partial open fracture after the contact of some asperities. The shaded area indicates contact regions where the apertures are set to zero; d m here denotes the mean separation of the open area (light blue zones)

Figure 2

 2 Figure 2.6: Schematic diagram of the two representations of a rough fracture with L = 64 m. (a) Volume-represented fracture with 52,128 prism elements (8 layers along z); (b) surface-represented fracture with 4,096 quadrilateral elements. Both representations share the same mesh size of 1 m × 1 m along x and y

  Figure 2.7: Schematic of the initial conditions and boundary conditions applied in the x-and ydirections for each rough fracture. (a) Flow along the x-axis; (b) flow along the y-axis. The inlet and outlet boundaries are marked as red and blue, respectively

Fig. 2

 2 Fig. 2.8 shows the comparison of the fluid flow between volume-represented fractures and surface-represented fractures for steady-state flow. The results are given in terms of the ratio of the hydraulic aperture d h to the mean aperture d m as a function of the h rms -normalized mean aperture.As the roughness increases, a more significant deviation from the parallel plate flow is observed. The same simulations performed in the y-direction indicate anisotropy of the fracture flow. The good consistency between the two kinds of flow illustrates the effectiveness of using 2D surface elements to approximate 3D volumetric elements in fluid flow simulations at low Reynolds numbers (Re ≪ 1).

Figure 2

 2 Figure 2.8: Comparison of the fluid flow between volume-and surface-represented rough fractures. The aperture evolution is given by the ratio of the hydraulic aperture d h to the mean aperture d m as a function of the h rms -normalized mean aperture d m /h rms . The horizontal dashed line denotes the parallel plate model in which d h /d m = 1, while the solid vertical line indicates the divide where the fracture fully opens (right part) and comes into contact (left part)

Figure 2 .

 2 Figure 2.10: Map of the aperture distribution with H = 0.8, h rms = 0.1 m and d m = 0.3 m at the 512 m scale. The fluid flows from left to right in the x-direction and from bottom to top in the y-direction

Figure 2 .

 2 Figure 2.11: Schematic diagram of the pressure field evolution in a rough fracture and the parallel plate model with the same size (512 m) and the same mean aperture (0.3 m) at t 1 = 5 × 10 -7 s, t 2 = 2 × 10 -6 s and t 3 = 5 × 10 -6 s. (a)-(c) Fluid flow along the x-axis for the parallel plate model; (d)-(f) fluid flow along the x-axis for the rough fracture; and (g)-(i) fluid flow along the y-axis for the rough fracture

  Figure 2.13: Pressure diffusion along the x-axis (averaged along y) with d m /h rms = 3 at t 1 = 5 × 10 -7 s, t 2 = 2 × 10 -6 s and t 3 = 5 × 10 -6 s. (a) Comparison between the rough fracture and the parallel plate model with the same d m = 0.3 m; (b) comparison between the rough fracture and the parallel plate model with the best fitting hydraulic diffusivity D e = 0.423D m

  Figure 2.15: (left) Evolution of the contact area as a function of the normal closure (defined as the fracture's vertical displacement, which is positive along negative z direction) and mean aperture d m . (right) Evolution of the hydraulic aperture d h as a function of the mean aperture d m

Figure 2 .

 2 Figure 2.16: Log-log graph of the effective hydraulic diffusivity D e as a function of the hydraulic aperture d h . The dotted line corresponds to the parallel plate model (Eq. 2.19)

  Figure 2.17: Schematic diagram of the pressure diffusion evolution along the x-axis of the rough fracture with different contact areas. (a)-(c) t 1 = 2 × 10 -6 s, t 2 = 4 × 10 -6 s and t 3 = 1 × 10 -5 s with 5.43% contact; (d)-(f) t 1 = 2 × 10 -4 s, t 2 = 4 × 10 -4 s and t 3 = 0.001 s with 33.31% contact; and (g)-(i) t 1 = 0.002 s, t 2 = 0.004 s and t 3 = 0.01 s with 50.51% contact

Figure 2 .

 2 Figure 2.19: Ratio of the effective hydraulic diffusivity and the parallel plate fracture diffusivity D e /D m as a function of the h rms -normalized mean aperture d m /h rms . The vertical black line corresponds to the first asperity contacts when closing the fracture at d m /h rms ≈ 3

  Figure 2.20: Two examples of the evolution of the relative hydraulic diffusivity D e /D m for two other aperture fields with the same Hurst exponent H = 0.8 and the same RMS h rms = 0.1 m when changing the seed of the self-affine surface generator

Figure 2 .

 2 Figure 2.21: Pore pressure averaged along the y-axis as a function of t/t * at x = L/2 in the case of the parallel plate model with d m = 0.3 m

Fig. ( 2

 2 Fig. (2.22) shows an example of the least square error when searching for the effective hydraulic diffusivity D e = 0.423D m for d m /h rms = 3 and an imposed pressure drop along the x-axis. By stepwise selecting the resolution of the test diffusivity as ∆D = 0.1D m , ∆D = 0.01D m and ∆D = 0.001D m , D e1 = 0.4D m , D e2 = 0.42D m and D e3 = 0.423D m are obtained, respectively.

  Figure 2.22: Least square error as a function of D e /D m . (a) D e1 = 0.4D m with the step of the test diffusivity ∆D = 0.1D m ; (b) D e2 = 0.42D m with ∆D = 0.01D m ; and (c) D e3 = 0.423D m with ∆D = 0.001D m

Figure 3 .

 3 Figure 3.1: Sketch of the model: a partially open fault (colored horizontal cross-section) is embedded in a 3D block (in dark gray) of size 256 × 256 × 256 m 3 loaded by a normal stress σ n . (a) The contact areas (in light gray) along the fracture support the applied normal stress with variable local amplification of the stress (reddish zones). The open part of the fracture where fluid circulates is shown in blue. (b) 3D view of the fluid circulation in the fracture (colored horizontal cross-section) and in the pore space of the matrix (white arrows) when a pressure drop is applied in the x-direction (the vertical thick, blue layers indicate the pressure field in the matrix). The contact areas are transparent here.

Figure 3 . 2 :

 32 Figure 3.2: Generation of the fault topology and aperture map: (a) aperture map with only a few contact points; (b) power density function (pdf) of the aperture distribution with a Gaussian fit (average aperture is h = 0.27 m; RMS is h rms = 0.09 m).

Figure 3 . 3 :

 33 Figure 3.3: Cross-section through the block along the x-direction with a large vertical magnification to show the fracture closure process due to two different mechanisms. (a) Initial unloaded fracture aperture: fracture surfaces are just touching one another at a few contact points. (b) Phase 1 -mechanical closure (contact area < 20%): the contact area increases because of a normal loading (imposed normal displacement) applied on the top of the host rock; (c) asperities that overlap following a 'rigid-plastic' behavior. (d) Phase 2 -sealing closure due to mineral deposits (contact area > 20%). Successive layers of sealing material are assumed to grow homogeneously with a constant thickness from the border to the center of the open space (syntaxial mode). The open space for fluid flow is shown in light blue, whereas the host rock is shown in light gray.

Figure 3 . 4 :

 34 Figure 3.4: Definition of the fracture displacement u f . u f = u b -u m , where u b is the displacement of the fracture-matrix system, while u m is the displacement of the intact rock matrix under the same applied normal loading σ n .

Figure 3 . 5 :

 35 Figure 3.5: Evolution of the normal stress over the fracture closure (imposed displacement). The slopes of the red lines at each step of the mechanical/sealing closure correspond to its normal stiffness.

Figure 3 .

 3 Figure 3.6: Evolution of the relative contact area A c /A as a function of the normalized normal displacement imposed on the fractured block (u f -h)/h rms , where the mean aperture is h and the RMS is h rms . The pink domain corresponds to the mechanical closure of the fracture, assuming an overlapping model of contact (i.e., purely plastic rheology of the asperities). The light blue domain corresponds to the sealing closure due to homogeneous mineral deposits in the open fracture space. The evolution of the contact area can be fitted by the cumulative distribution function of a normal distribution y = (1 + erf (x/ √ 2))/2. The transition from a mechanical to a sealing closure is defined at 20% of contact, which corresponds to (u f -h)/h rms ∼ -1.

Figure 3 .

 3 Figure 3.7: Cross-section view of the initial conditions and boundary conditions for (a) the rock matrix and (b) the fracture. Transfer functions at the fracture-rock interface are also illustrated.

Figure 3 .

 3 Figure 3.8: (a) Finite elements of the open fracture space (in green) and of the rock mass (in brown). Coarse-graining of the elements develops when moving away from the fracture. (b) Details of the meshing of the open space of the fracture by prism elements.

Figure 3

 3 Figure 3.9: a) Fitting of the fracture normal stiffness as a function of fracture closure during the mechanical stage for F1 using the stress-linear stiffness approach. b) Linear regression of the stiffnessstress data of the stress-linear stiffness approach.

  m a l s t r e s s ( M P a ) C l o s u r e ( m ) R P m o d e l ( s t r e s s -l i n e a r a p p r o a c h ) R P m o d e l ( y i e l d -s t r e s s s p = 4 5 . 5 M P a ) E P m o d e l ( s p = 4 6 . 7 M P a )

  .16) with κ 0 = e 6 • σ ref n /h rms , the stiffness at full closure of the fault (V f = 0), and V c f = A • h rms a characteristic volume of the fault. The reference normal stress σ ref n obtained here is approximately 0.2 MPa.

Figure 3 .

 3 Figure 3.11: Log-linear plot of the fracture normal stiffness against the fault volume for all studied faults. The pink curve denotes the prediction from Eq. (3.16): κ ≈ κ 0 exp (-V f /V c f ).

Figure 3 .

 3 Figure 3.12: Map view of the flow inside fault F1. (a) Pressure drop is applied along the x-direction. (b) Pressure drop is applied along the y-direction. The white zones indicate the contact regions.

R

  h a n i c a l c l o s u r e -X S e a l i n g c l o s u r e -X M e c h a n i c a l c l o s u r e -Y S e a l i n g c l o s u r e -Y k = ( 1 . 2 4 3 e -4 ) / A 2 * ( V f -V f p e r c ) 2

Figure 3 .

 3 Figure 3.13: Evolution of the permeability as a function of the percolation fault volume (the fault volume minus the fault volume at the percolation threshold (V f -V perc f

Figure 3 .

 3 Figure 3.14: Fracture stiffness as a function of the square root of permeability. The pink line denotes the prediction from Eq. (3.19) with a representative fault volume of 3000 m 3 at the percolation threshold.

Figure 3 .

 3 Figure 3.15: Evolution of the permeability ratio (k y /k x ) as a function of the fault volume during the mechanical & sealing closure process for the three studied faults.

Figure 3 .

 3 Figure 3.18: (a) Power spectral density C(q) of the fault topography F1 along both directions (X and Y). The dashed lines indicate the domain of measured surfaces from the field and laboratories[START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF] (top line is C(q) = 2 × 10 -2 q -3 , medium line is C(q) = 2 × 10 -4 q -3 and bottom line is C(q) = 2 × 10 -6 q -3 ); (b) RMS h rms evolution with scale L along X and Y for the same synthetic fault.

Figure 3 . 3 RFigure 3 .Figure 3 .

 3333 Figure 3.19: Initial aperture maps for all 6 faults F1, F2, F3, F4, F5 and F6. The transparent area represents initial contact points without any normal loading.

Figure 3 .Figure 3 .

 33 Figure 3.22: Linear behavior of the normal stiffness against the normal stress. (a) linear regression of stiffness-stress data of the stress-linear stiffness approach (χ = 23.2 m -1 ) for the rigid-plastic (RP) model; (b) linear regression of stiffness-stress data of the asperities yield approach (χ = 23.1 m -1 ) for the rigid-plastic (RP) model; (c) linear regression of stiffness-stress data of the elastoplastic (EP) model.

Figure 3 .

 3 Figure 3.24: Height steps of the progressive closure of the fracture opening for fault F1 due to 1) increase of applied normal stress related to (a) -(d): case c0 to case c4; 2) precipitation of fracture filling material as an homogeneous coating of the surfaces corresponds to (e) -(l): case c4 to case c11.

Figure 3 .

 3 Figure 3.25: Evolution of the permeability in the x-and y-direction as a function of the fault volume for the studied 6 faults.

Figure 4 . 1 :

 41 Figure 4.1: Fluid flow model of the fracture

Figure 4 .

 4 Figure 4.2: Linear elastic traction-separation cohesive model

Figure 4 . 3 :

 43 Figure 4.3: (a) Geometry of the conceptual model. The red dot represents the injection point. The blue dashed line and blue solid line indicate hydraulic fracturing paths and pre-existing faults, respectively. Each segment is 10 meters long. The two parallel black dashed lines denote the permeable fracture zone. (b) Schematic diagram of induced fault slip during and after injection. The point Q s represents the stop injection position. d is the quiescence space between Q s and the fault F2. The red and blue color denotes the damage variable D = 1 and D = 0, which indicates fully damage and zero damage respectively.

Figure 4 . 4 :Figure 4

 444 Figure 4.4: Pore pressure distribution during injection and after termination. (a) onset of the injection; (b) F1 activation; (c) the second hydraulic fracture initiation; (d) stop injection; (e) F2 activation; (f) pore pressure diffusion in the porous rock matrix

Figure 4 .

 4 Figure 4.6: Evolution of (a) pressure and (b) stress of F2 during shut-in stage

Figure 4

 4 Figure 4.7: Pore pressure evolution of Node 2 in the compressional zone

Figure 4 . 8 :

 48 Figure 4.8: Evolution of (a) total slip and (b) slip rate of F1 and F2 over time

Figure 4 . 9 :Figure 4 .

 494 Figure 4.9: Hydraulic fracturing front and pressure front propagation. (a) case A2; (b) case A8; (c) case A12; (d) case with no leak-off

Figure 4 .

 4 Figure 4.11: Seismic slip and slip rate of (a) F1 and (b) F2 with different friction coefficient

Figure 4 .

 4 Figure 4.12: Seismic slip and slip rate of (a) F1 and (b) F2 with different fault angle

Figure 4 .Figure 4 .

 44 Figure 4.13: Injection scenarios with a constant injection volume: (a) constant injection rate and (b) cyclic injection rate

Figure 4 .

 4 Figure 4.15: Earthquake parameters scaling for fault slip, fault size, stress drop and magnitude, modified from[START_REF] Zoback | Earthquake triggering and large-scale geologic storage of carbon dioxide[END_REF]

Figure 4 .

 4 Figure 4.16: Moment magnitude as a function of rupture area, modified from[START_REF] Yoon | Discrete bonded particle modelling of fault activation near a nuclear waste repository site and comparison to static rupture earthquake scaling laws[END_REF]

Figure 4 .

 4 Figure 4.18: Comparison between analytical solutions and numerical simulations and the relative error

  1.1 Enhanced Geothermal System (EGS); source: https://www.egs-energy.com/projects/ 1.2 Simplified flow chart on harnessing geothermal energy by enhancing geothermal reservoirs properties and the associated techniques/issues, modified from (Kluge, 2021) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3 3D view of the main fault zone in Soultz-Sous-Forêts, modified from

Table 1 .

 1 1: Comparison of common contact models in rock engineering.

	Contact models	Descriptions	Asperity interactions
	HGW model	Hertzian/Greenwood-Williams contact theories	No
	Hopkins model	Elastic deformation; Cylindrically shaped asperities	Yes
	Interpenetration model Unconditionally deformable asperities that are merged	No
	FFT-based model	Elastic or elastic-plastic deformation	Yes

Table 3 .

 3 1: Material properties used in the numerical simulations.

	Property	Symbol	Value	Unit
	Young's modulus	E	25.0× 10 9	Pa
	Poisson ratio	ν	0.25	-
	Matrix permeability	k m	2.0 × 10 -16	m 2
	Porosity	ϕ	0.085	-
	Solid bulk modulus	K s	4.15 × 10 10	Pa
	Fluid modulus	K f	2.2 × 10 9	Pa
	Fluid viscosity	η	0.001	Pa•s
	Solid density	ρ s	2500	kg/m 3
	Fluid density	ρ f	1000	kg/m 3

Table 4 .

 4 1: Initial input material properties

	Property	Symbol Value	Unit	Notes
	Young's modulus	E	37.5	GPa	(Meyer et al., 2017)
	Poisson ratio	v	0.25	-	(Meyer et al., 2017)
	Porosity	ϕ	0.01	-	-
	Rock permeability	k	1.0e-16	m 2	-
	Permeability near fault	k f	1.0e-15	m 2	-
	Fracture Energy	G c	80	Pa•m	-
	Tensile strength	T c	2	MPa	-
	Penalty stiffness	K 0	800	GPa/m	(Turon et al., 2007)
	Biot's coefficient	b	1.0	-	-
	Friction coefficient	µ 1 , µ 2	0.2	-	-
	Fault orientation	α 1 , α 2	30	°-
	Initial aperture	w	0.4	mm	(Meyer et al., 2017)
	Fluid viscosity	η	0.001	Pa•s	water
	Leak-off coefficient	c	1.0e-11 m 3 /kPa•s	(Yao, 2012)
	Saturated degree	s	100	%	-
	Initial maximum stress	σ 1	36	MPa	(Meyer et al., 2017)
	Initial minimum stress	σ 3	29	MPa	(Meyer et al., 2017)
	Initial pore pressure Injection rate	p 0 Q	23.7 0.001	MPa m 3 •s -1	(Meyer et al., 2017) -
	Injection time	t i	160	s	-
	Shut-in time	t s	600	s	-

Table 4 .

 4 2: Summary of key parameters and results with different matrix permeability

	Case * Injection	Permeability	Slip onset of	Slip onset of
		time (s)	(m 2 )	F1 (s)	F2 (s) **
	A1	150	1.0e-15	38.12	-
	A2	160	1.0e-15	38.12	360.1
	A3	165	1.0e-15	38.12	188.1
	A4	170	1.0e-15	38.12	74.5
	A5	120	5.0e-15	42.01	-
	A6	130	5.0e-15	42.01	199.4
	A7	140	5.0e-15	42.01	86.0
	A8	160	5.0e-15	42.01	34.7
	A9	100	1.0e-14	35.83/47.81	-
	A10 110	1.0e-14	35.83/47.81	141.4
	A11 130	1.0e-14	35.83/47.81	25.2
	A12 160	1.0e-14	35.83/47.81	150.8 (in)
	* all cases have fault orientation α = 30°, friction coefficient µ 1 = µ 2 = 0.2, injection rate Q = 0.001 m 3 /s;
	** the slip time of F2 specifies the time after termination; F2 slip during injection
	is marked as '(in)'. All the same below.

Table 4 .

 4 3: Summary of key parameters and results with different friction coefficients

	Case * Injection	k f (m 2 )	µ 1 µ 2 Slip onset of	Slip onset of
		time (s)		F1 (s)	F2 (s)
	A2	160	1.0e-15	0.2 0.2 38.12	360.1
	B1	160	1.0e-15	0.4 0.4 38.12	-
	B2	160	1.0e-15	0.6 0.6 38.12	-
	A12 160	1.0e-14	0.2 0.2 35.83	150.8 (in)
	B3	170	1.0e-14	0.2 0.4 35.83	-
	B4	180	1.0e-14	0.2 0.4 35.83	75.7
	B5	160	1.0e-13	0.2 0.4 25.62	145.9 (in)
	B6	175	1.0e-13	0.2 0.6 25.62	-
	B7	180	1.0e-13	0.2 0.6 25.62	10.2

* all cases have fault orientation α = 30°, µ 1 = 0.2, injection rate Q = 0.001 m 3 /s;

Table 4 .

 4 5: Summary of key parameters and results with different injection schemes all cases have k f =1e-15 m 2 , µ 1 = µ 2 = 0.2, fault orientation α = 30°.

	Case * Injection	Injection rate	Slip onset of	Slip onset of
		time (s)	(m 3 /s)	F1 (s)	F2 (s)
	A2	160	0.001	38.12	360.1
	E1	100	0.0016	21.85	97.7
	E2	200	0.0008	47.24	410.0
	E3	320	0.0005	73.17	-
	E4	40×4	0.001	38.12	-
	E5	30+50+40×2 0.001	47.83/80.27	-

*

Table 4 .

 4 6: Input material properties for the KGD model

	Property	Symbol	Value	Unit
	Young's modulus	E	17	GPa
	Poisson ratio	v	0.2	-
	Fracture Energy	G c	120	Pa•m
	Tensile strength	τ c	1.25	MPa
	Biot's coefficient	b	0.75	-
	Biot's Modulus	M b	68.7	MPa
	Permeability	k	1.0×10 -16	m 2
	Porosity	ϕ	0.2	-
	Fluid viscosity	η	0.0001	Pa•s
	Injection rate	Q	0.001	m 3 •s -1

Acknowledgements

First of all, I would like to thank the China Scholarship Council (CSC) for supporting my PhD thesis, which has completely changed the trajectory of my life. I still vividly remember the various efforts I made several years ago to prepare my application for the PhD project. ACKNOWLEDGEMENTS I really miss the time spent with the following guys: Minghe, Xianyang, Peng, Guodong, Pan and Shangyao. We had a lot of fun and memorable experience together, especially considering that most of them happened during the period of lockdown and curfew. Thanks for our mutual support.

I would like to thank my colleagues at EOST/ED413 who are listed below but not limited to:

Acknowledgements

QD is funded by the China Scholarship Council (Grant No. 201808510128). This work was performed under the framework of the Laboratory of Excellence LABEX ANR-11-LABX-0050-G-EAU-THERMIE-PROFONDE and the Interdisciplinary Thematic Institute GeoT, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm, supported by IdEx Unistra (ANR 10 IDEX 0002) and by SFRI STRAT'US project (ANR 20 SFRI 0012) under the framework of the French Investments for the Future Program. Acknowledgements Funding: This work was supported by the European Commission Horizon 2020 research and innovation program (DESTRESS) [grant number 691728] and the China Scholarship Council [grant number 201808510128]. Special thanks to the free platform 'contact.engineering' for comput-

Abstract

We develop a fully coupled hydro-mechanical model to simulate fault slip due to fluid injection. We consider the interaction between a hydraulic fracture and pre-existing faults as well as the fluid exchange between the fracture/fault and the porous matrix. In order to consider a pressure diffusion mechanism, we set a relatively high permeability around the stimulated path. Our parametric study shows that a couple of factors affect the fault activation and its slip behavior such as fault properties, friction properties and injection scenario. We observe that pore pressure diffusion induces poroelastic stress change, which are able to produce shut-in events with a time and space lag. This mechanism also affects the slip behavior during injection in particular when the surrounding permeability is high (e.g., up to 1e-13 m 2 ), and provides a new insight into understanding the occurrence of stronger seismic events after shut-in compared to the injection phase. In addition, we show that small perturbations may trigger large seismic fault slip which highlights the key role of the initial fault stress state. The results have profound implications for deep fluid injection related engineering as well as for soft cyclic injection strategies aiming to mitigate the risk of large earthquakes. that this is because the pre-existing fault modifies local stress distributions and the latter depends primarily on the fault orientation. In addition, as the fault dip decreases, the triggered slip and slip rate decrease (Fig. 4.12(a)), in particular, when the fault angle is 10°(case C4), the slip exhibits mainly aseismic behavior and the averaged slip rate is only about 2.3 mm/s.

During the shut-in stage, except case A2, only case C1 with α = 45°induces F2 activation at 85.7 s after termination in comparison with 360.1 s in case A2. For the 20°model, the shut-in effect of F2 appears when the surrounding permeability is enhanced from k f 1 = 1e-15 m 2 to k f 3 = 1e-14 m 2 (case C3). For the 10°model, F2 can be activated not only with even a higher surrounding permeability (1e-13 m 2 ), but a higher injection rate 1.6 m 3 /s in case C3 (we will show that a higher injection rate is more likely to trigger shut-in slip in the next section). The results demonstrate that a higher fault orientation is more sensitive to F2 reactivation after shut-in. Furthermore, it tends to induce a larger shear slippage with a higher slip rate (Fig. 4.12(b)).

Injection scheme

In this section, we investigate the effect of injection scheme, by considering case A2 as the reference case (Tab. 4.5).

In a first step, we change the injection rate to 0.0016 m 3 /s, 0.0008 m 3 /s and 0.0005 m 3 /s , corresponding cases E1, E2 and E3 (Fig. 4.13a), respectively. Worth noting is that we adjust the duration of the injection so as to maintain the same injection volume for all experiments. As expected, lowering the injection rate extends the time to reactivate F1. When Q = 0.0016, 0.001 and 0.0008 m 3 /s, F2 can be reactivated but with a longer delay at 97.7 s, 360.1 s and 410 s after shut-in, respectively. No slip of F2 is observed when the injection rate further reduces to 0.0005 m 3 /s. Either during in-