Use of Mohr Diagrams to Predict Fracturing in a Potential Geothermal Reservoir

: Inferences have to be made about likely structures and their effects on ﬂuid ﬂow in a geothermal reservoir at the pre-drilling stage. Simple mechanical modelling, using reasonable ranges of values for rock properties, stresses and ﬂuid pressures, is used here to predict the range of possible structures that are likely to exist in the sub-surface and that may be generated during stimulation of a potential geothermal reservoir. In particular, Mohr diagrams are used to show under what ﬂuid pressures and stresses different types and orientations of fractures are likely to be reactivated or generated. The approach enables the effects of parameters to be modelled individually, and for the types and orientations of fractures to be considered. This modelling is useful for helping geoscientists consider, model, and predict the ranges of mechanical properties of rock, stresses, ﬂuid pressures, and the resultant fractures that are likely to occur in the sub-surface. Here, the modelling is applied to folded and thrusted greywackes and slates, which are planned to be developed as an Enhanced Geothermal System beneath Göttingen. are drawn approximately perpendicular to bedding. We predict the following sequence as is increased.


Introduction
A study is being undertaken to predict the structures and fluid flow behaviour in the sub-surface as part of a proposed geothermal project at Göttingen in Germany [1][2][3]. The reservoir rocks are the Devonian and Carboniferous metasedimentary sequence thought to occur at a depth of~1.5 km, beneath a cover of Permian and Mesozoic sedimentary rocks [4]. Göttingen is one of the sites within the European funded project Multidisciplinary and multi-context demonstration of EGS exploration and Exploitation Techniques and potentials (MEET) to develop enhanced geothermal systems (EGS) across Europe (e.g., [5][6][7]). The Triassic Bunter Sandstone is also being considered as a potential medium-deep geothermal reservoir (i.e., at depths between about 200 m and 1000 m; e.g., [8]), including heat storage options [2], and this is also discussed here.
Only two wells have yet penetrated the metasedimentary rocks in the Göttingen area, these not extending far below the base Zechstein, and so limited well data are available. Additionally, the limited seismic data available show poor resolution beneath the Zechstein salt (e.g., [9,10]). Two seismic lines were shot in 2015, these being 10 km and 11 km long and to about 5 km deep. These seismic lines allow interpretation of the post-Carboniferous sedimentary rocks in the north-south striking Leinetal Graben [3], but do not enable reliable interpretation of the sub-Zechstein rocks. The metasedimentary rocks exposed in the Harz Mountains,~40 km NE of Göttingen, are being used to predict the rocks and structures that may occur in the sub-surface at Göttingen, and to gain data input for discrete fracture network (DFN) modelling. Exposures of the Bunter Sandstone occur in the area around Göttingen (e.g., [11][12][13]).

1.
Which lithologies are most likely to fracture? 2.
What stresses and fluid pressures are needed for the reactivation of pre-existing fractures or the development of new fractures? 3.
Which orientations and types of fractures are most likely to be reactivated? 4.
Will reactivated or new fractures show shear or extension? 5.
What effects do heterogeneities (veins, joints, cleavage, bedding planes) have, and what are the different mechanical significances of veins vs. joints? 6.
What are the effects of Late Cretaceous and Tertiary exhumation and what amount of exhumation is needed to create joints?

Predictions about the Pre-Permian Geology beneath Göttingen
In the absence of well data or high-quality seismic data for the sub-Zechstein rocks, the rocks beneath the Permian Zechstein evaporites at Göttingen are predicted to be Devonian quartzitic sandstones and slates, and Carboniferous greywackes and slates (mainly Culm flysch deposits). Such rocks are exposed in the western Harz Mountains~40 km to the NE, in the Oberharz Anticline and the Culm Fold Zone, which belong to the par-autochthonous domain of the Harz Mountains (e.g., [3,20]) and in the Rhenish Massif~70 km to the SW (e.g., [21][22][23]). The fractures visible in these exposures include veins and joints [24], with cleavage being well-developed in the slates [25]. Small thrusts that appear to have displacements of up to a few metres are exposed [26], but larger thrusts are not wellexposed. The north-eastern boundary the Harz Mountains is marked by the Harznordrand  Table 1 shows features of these "basement" rocks that are likely to be important for geothermal energy. These predictions are expressed in terms of seven parameters that we consider important in describing a geothermal reservoir: lithology and inferred rheological behaviour, fluid type, stresses, fluid pressure, temperature, strains and existing structures, and geological history (cf. [30]).
More information is available about the Bunter Sandstone, which is exposed in the Göttingen area (e.g., [11,12]). Fractures exposed include normal faults [31], strike-slip faults [32], veins [33], sedimentary dykes [34], and joints [35]. There is limited well data, which suggests that the Bunter Sandstone has been affected by halokinesis of the Zechstein evaporites [1,36,37].  [24], with cleavage being well-deve to have displacements of up to a few m well-exposed. The north-eastern boun Harznordrand Fault, which was active d SE striking Mesozoic normal and obliqu of which contained economic Pb-Zn and  Table 1 shows features of these "ba geothermal energy. These predictions ar consider important in describing a geoth cal behaviour, fluid type, stresses, fluid p tures, and geological history (cf. [30]).
More information is available abou Göttingen area (e.g., [11,12]). Fractures faults [32], veins [33], sedimentary dyke which suggests that the Bunter Sandsto stein evaporites [1,36,37]. joints [24], with cleavage being well-developed in the slates [25]. Small thrusts that appear to have displacements of up to a few metres are exposed [26], but larger thrusts are not well-exposed. The north-eastern boundary the Harz Mountains is marked by the Harznordrand Fault, which was active during the Late Cretaceous and Tertiary [27]. NW-SE striking Mesozoic normal and oblique-slip faults occur in the Harz Mountains, some of which contained economic Pb-Zn and Ba-F deposits [28].  Table 1 shows features of these "basement" rocks that are likely to be important for geothermal energy. These predictions are expressed in terms of seven parameters that we consider important in describing a geothermal reservoir: lithology and inferred rheological behaviour, fluid type, stresses, fluid pressure, temperature, strains and existing structures, and geological history (cf. [30]).
) and extension fracturing ( Geosciences 2021, 11, x FOR PEER REVIEW 3 of 2 joints [24], with cleavage being well-developed in the slates [25]. Small thrusts that appea to have displacements of up to a few metres are exposed [26], but larger thrusts are no well-exposed. The north-eastern boundary the Harz Mountains is marked by th Harznordrand Fault, which was active during the Late Cretaceous and Tertiary [27]. NW SE striking Mesozoic normal and oblique-slip faults occur in the Harz Mountains, som of which contained economic Pb-Zn and Ba-F deposits [28].  Table 1 shows features of these "basement" rocks that are likely to be important fo geothermal energy. These predictions are expressed in terms of seven parameters that w consider important in describing a geothermal reservoir: lithology and inferred rheolog cal behaviour, fluid type, stresses, fluid pressure, temperature, strains and existing struc tures, and geological history (cf. [30]).
More information is available about the Bunter Sandstone, which is exposed in th Göttingen area (e.g., [11,12]). Fractures exposed include normal faults [31], strike-sli faults [32], veins [33], sedimentary dykes [34], and joints [35]. There is limited well data which suggests that the Bunter Sandstone has been affected by halokinesis of the Zech stein evaporites [1,36,37]. joints [24], with cleavage being well-developed in the slates [25]. Small thrusts that appear to have displacements of up to a few metres are exposed [26], but larger thrusts are not well-exposed. The north-eastern boundary the Harz Mountains is marked by the Harznordrand Fault, which was active during the Late Cretaceous and Tertiary [27]. NW-SE striking Mesozoic normal and oblique-slip faults occur in the Harz Mountains, some of which contained economic Pb-Zn and Ba-F deposits [28]. ing failure envelopes. Continuous line ("Cohesional") = Mohr failure envelope for intact rock, with tensile strength (T = 2 MPa), cohesion (SO = 4 MPa) and coefficient of internal friction (μ = 0.75). Dotted (parabolic) line ("Griffith") = Griffith criterion. Dashed line ("Cohesionless") = cohesionless failure (μ = 0.75) typical of slip on faults. Three stress states are shown, representing stable (), shear failure on cohesionless fractures () and extension fracturing (). (b) Mohr diagram showing the fields in which extension (), hybrid () and shear () fractures occur [29]. Table 1 shows features of these "basement" rocks that are likely to be important for geothermal energy. These predictions are expressed in terms of seven parameters that we consider important in describing a geothermal reservoir: lithology and inferred rheologi-cal behaviour, fluid type, stresses, fluid pressure, temperature, strains and existing struc-tures, and geological history (cf. [30]).
More information is available about the Bunter Sandstone, which is exposed in the Göttingen area (e.g., [11,12]). Fractures exposed include normal faults [31], strike-slip faults [32], veins [33], sedimentary dykes [34], and joints [35]. There is limited well data, which suggests that the Bunter Sandstone has been affected by halokinesis of the Zech-stein evaporites [1,36,37].  Table 1 shows features of these "basement" rocks that are likely to be important for geothermal energy. These predictions are expressed in terms of seven parameters that we consider important in describing a geothermal reservoir: lithology and inferred rheological behaviour, fluid type, stresses, fluid pressure, temperature, strains and existing structures, and geological history (cf. [30]).
) and shear ( joints [24], with cleavage being well-developed in the slates [25]. Small thrusts that appear to have displacements of up to a few metres are exposed [26], but larger thrusts are not well-exposed. The north-eastern boundary the Harz Mountains is marked by the Harznordrand Fault, which was active during the Late Cretaceous and Tertiary [27]. NW-SE striking Mesozoic normal and oblique-slip faults occur in the Harz Mountains, some of which contained economic Pb-Zn and Ba-F deposits [28].  Table 1 shows features of these "basement" rocks that are likely to be important for geothermal energy. These predictions are expressed in terms of seven parameters that we consider important in describing a geothermal reservoir: lithology and inferred rheological behaviour, fluid type, stresses, fluid pressure, temperature, strains and existing structures, and geological history (cf. [30]).
) fractures occur [29]. Table 1 shows features of these "basement" rocks that are likely to be important for geothermal energy. These predictions are expressed in terms of seven parameters that we consider important in describing a geothermal reservoir: lithology and inferred rheological behaviour, fluid type, stresses, fluid pressure, temperature, strains and existing structures, and geological history (cf. [30]).
More information is available about the Bunter Sandstone, which is exposed in the Göttingen area (e.g., [11,12]). Fractures exposed include normal faults [31], strike-slip faults [32], veins [33], sedimentary dykes [34], and joints [35]. There is limited well data, which suggests that the Bunter Sandstone has been affected by halokinesis of the Zechstein evaporites [1,36,37]. Table 1. Predictions about the geology of the potential geothermal reservoirs at Göttingen. See, for example, [38] for full information about deformation of the Bunter Sandstone in the Leinetal Graben. See [20,39] for further information about the geology of the Devonian and Carboniferous rocks of the Harz Mountains. Triassic sandstone (see Table 2 for mechanical properties) Devonian and Carboniferous greywackes and slates (see Table 2 for mechanical properties)

Fluid type
The chemistry and phase (liquid or gas) of the palaeoand present-day fluid(s)  MPa

Cohesion
The shear strength of a material when the stress normal to a shear surface is zero (e.g., [48,49]

Poisson ratio
The relationship between the tendency to shorten in one direction and the tendency to expand in another direction (e.g., [50]) 0.16 to 0.35 [51] 0.11 to 0.29 [41] 0.22 to 0.29 [42] Geostatic stress ratio * The ratio of the horizontal effective stress (σ H ) to the vertical effective stress (σ V ) (e.g., [52]). It gives the effect the overburden has on horizontal stresses. Influences the diameter of the Mohr circle. Values calculated using Equation (2) 0.19 to 0.54 0.125 to 0.41 [41] 0.28 to 0.41

Model Set-Up
A Microsoft Excel spreadsheet has been created to make the necessary calculations and to plot Mohr diagrams (e.g., Figure 2). For simplicity, we keep the Mohr diagram analysis two-dimensional, considering just the vertical stress and the horizontal stress in one direction. This approach is taken both because it simplifies the analysis and because the orientations and magnitudes of the horizontal stresses are presently unconstrained. We start with a reference state where the horizontal strains are zero, which is the uniaxial strain condition (e.g., [58]) and the horizontal effective stress (σ H ) is given by: where k 0 is the geostatic stress ratio [59] or coefficient of lateral earth pressure [60]. The geostatic stress ratio is the ratio of the horizontal effective stress to the vertical effective stress (i.e., k 0 = σ H /σ V ). For an isotropic elastic material: where ν is Poisson's ratio, which is generally in the range 0 to 0.5 (where 0.5 is incompressible). For most water-saturated rocks, is in the range 0.15 < ν < 0.4 (e.g., [61]). We consider the vertical stress to result from the weight of the overlying material and the horizontal stresses to result from the combined effects of the geostatic stress ratio and the applied tectonic stresses. We start the analysis with a base-case, in which fluid pressure is hydrostatic and there are no applied tectonic stresses. We then consider how changes in fluid pressure and/or apply tensile or compressional horizontal ("tectonic") stresses might lead to fracturing.

Mohr Diagrams, Stresses and Failure Envelopes
Mohr diagrams are commonly used to show the relationships between stresses, fluid pressure, and fractures (Figure 1a; e.g., [62]). They provide a convenient graphic representation of the effective stress states and of failure (e.g., [17]). Normal stress (σ N ) is plotted on the x-axis and shear stress (τ) on the y-axis, with the principal axes of stress (σ 1 = maximum compressive stress, σ 2 = intermediate compressive stress, σ 3 = minimum compressive stress) or principal axes of effective stress (σ 1 = maximum effective compressive stress, etc., where σ = σ − B.P F , where P F = fluid pressure and B is Biot's constant, which we will assume to be approximately 1). Since the magnitudes and orientations of the horizontal stresses are currently unknown, we consider stress in 2D, with σ 2 being ignored (e.g., [63]). We assume an Andersonian stress regime with one principal stress vertical (the overburden stress, σ V ) and the other horizontal stress (σ h ) [64]. In the base case models used here, the starting point is that there are no applied tectonic stresses, so σ V = σ 1 , and σ h = σ 3 ( Figure 2).  Table 3 Table 3   Any 2D state of effective stress can be represented by a circle that intersects the x-axis at σ 1 and σ 3 , with the centre of the circle representing the mean stress (σ Mean = (σ 1 + σ 3 )/2). Mohr diagrams are useful for illustrating the conditions under which fracturing may occur in a particular rock under specific conditions (i.e., the failure envelope), but only if the effective stresses are plotted.
Fracturing is classically attributed to conditions where the effective stress components exceed some critical value, usually termed the strength, which is thought to be a property of the material (e.g., [65]). Tensile failure occurs if: where T is the tensile strength. Compressive stresses are positive. Shear failure occurs if: where µ is the coefficient of (internal) friction and S 0 is the cohesion. These conditions are usually linked to produce a failure envelope (e.g., [66,67]), as in Figure 1a.
The composite failure criteria used in this work (e.g., Figures 1 and 2) combine a linear, Mohr-Coulomb envelope for shear fracture under compressive effective stresses, with a parabolic, plane Griffith envelope for tensile/hybrid fractures under tensile stress conditions. The two envelopes join at the τ-axis (σ = 0), where there is a discontinuity in the slopes of the failure envelopes. Continuity in the slope could be achieved by using the methods outlined in [68]. This type of composite failure envelope is consistent with a modification of the Griffith theory to account for closure and frictional behaviour under compression [69,70].
Following [71], the plane Griffith failure envelope in Mohr space (τ, σ ) is given by: (5) and the Mohr-Coulomb envelope by: where T = Tensile strength, S 0 = cohesion, and µ = coefficient of internal friction. Putting σ = 0 in both Equations (5) and (6) gives: Hence, the cohesion would be simply twice the tensile strength. The unconfined compressive strength (C 0 ), is a widely measured parameter representing the stress required for failure when σ 3 = 0 which for Mohr-Coulomb failure is given by: Substituting Equation (7) in Equation (8) gives: Using a value of µ ≈ 1, gives the following approximations for rocks: [71]. Sedimentary rocks typically show an exponential increase in C 0 with decreasing porosity [72].
We assume a lower bound to likely failure envelopes is given by the failure envelope shown by the "Cohesionless" line in Figure 1a, which is typical of slip on pre-existing cohesionless fractures (e.g., joints) (e.g., [73]). Such failure envelopes typically have coefficients of friction of 0.6 to 0.9 (average 0.75; e.g., [74]).
In Figure 1a, the effective stress state shown by Mohr circle  Table 1 shows features of these "basement" rocks geothermal energy. These predictions are expressed in t consider important in describing a geothermal reservoir cal behaviour, fluid type, stresses, fluid pressure, tempe tures, and geological history (cf. [30]).
is typical of a rock at depth of~2 km, with the maximum compressive stress being vertical (σ V ) and caused by the weight of the overlying material. The rock is under a hydrostatic fluid pressure (where P F /σ V =~0.4, based on the assumption that P F~2 0 MPa and σ V~5 0 MPa at 2 km) and a tectonic stress (−2 MPa, as may occur during regional extension). Note that this is in the stable region, but within~5 MPa of the failure envelope for rocks with cohesionless fractures. Figure 1a shows stress states that are stable ( to have displacements of up to a few metres are exposed [26], well-exposed. The north-eastern boundary the Harz Moun Harznordrand Fault, which was active during the Late Cretaceo SE striking Mesozoic normal and oblique-slip faults occur in th of which contained economic Pb-Zn and Ba-F deposits [28].  Table 1 shows features of these "basement" rocks that are geothermal energy. These predictions are expressed in terms of consider important in describing a geothermal reservoir: litholo cal behaviour, fluid type, stresses, fluid pressure, temperature, s tures, and geological history (cf. [30]).
), that would cause slip on a suitably orientated pre-existing cohesionless fractures ( joints [24], with cleavage being well-developed in the slates [25]. Small th to have displacements of up to a few metres are exposed [26], but large well-exposed. The north-eastern boundary the Harz Mountains is Harznordrand Fault, which was active during the Late Cretaceous and T SE striking Mesozoic normal and oblique-slip faults occur in the Harz M of which contained economic Pb-Zn and Ba-F deposits [28].  Table 1 shows features of these "basement" rocks that are likely to geothermal energy. These predictions are expressed in terms of seven pa consider important in describing a geothermal reservoir: lithology and i cal behaviour, fluid type, stresses, fluid pressure, temperature, strains an tures, and geological history (cf. [30]).
), and that would create tensile failure of the intact rock ( joints [24], with cleavage being well-developed in the slates [25]. Small thrusts that appear to have displacements of up to a few metres are exposed [26], but larger thrusts are not well-exposed. The north-eastern boundary the Harz Mountains is marked by the Harznordrand Fault, which was active during the Late Cretaceous and Tertiary [27]. NW-SE striking Mesozoic normal and oblique-slip faults occur in the Harz Mountains, some of which contained economic Pb-Zn and Ba-F deposits [28].  Table 1 shows features of these "basement" rocks that are likely to be important for geothermal energy. These predictions are expressed in terms of seven parameters that we consider important in describing a geothermal reservoir: lithology and inferred rheological behaviour, fluid type, stresses, fluid pressure, temperature, strains and existing structures, and geological history (cf. [30]).
In this paper, we use Mohr diagrams to illustrate the conditions under which a rock can go from a stable stress system (i.e., no active fracturing occurs) to an unstable stress system (i.e., fracturing occurs). The change from stable to unstable condition depends upon:
The stress state (σ), which is defined in terms of principal stresses, mean stress and differential stress [17]. Stresses are in turn controlled by factors, such as depth of burial (overburden), tectonic (horizontal) stresses, and other changes in the physical state of the material, such as expansion or contraction caused by temperature and volume change (e.g., [76]). Changes in stresses that lead to fracturing can either be by increasing [77] or reducing the [78] the applied compressive stresses; 3.
In the upper crust, fluid pressure in the pores and cracks combines with the applied stresses to produce an effective stress, where σ = σ − P F (e.g., [79][80][81]). In the absence of specific information, we use a Biot coefficient (B) of 1, where σ = σ − B.P F [82]. Changes in fluid pressure that can lead to fracturing can either be an increase in fluid pressure (e.g., [83]) or a reduction in fluid pressure, which can cause pore collapse (e.g., [84,85]). Pore collapse is not considered further in this paper.
The type and orientation of fractures developed can be predicted from the relationship between the Mohr circle for effective stress and the failure envelope (Figure 1b). Extension fractures typically develop perpendicular to the direction of σ 3 if σ 3 touches the failure envelope on the τ = 0 axis (Figure 1b  Table 1 shows features of these "basement" rocks that are likely to be important geothermal energy. These predictions are expressed in terms of seven parameters that consider important in describing a geothermal reservoir: lithology and inferred rheolo cal behaviour, fluid type, stresses, fluid pressure, temperature, strains and existing str tures, and geological history (cf. [30]).
), which generally requires a relatively low differential effective stress (σ Diff = σ 1 − σ 3 ). Shear fractures are predicted to develop if the Mohr circle touches the failure envelope in the compressive field (Figure 1b  Table 1 shows features of these "basement" rocks that are likely to b geothermal energy. These predictions are expressed in terms of seven par consider important in describing a geothermal reservoir: lithology and in cal behaviour, fluid type, stresses, fluid pressure, temperature, strains an tures, and geological history (cf. [30]).
), this typically requiring a relatively high σ Diff . Hybrid fractures develop by synchronous extension and shear (e.g., [29]), and occur if the Mohr circle touches curved part of the failure envelope within the tensile field of the Mohr diagram (Figure 1b  Table 1 shows features of these "basement" rocks that are likely to b geothermal energy. These predictions are expressed in terms of seven para consider important in describing a geothermal reservoir: lithology and in cal behaviour, fluid type, stresses, fluid pressure, temperature, strains and tures, and geological history (cf. [30]).
). Table 2 shows values from the literature for various rock mechanical parameters for the Bunter Sandstone, greywackes, and slates derived from triaxial tests. Table 2 also gives definitions and the significance of each parameter. Triaxial tests typically use intact hand specimen sized rock samples with no pre-existing fractures visible. As such, these values almost certainly over-estimate the strengths of the larger rock masses in the sub-surface, which generally have pre-existing fractures (e.g., [86]), and should be considered as upperlimit values of rock strength. We, therefore, use lower values for rock strength as inputs into the models ( Table 3). The failure envelope (e.g., Figure 2) for each rock type is controlled by the tensile strength (T), cohesion (S 0 ) and coefficient of internal friction (µ = tan ϕ, where ϕ = angle of internal friction) (Equations (3) and (4)). The most widely available rock strength parameter is the uniaxial compressive strength (UCS or C 0 ), determined directly from triaxial tests or estimated from geophysical log data [72,73].

Input Data, Assumptions, and Uncertainties
Note that the values presented in Table 3 are generally within the ranges of values for Devonian and Carboniferous greywackes and slates in the Harz Mountains and Belgium presented by [87]. Exceptions are the porosity (greywackes, mean value = 2.42%, range 0.04 to 7.35%, n = 88; slates, mean value = 2.21%, range 0.17 to 5.51%, n = 20), and the cohesion of the greywackes (mean value = 26.57 MPa, range 10.5 to 45.9 MPa, n = 5).
The Mohr circle for effective stress is controlled by the applied stresses and the fluid pressure. The vertical stress can be calculated from the densities of the overlying rocks (e.g., [73,88]). The magnitudes and orientations of these horizontal stresses are much harder to calculate (e.g., from well data; [89]) or predict (e.g., from tectonic stress tensors; [90]). Similarly, pore fluid pressure in rocks in the sub-surface is difficult to predict without well data. For example, [91] use sediment consolidation experiments and numerical models to predict fluid pressures.
In the absence of appropriate sub-surface data, the magnitudes of the horizontal stresses are considered in terms of the geostatic stress ratio of the rock (e.g., [92]). The pore pressures are discussed in terms of the hydrostatic pressure and any assumed overpressure.
The modelling presented here makes several assumptions. These assumptions are made to act as a starting point, to simplify the analysis and because of the various uncertainties: 1. An Andersonian stress system is assumed, i.e., with one of the principal axes of stress being vertical and the other two being horizontal [64]; 2.
The analysis is carried out in two-dimensions, considering just vertical stress and horizontal stress. This simplifies the analysis and is, we argue, justified at the predrilling stage of analysis because of the magnitudes and orientations of the horizontal stresses are currently unknown. Hydrofracture data from three wells in the region suggest a thrust regime with a maximum horizontal stress orientated~WNW-ESE [93]; 3.
The vertical stress is produced by the weight of overburden; 4.
The failure parameters used in the modelling (Table 3) are assumed to be representative of the rock properties in the sub-surface.

Base Case Models
We start with base case models ( Figure 2) to derive threshold conditions for the reactivation of fractures or the generation of new fractures. These base cases use the uniaxial strain condition (e.g., [58]), with no applied tectonic stresses. The fluids in pore spaces and open fractures are hydrostatically pressured, with the water table being at the ground surface. In such circumstances, the vertical effective stress creates a horizontal effective stress that is given by the geostatic stress ratio (Figure 2), with the input parameters used shown in Table 3. The base case models are used as a starting point for experiments in which fluid pressure is increased (as would occur either by natural overpressure or by hydraulic stimulation), and the horizontal (tectonic) stress either decreased or increased, until the model predicts either the reactivation of existing cohesionless fractures (hydraulic stimulation) or the generation of new fractures.
For most rocks, T << 100 MPa, with poorly-consolidated sediments having T → 0, whereas for shear failure 0.5 < µ < 1.5 (generally) and 0 < S 0 < 50 MPa. In this paper, we use ranges of T, S 0 , and µ that are appropriate for the rocks being considered.
An important point illustrated by Figure 2 is that, at typical K 0 values (0.3-0.5), the stress state is stable and no fracturing would be expected. This means that low K 0 values are needed to allow failure in the base case model, which simulates simple burial, with hydrostatic fluid pressure and no applied tectonic stresses. Any failure would largely be by shearing of pre-existing fractures, which may occur in the Devonian and Carboniferous rocks that are considered to occur beneath Göttingen.

A Range of Stress States for Fracturing
We consider a range of states of effective stress, seven of which are shown in Figure 3, together with their effect on fracturing. Note that the discussion is of changes in effective stresses, which can be caused by changes in applied stresses and/or fluid pressure. These states of stress have been modelled for the Bunter Sandstone and the Devonian and Carboniferous greywackes and slates by increasing the fluid pressure and by applying "tectonic" stresses to the base case models.

Effects of Key Parameters
A definition and the significance of each of the input parameters, along with representative values, are given in Table 2. For a given overburden stress, higher geostatic stress ratios produce larger horizontal stress and, hence, lower differential stress ( Figure 2). The effects of other key parameters are illustrated in Figure 4. The position and shape of the failure envelope is controlled by the tensile strength, cohesion, and coefficient of internal friction of the rock (Figure 4a). Vertical, and therefore horizontal, stress increases with depth as the overburden increases (Figure 4b). High values of compressional tectonic stress can cause the horizontal stress to exceed the vertical stress, tending to lead to shear failure, while tensile tectonic stress will increase the differential stress and may lead to the development of extension fractures (Figure 4c). Increases in fluid pressure move the Mohr circle for effective stress to the left, towards more tensile parts of the Mohr diagram (Figure 4d). Note, however, that a change in fluid pressure will cause a change in differential stress that is proportional to the geostatic stress ratio, and it is a common mistake to show a constant differential stress with changing fluid pressure (e.g., [94,95]). The base cases are close to cohesionless failure if lower K 0 values are used, which implies that fracturing will generally require some tectonic stress and/or overpressure, especially to initiate the development of new fractures. development of extension fractures (Figure 4c). Increases in fluid pressure move the Mohr circle for effective stress to the left, towards more tensile parts of the Mohr diagram (Figure 4d). Note, however, that a change in fluid pressure will cause a change in differential stress that is proportional to the geostatic stress ratio, and it is a common mistake to show a constant differential stress with changing fluid pressure (e.g., [94,95]). The base cases are close to cohesionless failure if lower K0 values are used, which implies that fracturing will generally require some tectonic stress and/or overpressure, especially to initiate the development of new fractures.    Table 3. (a) The base case Mohr circle (depth of 2 km, hydrostatic fluid pressure and no applied tectonic stress), with a representative cohesional failure envelope shown (Table 3, "Greywacke" and "Slate" columns). "Friction angle" refers to the slope the failure envelope. The Mohr diagrams presented in Figures 2-4 show two types of failure envelope. The "Cohesionless" lines on the Mohr diagrams are typical of slip on pre-existing unfilled fractures (e.g., joints; Figure 5a), while the "Cohesional" lines represent the behaviour of intact rock or fractures that are filled by cohesive material (e.g., veins or faults with cohesive gouge; Figure 5b). These different failure envelopes illustrate the different mechanical behaviours of such cohesionless fractures as joints and such cohesional fractures as fullyfilled veins and faults with clay gouge, highlighting the importance of distinguishing between different fracture types when using a field analogue for a geothermal reservoir.

Potential for Reactivating and Generating Fractures
Results of the models can be expressed in terms of the critical values of fluid pressure or horizontal (tectonic) stresses needed to reactivate existing cohesionless fractures, or to generate new fractures ( Table 4). Note that these are simply the results of the model, which is based on poorly constrained parameters and should only be taken as suggestions of what may occur under the conditions stated.  Table 3. (a) The base case Mohr circle (depth of 2 km, hydrostatic fluid pressure and no applied tectonic stress), with a representative cohesional failure envelope shown (Table 3, "Greywacke" and "Slate" columns). "Friction angle" refers to the slope the failure envelope. The Mohr diagrams presented in Figures 2-4 show two types of failure envelope. The "Cohesionless" lines on the Mohr diagrams are typical of slip on pre-existing unfilled fractures (e.g., joints; Figure 5a), while the "Cohesional" lines represent the behaviour of intact rock or fractures that are filled by cohesive material (e.g., veins or faults with cohesive gouge; Figure 5b). These different failure envelopes illustrate the different mechanical behaviours of such cohesionless fractures as joints and such cohesional fractures as fullyfilled veins and faults with clay gouge, highlighting the importance of distinguishing between different fracture types when using a field analogue for a geothermal reservoir.

Potential for Reactivating and Generating Fractures
Results of the models can be expressed in terms of the critical values of fluid pressure or horizontal (tectonic) stresses needed to reactivate existing cohesionless fractures, or to generate new fractures ( Table 4). Note that these are simply the results of the model, which is based on poorly constrained parameters and should only be taken as suggestions of what may occur under the conditions stated.

Bunter Sandstone
Critical values predicted by the modelling for the Bunter Sandstone at a depth of 1 km are presented in Table 4, from which the following predictions can be made. (1) Favourablyorientated cohesionless fractures in the Bunter Sandstone may be critically stressed without fluid overpressure or applied tectonic stresses (i.e., the base-case model) if the geostatic stress ratio is low (e.g., 0.2). (2) Shear may therefore occur on favourably-orientated cohesionless fractures (e.g., joints) if the fluid is hydrostatically pressured (i.e., fluid pressurẽ 10 MPa; e.g., Figure 3b)

Bunter Sandstone
Critical values predicted by the modelling for the Bunter Sandstone at a depth of 1 km are presented in Table 4, from which the following predictions can be made. (1) Favourably-orientated cohesionless fractures in the Bunter Sandstone may be critically stressed without fluid overpressure or applied tectonic stresses (i.e., the base-case model) if the geostatic stress ratio is low (e.g., 0.2). (2) Shear may therefore occur on favourablyorientated cohesionless fractures (e.g., joints) if the fluid is hydrostatically pressured (i.e., fluid pressure ~10 MPa; e.g., Figure 3b)

Devonian and Carboniferous Greywackes and Slates
Critical values predicted by the modelling for the Devonian and Carboniferous greywackes and slates at depths of 2 km and 4.5 km are presented in Table 4. These depths are used because they represent the proposed depths for an initial research well and an exploitation well, respectively. The results of the modelling enable the following predictions to be made.
A fluid pressure of about 50 MPa would be needed to reactivate cohesionless fractures in both the greywackes and slates at a depth of 2 km if the geostatic stress ratio is high (e.g., 0.41; e.g., Figure 3b), with fluid pressures reaching lithostatic pressures; 3.
Gently-dipping extension fractures will start to develop in the greywackes at a depth of 2 km if the fluid pressure is about 72 MPa, but may develop in the slates at slightly lower pressures (about 68 MPa). The models predict that, in the absence of cohesionless fractures, increasing fluid pressure will initially create gently-dipping extension fractures. This is because of the assumption of uniaxial strain (i.e., that the rocks are laterally confined). Higher fluid pressures will be required to generate steeply-dipping extension fractures, in the absence of horizontal tensile stresses, such as those related to tectonic forces or to cooling (Section 6); 4.
Steeply-dipping extension fractures are predicted to develop in the greywackes at 2 km depth if there is a tensile tectonic stress between about −24 and −40 MPa, but are likely to develop in the slates at tensile tectonic stress between about −24 to −33 MPa (e.g., Figure 3d); 5.
Extension fractures in all directions may develop in the greywackes at 2 km depth if the fluid pressure is between about 100 and 210 MPa, but would develop in the slates if the fluid pressure is between about 90 and 105 MPa (e.g., Figure 3e); 6.
Shear (normal faulting) may begin on favourably-orientated cohesionless fractures in the greywackes at a depth of 2 km if there is a tensile tectonic stress between about 0 and −7.5 MPa, and in the slates if there is a tensile tectonic stress between about −2 and −6.5 MPa (e.g., Figure 3b); 7.
Favourably-orientated cohesionless fractures may be reactivated as thrusts in the greywackes at 2 km depth if there is a compressive tectonic stress of about 170 MPa, and in the slates at about 140 MPa (e.g., Figure 3f); 8.
New thrusts will begin to develop in the greywackes at a depth of 2 km if there is a compressive tectonic stress of about 360 MPa, but will develop in the slates if there is a compressive tectonic stress of about 270 MPa (e.g., Figure 3g).
The different critical values for the greywackes and slates predicted by the modelling are caused by the different physical properties that have been used for these rock types. The different geostatic stress ratios and failure criteria (especially cohesion and tensile strengths) are particularly important. The different geostatic stress ratios of different rock types may explain variations of fractures with lithology. This is commonly seen in greywacke/slate and sandstone/mudstone sequences (e.g., [96]). The tensile strengths of slates tend to be lower than that of greywackes, so extension fractures may develop in the slates at lower fluid pressures and lower tensile tectonic stresses than in the greywackes. Although this suggests that fractures will be created in slates before fracturing occurs in greywackes during stimulation (Section 6.1), it does not necessarily mean that there will be more, wider or better-connected fractures in the slates than in the greywackes. The model cannot predict the fluid flow properties of the slates or greywackes.

Possible Effects of Exhumation on the Bunter Sandstone
Here, we use simple modelling to determine at what depth joints are likely to have developed in the Bunter Sandstone of the Leinetal Graben, based on the changes in stresses and temperatures the rocks are likely to have experienced during exhumation. The Bunter Sandstone is currently exposed on the flanks of the Leinetal Graben but at depths of between ≈300 m and ≈850 m in Borehole Sudhein II in the Graben [37]. Note that the top of the Bunter Sandstone has been faulted out and the base is mixed with Zechstein Salt in Borehole Sudhein II [37]. The values used in the modelling are shown in Table 5, with the Mohr diagrams shown in Figure 6. The vertical stress changes as the thickness and, therefore, weight of the overlying material changes during burial or exhumation, and this causes a change in the horizontal stresses that is proportional to the geostatic stress ratio (Figure 4b). A change in depth typically causes a change in temperature that is related to the geothermal gradient (e.g., [97]). A decrease in temperature will tend to cause contraction of rock, determined by the coefficient of thermal expansion, and this will tend to cause a reduction in compressional stresses in the rock (e.g., [98]).

Possible Effects of Exhumation on the Bunter Sandstone
Here, we use simple modelling to determine at what depth joints are likely to hav developed in the Bunter Sandstone of the Leinetal Graben, based on the changes i stresses and temperatures the rocks are likely to have experienced during exhumation The Bunter Sandstone is currently exposed on the flanks of the Leinetal Graben but a depths of between ≈300 m and ≈850 m in Borehole Sudhein II in the Graben [37]. Note tha the top of the Bunter Sandstone has been faulted out and the base is mixed with Zechstei Salt in Borehole Sudhein II [37]. The values used in the modelling are shown in Table 5 with the Mohr diagrams shown in Figure 6. The vertical stress changes as the thicknes and, therefore, weight of the overlying material changes during burial or exhumation, an this causes a change in the horizontal stresses that is proportional to the geostatic stres ratio (Figure 4b). A change in depth typically causes a change in temperature that is re lated to the geothermal gradient (e.g., [97]). A decrease in temperature will tend to caus contraction of rock, determined by the coefficient of thermal expansion, and this will ten to cause a reduction in compressional stresses in the rock (e.g., [98]).  The model presented in Figure 6 includes reductions in horizontal stresses caused by decreases in overburden and temperature to predict at what depth extension fractures (e.g., joints) may develop in the Bunter Sandstone. We use a geothermal gradient of 30 • C per km, which is similar to that modelled by [99]. The model assumes, for simplicity, an initial 1.4 km depth for the top of the Bunter Sandstone, and an initial depth of 2 km for the base of the Bunter Sandstone (Figure 6a). Vertical and horizontal stresses reduce during exhumation, with the model predicting the development of joints beginning after 0.9 km exhumation, with the top of the Bunter at a depth of 0.5 km (Figure 6b). Note, however, that this base case model assumes hydrostatically-pressured pore fluids and no applied tectonic stresses. Joints would form at greater depths if the fluids are over-pressured, if there is a tensile tectonic stress, or if the geothermal gradient were higher. Here, we consider the possible effects of stimulation by increasing fluid pressure on the folded, cleaved and veined Devonian and Carboniferous greywackes and slates, based on the modelling discussed in Section 5 and using a simple schematic figure for these rocks (Figure 7a). The modelling suggests the following sequence of events as fluid pressure is gradually increased (Figure 7):

1.
Favourably-orientated cohesionless fractures (joints) may be critically-stressed in greywackes with a very low geostatic stress ratio, so may undergo normal faulting even at hydrostatic fluid pressure (Figure 7a). Extension fractures with any orientation may develop in the greywackes if the fluid pressure is between about 100 and 210 MPa (Figure 6f).

Possible Effects of Stimulation on Devonian and Carboniferous Rocks at 4.5 km Depth
The critical values for deformation predicted by the modelling for the greywackes and slates at a depth of 4.5 km are shown in Table 4. These critical values are higher than predicted for a depth of 2 km, but suggest that stimulation would have similar effects as at a depth of 2 km. One notable difference is that higher differential stresses at depth would mean that stimulation is more likely to reactivate fractures in shear, or to create shear fractures, at greater depths. The modelling predicts that the slates will be particularly likely to develop shear fractures during stimulation. Figure 8 highlights the effects of the different fluid pressures needed to stimulate reservoirs at different depths. The difference in fluid pressure between the top and bottom of a 2500 m column of water would be approximately 25 MPa. If fluid pressure is increased at a depth of 4.5 km such that it exceeds the horizontal stress, this may cause the fluid pressure at 2 km to exceed the vertical stress, if these different depths are hydraulicallyconnected ( Figure 8). This suggests that care will be needed during stimulation if it is not intended to fracture rocks at shallower levels. We note that the Zechstein salts are likely to act as an effective top-seal (e.g., [101]) during stimulation.

Possibility of Open Fractures in the Devonian and Carboniferous Rocks below Göttingen
Although we have carried out simple modelling of joint development in the Bunter Sandstone, we have not performed so for the Devonian and Carboniferous rocks that are inferred to occur below about 1.5 km under Göttingen. This is because we consider there to be too many uncertainties for meaningful analysis of the effects of Late Cretaceous to Tertiary exhumation and cooling on these deeper and older rocks. There are, however, two arguments for the occurrence of open fractures in the Devonian and Carboniferous rocks. Firstly, joints will have developed in those rocks as they were exhumed at the end of the Variscan Orogeny, before deposition of the Permian and Mesozoic rocks unconformably above. It is possible, however, that these late-or post-Variscan joints will have subsequently been mineralised. Secondly, there is evidence that joints and other open fractures can occur at the depths of the proposed geothermal reservoir at Göttingen. For example, joints have been reported at depths of >2 km in tunnels (e.g., [102]) and in mines (e.g., [103,104]). Similarly, open fractures have been reported in fractured "basement" petroleum fields at depths of >4 km (e.g., [105]).    Table 4. The fold has an amplitude of a few metres to tens of metres. The model is for a depth of 2 km, at which hydrostatic fluid pressure would be 19.62 MPa. Veins are shown in the greywackes in the outer arcs of folds, and a fanning axial planar cleavage is shown. Joints are drawn approximately perpendicular to bedding. We predict the following sequence of events as fluid pressure is increased. (a) Favourably-orientated cohesionless fractures (joints) may be critically-stressed in greywackes with very low geostatic stress ratios, so may undergo normal faulting even at hydrostatic fluid pressure.  Table 4. The fold has an amplitude of a few metres to tens of metres. The model is for a depth of 2 km, at which hydrostatic fluid pressure would be 19.62 MPa. Veins are shown in the greywackes in the outer arcs of folds, and a fanning axial planar cleavage is shown (e.g., Stephan et al., 2016). Joints are drawn approximately perpendicular to bedding. We predict the following sequence of events as fluid pressure is increased. (a) Favourably-orientated cohesionless fractures (joints) may be critically-stressed in greywackes with very low geostatic stress ratios, so may undergo normal faulting even at hydrostatic fluid pressure.

Potential Benefits, Problems, and Improvements
Although the modelling presented in this paper is simple, it has several uses when attempting to predict what is likely to occur in the sub-surface. Firstly, constructing the models requires consideration of the range of the conditions that are likely to occur, including the lithologies and their mechanical properties, stresses, fluid types and fluid pressures, and the presence of different types of fractures. Secondly, predictions can be made about whether natural open fractures (e.g., joints) may occur, and whether exhumation might have created joints at reservoir depths. Thirdly, it enables predictions to be made about the orientations and types of induced fractures that are likely to occur under a range of possible conditions (Figure 7). It may also enable predictions to be made about whether there will be induced seismicity during stimulation.
There are various potential problems with the simplicity of the approach presented here, but these suggest ways in which the models can be improved. For example: 1. The Mohr diagram models used give little direct information about potential fluid flow in the sub-surface. The approach could, however, be used in combination with

Potential Benefits, Problems, and Improvements
Although the modelling presented in this paper is simple, it has several uses when attempting to predict what is likely to occur in the sub-surface. Firstly, constructing the models requires consideration of the range of the conditions that are likely to occur, including the lithologies and their mechanical properties, stresses, fluid types and fluid pressures, and the presence of different types of fractures. Secondly, predictions can be made about whether natural open fractures (e.g., joints) may occur, and whether exhumation might have created joints at reservoir depths. Thirdly, it enables predictions to be made about the orientations and types of induced fractures that are likely to occur under a range of possible conditions (Figure 7). It may also enable predictions to be made about whether there will be induced seismicity during stimulation.
There are various potential problems with the simplicity of the approach presented here, but these suggest ways in which the models can be improved. For example: 1.
The Mohr diagram models used give little direct information about potential fluid flow in the sub-surface. The approach could, however, be used in combination with other modelling approaches. For example, it would be useful to compare predictions of critically stressed fractures from Mohr diagrams with distinct element analysis of fracture networks (e.g., [106,107]); 2.
The values for rock properties used are based on triaxial tests, which are probably over-estimates because small, unfractured samples are generally used (e.g., [86]). More accurate methods for estimating the material properties of rock masses are available (e.g., [108]), and these methods could be used when more detailed information becomes available about the fracture patterns in the rock mass; 3.
Similarly, we have used rock mechanical properties from the literature and have made various simplifying assumptions (e.g., no applied tectonic stresses, Biot coefficient = 1). The modelling can be improved and the assumptions properly tested as the input parameters become better constrained, for example as borehole data become available; 4.
The anisotropy of the slates has not been modelled in a sophisticated way here, and this can be improved using more detailed information about the relationships between the angle between in situ stresses and cleavage (e.g., [109]); 5.
The Mohr diagram analysis used here is two-dimensional, mainly because the magnitudes and orientations of the horizontal stresses are unknown. The analysis could be expanded to three dimensions when such information becomes available, for example from well data. Although predictions can be made about the stresses involved in the Variscan Orogeny and the formation of the Leinetal Graben, those predictions do not help with making predictions about the present-day stresses.

Conclusions
This paper shows how simple mechanical modelling, using Mohr diagrams and reasonable ranges of values for rock properties, stresses, and fluid pressures, can be used to predict fracturing in a potential geothermal reservoir. Inferences can be made about the range of structures likely to exist in the sub-surface and that may be generated during stimulation. Critical values of fluid pressure and applied tectonic stresses determine under what conditions different types and orientations of fractures are likely to occur or be generated during stimulation.
A model is presented for the development of shear and extension fractures in the Devonian and Carboniferous greywackes and slates as fluid pressure increases, as would occur during hydraulic stimulation (Figure 7). For a depth of 2 km, this involves: (a) shear on favourably-orientated cohesionless fractures in greywackes with a very low geostatic stress ratio under hydrostatic fluid pressure (  (Figure 7f). A similar sequence is predicted for a modelled depth of 4.5 km, although higher differential stresses at greater depths imply that hybrid or shear fractures are more likely to form than are extension fractures as fluid pressure increases.
The modelling addresses many key questions asked in Section 1. (1) The tensile strengths of slates tend to be lower than those of greywackes, suggesting lower fluid pressure or lower tensile tectonic stresses are needed to hydraulically fracture slates than greywackes. (2) Critical values of fluid pressure and tectonic stresses for fracturing can be predicted (Table 4). (3) Steeply-dipping cohesionless fractures (e.g., joints) are most likely to be reactivated first, as shear fractures. (4) Existing cohesionless fractures (e.g., joints) may initially reactivate in shear while cohesional fractures (e.g., veins) will tend to reactivate as extension fractures. (5) Shear fracturing generally requires the reactivation of cohesionless fractures as tensile tectonic stresses are applied, or the application of compressive tectonic stresses to generate new shear fractures. (6) Exhumation and related cooling may be responsible for the creation of joints in the Bunter Sandstone and the Variscan rocks beneath Göttingen. It is likely that joints developed in the Variscan rocks prior to deposition of the Permian rocks, but these may have subsequently been mineralised.
The relationships between stresses, fluid pressure and the types and orientations of fractures are important inputs for more detailed modelling techniques (e.g., DFN modelling [110]), and can be used to test the mechanical implications of those detailed modelling techniques. Not only does this modelling help geoscientists consider and model the ranges of mechanical properties of rock, stresses, fluid pressures and the resultant fractures that are likely to occur in the sub-surface, it encourages them to consider the ranges of key parameters and their effects. Funding: This work has received funding from the European Union's Horizon 2020 research and innovation programme (grant agreement № 792037-MEET Project).

Data Availability Statement:
No new data were created or analysed in this study. Data sharing is not applicable to this article.