Multiscale Characterisation of Fracture Patterns of a Crystalline Reservoir Analogue

: For an accurate multiscale property modelling of fractured crystalline geothermal reservoirs, an enhanced characterisation of the geometrical features and variability of the fracture network properties is an essential prerequisite. Combining regional digital elevation model analysis and local outcrop investigation, the study comprises the characterisation of the fracture pattern of a crystalline reservoir analogue in the Northern Odenwald, with LiDAR and GIS structural interpretation. This approach provides insights into the 3D architecture of the fault and fracture network, its clustering, and its connectivity. Mapped discontinuities show a homogeneous length distribution, which follows a power law with a − 2.03 scaling factor. The connectivity of the fracture network is heterogenous, due to a fault control at the hectometric scale. Clustering is marked by long sub-vertical fractures at the outcrop scale, and strongly enhance heterogeneity around weathered fracture and fault corridors. The multi-variable dataset created within this study can be used as input data for accurate discrete fracture networks and ﬂuid-ﬂow modelling of reservoirs of similar type.


Introduction
Natural permeable fracture zones are primary exploration and exploitation targets in granitic deep geothermal systems [1]. In the Northern Upper Rhine Graben, no clear quantitative overview of the structural network of these potential reservoirs is available yet. Such a dataset is crucial to estimate rock property variability and compartmentalisation by faults and fracture networks for deep geothermal and heat storage projects in crystalline units [2][3][4][5][6].
The crystalline basement of the northern Odenwald may be a potential location for the geothermal underground research laboratory GeoLaB to implement high flow rate experiments to understand, design, and increase the acceptability of Enhanced Geothermal Systems (EGS) [7]. Thus, the structural investigation of such units is essential. The hydraulic yields of geothermal reservoirs in crystalline rocks are controlled by the structural architecture of fault and fracture corridors and its associated hydrothermal alteration history controlling the present-day permeability [8][9][10]. Economically, the hydrothermalised and fractured top granitic basement is more promising due to naturally enhanced flow properties [11]. Thus, the structural architecture needs to be investigated, as the permeability of fracture zones is the highest in these features [3,9,12,13]. A quantification of the fault and fracture network, which can be achieved by analogue investigation, is required to predict such rock media behaviour.
This study aims to reach a quantitative knowledge basis of faulted crystalline reservoirs structural architecture in the Upper Rhine Graben (URG) vicinity, especially at the sub-seismic scale (objects that cannot be resolved with conventional seismic techniques). This near-surface dataset is thus of importance, as it provides the basis to estimate the input parameters of 3D fractured reservoir models. The dataset is used to build a conceptual model with sub-domain definitions and quantification of the fracture network geometry. Such a model predicts ranges of fracture density, intensity and connectivity, and improves the understanding of the interconnections between structural features at different scales.
This study identifies the scale dependence of the clustering of the fracture network around fault zones and their control by the primary fault system orientation. In conceptual models of granitic and granodioritic reservoirs, structure-dependent, asymmetrical fracturing and weathering influence the reservoir architecture and the localisation of natural drains [2,14].

Geological Context
The Odenwald Crystalline Complex is part of the Mid-German Crystalline High, located south of the Rheic suture between Armorica and Laurussia ( Figure 1a).
The western Odenwald, also called Bergsträsser Odenwald, is by definition divided into three units [15]: the Frankenstein Massif in the North-West (Unit I), the Flasergranitoid zone (Unit II), and the South Odenwald (Unit III). The Western Odenwald is separated from the Eastern Odenwald (also called Böllsteiner Odenwald) by the N010 • E striking Otzberg Fault Zone (OTZ). Unit (I) and Unit (II) are tectonically independent, as the age [16] of the Frankenstein complex differs from the age of crystalline units in the Bergsträsser granitic units, and as a Carboniferous [16] tectonic fault zone delimits the two units. The crystalline Odenwald is a deeply eroded paleo subduction-related volcanic arc of the active Armorican margin, as suggested by the age and nature of the granitic plutons in the Frankenstein massif [17]. The Variscan orogeny is subdivided into four deformation phases [18]. Compressional phases D1 and D2 during the Variscan orogeny initiated thrust sutures [18,19]. In the late Carboniferous, the area underwent two phases of deformation, namely, D3 with normal to oblique sinistral strike shearing and D4 with normal shearing. D3 has a NE-SW striking direction, and D4 strikes NNE-SSW. Furthermore, D4 expressed localised strain zones and contributed to the formation of brittle deformation structures. During D4, the emplacement of syntectonic plutons in the northern Odenwald was facilitated by a transtensive to extensive stress field [18]. This extensional regime continues during the Permian, with the development of large intra-mountainous basins filled with coal and siliciclastic deposits (Saar-Nahe, Lorraine) [20,21]. These sandstones are locally associated with basaltic and rhyolitic volcanic episodes. During the Mesozoic, the area entered a subsidence regime, with continuous sedimentary deposition from Triassic up to Early Cretaceous [22]. From the Late Cretaceous until the Early Eocene, a regional uplift caused the erosion of a large section of the sedimentary units down to the Lower Triassic [20,23,24].
The Cenozoic tectonic history of the Odenwald is marked by the URG rift opening [25][26][27][28]. The URG development is characterised by passive rifting in the alpine foreland [23,27,[29][30][31] and is associated in the Northern URG (NURG) with reactivation of faults systems of Variscan age [24,32,33]. This extensive regime is associated locally with basaltic and trachytic volcanism [25]. Granodioritic and granitic dykes intrude the diorites of the northern Frankenstein Complex [15,25]. In the southern part of the Frankenstein complex, gabbroic intrusions that emplaced at 362 ± 9 Ma [16] are the predominant rock type. Amphibolites and hornfels compose the host rock of these intrusions [16,17]. Several fault systems structure the Frankenstein Massif. The major one is the N000-N020 striking OTZ on the east, delimitating it from Unit IV. A second system of fault zone delimitates the western border of the Frankenstein Massif from the URG.
The Mainzer Berg quarry is located in the northern part of the Frankenstein Massif, between Darmstadt and Dieburg ( Figure 1a). The area is structured by different fault systems, striking N000-N020 • E and N110-130 • E. Quartz monzodiorite, Bt-Hbl granodiorite, and Bt-Hbl granite are identified in the nearby GA1 and GA2 wells [25,34,35]. Intrusive bodies inside the diorite unit are composed of granite and granodiorite. These intrusions, reported 1 km north of the quarry, are orientated N160-N170 • E. The quarry is mainly composed of granite and granodiorite, in which mafic unit remnants crop out in the NE part. Locally, a multi-kilometric fault system crosses the Messel pit with an N030-N040 • E direction. Several secondary faults, striking mainly N010-N020 • E, and secondly N040 • E, N090 • E, N110-N120 • E, and N160 • E, are exposed in the Mainzer Berg area. This study describes these fault systems and their damage zone architecture in detail.
The Mainzer Berg quarry is located in the northern part of the Frankenstein Massif, between Darmstadt and Dieburg ( Figure 1a). The area is structured by different fault systems, striking N000-N020° E and N110-130° E. Quartz monzodiorite, Bt-Hbl granodiorite, and Bt-Hbl granite are identified in the nearby GA1 and GA2 wells [25,34,35]. Intrusive bodies inside the diorite unit are composed of granite and granodiorite. These intrusions, reported 1 km north of the quarry, are orientated N160-N170° E. The quarry is mainly composed of granite and granodiorite, in which mafic unit remnants crop out in the NE part. Locally, a multi-kilometric fault system crosses the Messel pit with an N030-N040° E direction. Several secondary faults, striking mainly N010-N020° E, and secondly N040° E, N090° E, N110-N120° E, and N160° E, are exposed in the Mainzer Berg area. This study describes these fault systems and their damage zone architecture in detail.

Materials and Methods
This outcrop is a suitable analogue of crystalline geothermal reservoirs. It is constituted by granodiorite and biotite-rich granite, representing the magmatic basement of this section of the NURG [37,38]. Before the URG Cenozoic rifting, the NURG area and the Unit (I) of the Odenwald belonged to the same unit, and underwent a similar deformation history. Thus, it is expected that localities of Mainzer berg would show a similar fracture pattern as what can be encountered on URG rift shoulders. A multi-methodological and multi-scalar approach allowed us to access the fracture network geometry and dimension.

Materials and Methods
This outcrop is a suitable analogue of crystalline geothermal reservoirs. It is constituted by granodiorite and biotite-rich granite, representing the magmatic basement of this section of the NURG [37,38]. Before the URG Cenozoic rifting, the NURG area and the Unit (I) of the Odenwald belonged to the same unit, and underwent a similar deformation history. Thus, it is expected that localities of Mainzer berg would show a similar fracture pattern as what can be encountered on URG rift shoulders. A multi-methodological and multi-scalar approach allowed us to access the fracture network geometry and dimension.
Digital Elevation Model (DEM) analysis focuses on the regional-scale fault and discontinuity patterns investigation, with 25 m and 1 m resolution, respectively. The methodology for lineament analysis, followed in this study, is described in previous work [3,6,39,40]. A lineament is traced when a discontinuity is identified on the hill shaded DEM. Lineaments are depicted by X and Y coordinates of start and endpoints. From these coordinates, length and strike are computed. For poly-segmented lineaments, a mean strike value is given, calculated from the straight line separating the start and endpoint of the full lineament. Misinterpretations are avoided by simultaneous observation and digitalisation on four shade orientations: N000 • E, N045 • E, N090 • E, and N135 • E. The regional-scale investigations with DEM allow the identification of discontinuities with lengths ranging from several tens of kilometres down to several tens of meters, e.g., seismic scale in the frame of reservoir characterisation.
At the outcrop scale, a quantitative approach with remote sensing technologies such as LiDAR (Light imaging, detection and ranging) has the advantage to allow a structural analysis, thanks to high-resolution point clouds [41,42]. The reliability and usability of LiDAR-derived data for discontinuity analysis was the focus of several previous publications [43][44][45].
The structural data are approached by a classical 2D GIS interpretation combined with LiDAR point cloud acquisition for semi-automatised fracture plane recognition. The LiDAR device used is a RIEGL VG 400, equipped with a Nikon Camera for RGB acquisition. Raw LiDAR datasets have been processed using RiSCAN PRO and further treated with the open-access CloudCompare software.
The point clouds of the chosen quarry walls ( Figure 2a) were treated under the following workflow. (1) Cleaning the data using the Noise and Canupo filtering to remove vegetation and non-interpretable areas. (2) Norm calculation, conversion into dip direction.
(3) Automatic plane recognition via a RANSAC shape detection plugin [46] with settings that vary depending on the distance to the fitting plane and the distance of investigation in the point cloud. For long-range acquisitions, the distance feature range was selected up to 20 cm. For close-ups (spots A, B and C), a limit to 5 cm was fixed to avoid merging several fracture planes into one, thus reducing the plane recognition errors. Angle variation maximum was set to 20 • . (4) Manual quality checks to remove planes that do not represent natural fractures. (5) Export of the detected fracture planes properties to calculate volumetric fracture density and intensity (noted P30 and P32, respectively [47,48]).
The cluster analysis of the LiDAR extracted planes is performed using Stereonet 11 [49]. Parameters for the clustering are the trend, the plunge, the number of elements included in the cluster, and Fisher distribution statistical parameters a95, a99, and kappa [50].
The (x,y,z) points clouds from the outcrop surface were also rasterised and incorporated in QGIS to extract a hill shade image from the LiDAR data using the same orientation panel as for the DEM investigation. For each profile, two views were projected, e.g., the top view, corresponding to a horizontal plane, and the side view, corresponding to a vertical plane. Fault and fracture network features were digitised and interpreted on these images to extract properties such as fracture position, apparent fracture length, apparent fracture orientation from top-view projection, and apparent dip from the side-view shot. The classification of the different fault zone compartments follows [51,52]. The number of connections per line (noted C L [53]), areal fracture density and intensity (noted P20 and P21, respectively [47,48]) were calculated from the digitised layers. Fracture connectivity and node topology were assessed using the MATLAB Plugin FracPaQ ® [54], following the I-X-Y node classification detailed in [53,55,56]. On specific profiles (C, D, and E), spacing and coefficient of variation, noted C v , were calculated along artificial scanlines to determine if the spatial distribution of the fracture network is clustered (C v > 1) or randomly distributed (C v < 1) [48]. The scanline methodology is a classical approach to characterise the fracture variability in 1D [57,58]. Artificial scanline sampling provides an overview of the fracture clustering and spacing, independently from the scanline orientation. Spacing was calculated between each fracture crossing the scanline and expressed as distance-occurrence frequency diagrams [59]. This approach allows more insight into the clustering aspect of the fracture network and quantifies the deviation to a uniformly distributed fracture pattern [59]. Cumulative length distributions of discontinuities were then analysed and fitted to a power law. On a local scale, log-normal and exponential fits may also be valid [60][61][62][63]. However, the multiscale character of the fracture network investigated here requires a power law. For crystalline rocks, this approach is commonly used [2,3,6,40,64]. Power law exponent b gives insights into the structural architecture and connectivity [40,60,65]. When b is equal to 2, the analysed property has a fractal character. When b < 2, the length population is mainly composed of long fractures, whereas when b > 2, the fracture network comprises short size elements [64]. Fracture network characterisation was also amended by field observations, such as the fracture aperture measurements approach [58], and visual characterisation of fracture mineralisation.
This multi-method approach allows an extensive, almost contact-free, and multiscale investigation of structural objects.

Regional Lineaments
Overall, P20 and P21 on 25 m resolution DEM reach 1.09 × 10 −6 lineament.m −1 and 1.46 × 10 −3 m·m −2 , respectively (Figure 3a, Table 1). Density and intensity increase, with P20 of 2.47 × 10 −6 lineament.m −1 and P21 of 2.16 × 10 −3 m·m −2 when the DEM resolution increases to 1 m ( Figure 3b, Table 1). The lineament analysis in the area of interest around the Mainzer Berg granodiorite (Figure 3c, Table 1) shows P20 4.7 × 10 −5 lineament.m −1 and P21 of 6.97 × 10 −3 m·m −2 . In the Messel pit, identified lineaments have P20 of 3.14 × 10 −5 lineament·m −1 and P21 of 7.06 × 10 −3 m·m −2 . On satellite imaging of the Mainzer Berg quarry, lineament density and intensity reach 1.89 × 10 −4 lineament·m −1 and 1.26× 10 −2 m·m −2 , respectively. The lineament density and intensity are higher in crystalline rocks than in sedimentary and volcanic Permian deposits. In all the observed scales, lineaments are mainly orientated N010 • E and N080 • E. In the vicinity of the Mainzer Berg quarry, the multi kilometre-scale lineaments are orientated N010 • E, N050 • E, and N100 • E-120 • E. In the Mainzer Berg quarry walls, identified discontinuities are of smaller size (<100 m). Local lineaments have a preferred orientation of N020 • E and N100 • E. The length varies between 4 and 356 m, having a mean of 66 m (interpretation of satellite imaging, Table 1). In comparison, the Messel Maar is affected by an N045 • E striking strike-slip fault and exhibits sub orthogonal lineaments and secondary faults orientated N000 • E-N010 • E and N150 • E-N170 • E as recorded in a previous study [25]. Length of lineaments in the Messel Maar range from 37 to 1918 m, with a mean length of 224 m. Table 1. Lineament analysis statistics, with dimension features, power law parameters a, b, and rsq (r 2 ), and areal fracture density and intensity.

LiDAR Imaging at the Outcrop Scale of the Structural Network
The five sections analysed in Mainzer Berg quarry exhibit the variability of the fracture pattern in each profile (A to F). A is orientated N010 and exposes several zones of high fracture intensity with 50 cm to 1 m thickness. The fractures apertures vary from 1 to 5 mm. Several of those fracture corridors present fractures filled with clay minerals. Profile B, orientated N110, crosscuts a high-intensity fracture corridor. Clays and hematite infills were identified in these fractures with apparent apertures varying from 0.5 to 8 mm. Several measured planes also exhibit cataclastic infills, especially in large opening structures. Profile C is located in the southern part of the quarry with an orientation of N100 • E and crosses a clearly defined fault zone with two fault cores of variable fracture intensity associated to clay-rich alteration. The fault zone here reaches the boundary with the overlaying Rotliegend deposits. Profile D, orientated N170 • E, presents large fractures (length above 5 m), but no clearly defined corridors. Profile E is orientated N090 • E and presents several fracture corridors with a 50 • apparent dip angle. Some of the fracture planes exhibit clay infills. However, the alteration in clays in fracture planes is not as pronounced as for profile C ( Figure 2).
LiDAR extracted datasets of the five profiles show a complex fracture pattern ( Figure 4 and Table 2). On profile A, fracture density reaches 0.02 m −1 , and fracture intensity,   On the five profiles, power law trends can be extracted (Table 2), with coefficient a varying from 32.9 to 8537.1, with an average of 911.8; and coefficient b varying from −3.8 to −2.69, with an average of −3.29.
Eleven clusters describe 55% of the fracture orientation distribution (Table 3). Kappa here represents the parameter to estimate the distribution of fracture orientation [50] and varies from 42.2 to 190.

GIS Structural Analysis of Profiles
The digitalisation of fractures and faults on the profiles adds complementary information to the fracture pattern characterisation (Figures 5 and 6, Table 4 Artificial scanlines were traced on three sections (C, D, and E) to extract P10 and Cv values from the GIS fracture maps (Table 5)

Fault Zone Domains
The Mainzer Berg southern quarry wall offers good outcrop conditions to analyse the structural architecture of one of the secondary faults affecting the granodioritic intrusions. Four different subdomains can be defined in this outcrop section (Figure 7a,b). From east to west, a damage zone with increasing fracture density towards the west. The minimal thickness of this damage zone DZ1 reaches more than 15 m. The internal fracture network is characterised by a few long sub-horizontal fractures (length > 10 m). Fracture orientations follow an N100-115 • E main strike, with a secondary direction striking N030 • E. Horizontal connectivity is mainly represented by I and X nodes (Figure 7c), which characterise a rather stochastic network with few cross-joints [66,67]. The random distribution observed on the horizontal plane is absent on the vertical plane (Figure 7d). In the vertical direction, an increased number of Y nodes is observed. The structural network presents slight compartmentalisation by the large sub-horizontal fractures, on which secondary fractures abduct with a 65 • to sub-orthogonal angle. Towards the west, the increasing amount of sub-vertical long fractures marks the transition to the next subdomain FC1. DZ1 is adjacent to the strongly deformed fault core. The thickness of this first fault core FC1 varies between 30 cm and 80 cm, and the whole entity is orientated N305-315 • E. Vertically, sub-vertical fractures of more than 5 m in length delimit the borders of FC1. In between these long delimitating structures, most of the fracture network is composed of more minor fractures (length < 1 m) crosscutting the long vertical fractures, with a 40 • to sub-orthogonal angle. The fractures composing this fault core are challenging to identify by the LiDAR workflow due to their small size. On the horizontal plane (Figure 7c), I and Y node networks mainly control the connectivity, typical for fault systems [52,55]. Vertical connectivity is primarily present as Y nodes, resulting from the geometrical relationship between compartmentalising long fractures and smaller ones (Figure 7d). At the western limit of FC1, fracture intensity reduces, changing to another subdomain of the fault system. An 8 m thick second damage zone (DZ2) constitutes the third subdomain (Figure 7a,b). The structure of DZ2 is described by the presence of two sets of long sub-vertical fractures delimiting fracture corridors. The fracture corridors strike N310-320 • E, sub-parallel to FC1 main strike. A secondary set, with strike varying from N135 • E to N160 • E, with lower dip of 55-65 • E and associated conjugate fractures with apparent dip toward the west, crosscuts the large fracture corridors. The substantial reduction of Y nodes number, associated with a higher number of I and X nodes on the horizontal map, marks the transition to FC1 (Figure 7c). This reduction characterises a switch to a more stochastic topology of the facture network. Vertical connectivity also shows a high amount of I nodes (Figure 7d), suggesting a lower connection rate between fractures in DZ2 compared to DZ1 and FC1.
DZ2 is bordered on the west by a second fault core FC2, striking N010-N020 • E. It comprises secondary fractures, mainly striking N300 • E, delimited by long fractures, striking N010 and N020 • E. Despite their clear presence on GIS imaging, these major fractures do not have an intense relief and are not detected by the LiDAR. The average fracture length is 1.5 m, and sub-vertical fractures, extending up to above 5 m, structure FC2. Horizontal and vertical node networks in FC2 show fewer Y nodes than in FC1, suggesting lower connectivity between fractures in this fault core.

Clustering of the Fracture Network
As mentioned before, the fracture network in the Mainzer Berg quarry southern wall is clustered by clearly defined subdomains (Figure 7). Relative spacing analysis on profiles C, D, and E (Figure 8a,b, Table 5) show that this clustering is also present on profiles that do not exhibit a clearly defined fault system. Clustering of the fracture network exists, even in the background fracture pattern of the granodioritic pluton (profiles D and E, Figure 8b). Clustering here depends on the direction of analysis, as the C v value is, on average, close to 1 (with variation from 0.54 to 1.89) (Figure 8). P10 values computed from the artificial scanlines on profiles C, D, and E (Figure 5d-f) show two configurations (Figure 8b). On profile C, P10 values are high, and the coefficient of variation is varying between 0.65 and 1.4, typical of a highly fractured media with a relatively uniform distribution of the fracture network. In contrast, profiles D and E have a smaller P10, but the high C v suggests a stronger fracture network compartmentalisation. Spacing expressed as occurrence along with the relative distance [59] allows us to detect potential clustering of the fracture distribution ( Figure 9). The curves that express the most substantial clustering have the highest deviation to the random distribution (black lines in Figure 9). This enhanced clustering is specific to few directions on the three profiles. For profile C, as illustrated by Figures 2 and 7, clustering is the strongest for scanlines orthogonal to the several observed fault cores (lines 5, 6, 9 on vertical view, Figure 9a, and lines 9, 10, 16, 17 on horizontal view, Figure 9b). On the other hand, scanlines sub-parallel to the strike of the fault core have a rather random distribution pattern (lines 7 and 10, Figure 9a). Parallel to the central axis of fault cores, the fracture network is randomly distributed in spacing.  On the two other profiles, D and E, clustering is identified on lines 1, 4, and 9 on the side view (Figure 9c), and lines 6, 9, and 10 on the top view ( Figure 9d). However, here, the lower fracture density nuances the real impact of this clustering on flow anisotropy. On the vertical projection of profile E, clustering is the highest in scanlines crosscutting long fracture planes (lines 5, 8, 9, 14, 15, 16; Figure 9e). Clustering also increases towards the eastern side of the profile, as observed on lines 2 and 5, which show the highest horizontal clustering (Figure 9f).

Topology of the Fracture Network, Influence on Connectivity and Flow Properties
The analysis of the fracture network topology is helpful to characterise relationships between fracture and their mechanical interactions [56].
IXY topology from the fracture network is interpreted from the digitised profiles ( Figure 10), and the number of connections per line (noted C L [53]) is calculated from this topological analysis. The horizontal connectivity of lineaments is estimated on the three digitised layers.
All have a C L under 2, suggesting poor connectivity.
However, the effect of censoring could explain these restrained connectivity patterns, as the shape of the scanned area has a limited (x,y) extension, and the length of several fractures is equal to the extension of the investigation area. For areas where the investigation surface is far more extensive than the size of investigated elements (example of profile D), C L increases to 3.58, reaching a highly connected network level.
By comparing horizontal C L and vertical C L , obtained from top and side views of the profiles, respectively, the difference between connectivity values highlights its anisotropic character.
The anisotropic character of fracture connectivity is low. For profiles A and B, the connectivity is relatively low (CL < 2) and slightly higher in the horizontal plane. For these two small size profiles, which height is below 5 m, the connectivity is lower than for profiles C, D, and E, which have a broader angle of observation and consider larger size objects.
Profile C has a C L between 2 and 3.57, characterising a relatively well-connected fracture network. The horizontal connectivity extracted from the top view is characterised by the prominence of Y nodes, whereas a higher number of X nodes depicts the vertical topology of the network. In the damage zones of the clearly defined fault zone (profile C, Figure 7), the angle between abutting minor fractures and the corridor bounding long fractures is relatively high (above 60 • ), suggesting clustering of the mechanical deformation. Inside fault cores, however, the angle between the different families of fractures is lower than 45 • . This angle difference can be associated with a synthetic shear component, as observed for the Messel Fault zone [25,34] and associated minor shear faults.
Profile D has the highest observed C L , with vertical C L around 2, and horizontal C L above 3.57, suggesting a very well-connected network in the horizontal direction. C L is higher on horizontal than on vertical planes on profile D. Profile D also exhibits a high amount of X nodes, suggesting a different spatial organisation of the fracture topology than the other profiles. A high amount of very long and steeply dipping fractures intersect each other and favour these X nodes.
The opposite trend with higher C L is observed for profile E on the vertical plane (C L around 2) than on the horizontal plane (C L largely below 2). The higher Y nodes number suggests a topology on the vertical planes, where certain fractures abduct others, consistent with the observed corridor geometries. These corridors are sub-parallel to each other and, thus, not allowing high horizontal connectivity. Long, 45-50 • dipping fractures, associated with an enhanced clustering, favour vertical planes connectivity and accentuate the connectivity anisotropy.

Multiscale Laws of Fault and Fracture Network in Granitic and Granodioritic Units
The multiscale character of the discontinuity network is demonstrated here on five orders of magnitude from the DEM, LiDAR, and 2D GIS structural analysis ( Figure 11). The cumulative distribution of element lengths shows a good fit with a typical and continuous power law behaviour on 2D datasets, with an exponent of −2.03. The 3D datasets from the LiDAR have a similar slope. However, the LiDAR curve does not fit the 2D trend due to the higher intensity of elements detected. This artificial shift created by the 3D calculation of fracture plane dimensions directly affects the fracture intensity, over-estimated by the LiDAR compared to intensities measured on 2D GIS. Concerning lineament and fracture orientation, there is a global shift of the main directions of deformation depending on the scale of observation. Figure 11. Multiscale length distribution of lineaments in the northern Odenwald and associated strike changes. As LiDAR is a 3D dataset, they do not fit the general trend of the other 2D datasets. At the regional scale, primary directions are N045 • E, N080-090 • E, and N150 • E. These directions are associated with the boundaries between the different Northern Odenwald units [15] and lineaments of the Variscan orogeny in general [3,6,68]. N045 • E, N010 • E, N090 • E and N130 • E lineaments structure the Messel system [25]. NE striking elements were already reactivated in the transtensive shearing phase D3 [15,18,34]. The EW direction is specific to Messel, as it is minor in 25 m and 1 m resolution regional DEMs and almost unrepresented in local 1 m DEMs. This EW direction is compatible with antithetic conjugated structures. In the Mainzer Berg vicinity, main directions are N015-030 • E, N060 • E, N130 • E, and N150 • E. N015-N030 • E is the primary deformation direction, indicating two possible origins for these deformation features: (1) N015-030 • E structures to the shear direction of the Otzberg fault zone (N010 • E), and (2) influence of R (N060 • E) and R' (N015-N030 • E) conjugated structures to the Messel strike-slip (N045 • E) system. Directions N130 • E and N150 • E may be associated with the role of conjugated structures of the URG shearing zone [25], which were active since the Permian (faulting and weathering [69]). These directions are also present at the regional scale, linked with the shearing phase D4 [15,18], and thus could play a role in the orientation of favourable for permeable fault zones and fracture corridors. This trend could also extend to structural features of the buried reservoirs in the Upper Rhine Graben as in the SW-NE striking transfer zone [33].
Direction N060 • E also constitutes one of the outcrop scale main orientations, as identified on both LiDAR and GIS datasets. However, N060 • E is minor at the hectometric and kilometric scales. This orientation variability could suggest that small scale features striking EW can be decomposed into a sub-structure striking N060-N070 • E and that the N045 • E direction is only somewhat restricted to clearly defined fault systems such as in Messel (reactivated during Cenozoic) [25].
The observed multi-scalar behaviour of length distribution is typical for crystalline lithologies [3,67].
For DEM lineament topological analysis, a shift is observed from a Y-node prominent geometry to an X node geometry, with the resolution increase from 25 m to 1 m ( Figure 10). Y-node geometry is typical for fault-related lineaments with a clustering behaviour [48,52,56]. Thus, elements identified at this scale may be considered as tectonic structures. The topology extracted from the network is dependent on the resolution of the input data, regardless of the size of the investigated area. Thus, there is no clear multiscale continuity in the behaviour of discontinuity topology.

Conceptual Granitic Reservoir Model from the Analogue
The investigation of the Mainzer Berg quarry suggests a structural architecture of the granodioritic reservoir, which is comparable to previous studies on crystalline basement [14,67,70,71].
The structural alignment of crystalline lithological boundaries controlled by granitic plutons and associated faults constitutes the first level of heterogeneity.
In the granitic pluton, the pluton magmatic fabric can control the structural architecture [72][73][74] and by fault systems and their associated clustering. The background fractures show strikes that follow the regional trend and a relatively random spacing. In fault zones, the fracture network is clustered into fracture corridors. The transition between these fracture geometries is marked by the modification of the fracture network topology, which compartmentalises the reservoir with an increased role of metric-size fractures. This compartmentalisation is supported by long sub-vertical fractures delimitating fault corridors in the fault zone vicinity.
These large-sized fracture corridors and secondary fault systems can reach the sedimentary cover, enhancing mineral alteration in high intensity subdomains (FC1, profile C). Such configurations are typical for permeable fracture zones targeted in deep geothermal systems [1,14]. The understanding of the spatial configuration of these features is thus essential for crystalline reservoir exploration. From this analogue study, the fracture net-work connectivity and topology are closely linked to the architecture of corridors and the intensity of deformation concentration within these corridors.
The transition to the buried geothermal reservoir structural analysis reveals that N-S, NW-SE fault systems are major. In contrast, NE-SW systems have secondary importance on localisation of dilation and are thus more eager to present permeable fracture corridors [5]. With the orientation of the present in situ stress field with the maximum horizontal stress being NNW-SSE orientated, the fracture and fault systems are preferably orientated for shear and tensional reactivation and thus act as permeable corridors. Mapped lineaments intensities are higher in granitic and granodioritic plutons than in sedimentary and paragenetic metamorphic domains. Similar variability is reported in the crystalline basement regionally [3,6]. It can be explained by (1) a more substantial relief in granitic domains than in the Sprendlinger Horst area and (2) shallower dip angles of structures in sedimentary and gneisses than in granitic and granodioritic units. Because of more sub-vertical features present in batholiths, these items are more easily identified in lineament mapping [6].

Conclusions
The Mainzer Berg quarry offers a 3D live overview of the structural architecture of a granodioritic reservoir analogue, exposing a heterogenous fracture pattern with clearly compartmentalised faults and fracture corridors, as well as less deformed areas. This study provides a large dataset of quantitative estimation of fracture network dimensions and architecture related to large-scale structural features. The need to quantify the fracture pattern parameters implies the statistical estimation of the orientation, length, and spacing of fractures and faults on five orders of magnitude, which are provided here. A global power law with an exponent of −2.03 fits the length distribution. Thus, it could be considered a homogeneous value for the different potential discrete fracture network modelling scales.
However, orientation variation is scale-dependent. At the regional scale, major elements structuring the lineament network are the Messel (N045 • E) and Otzberg (N010 • E) faults, and secondly, the URG boundary fault, striking NS. The lithological boundaries and associated fault systems orientated ENE-WSW also describe the structural architecture. These regional directions are not expressed with the same proportion at all scales, suggesting compartmentalisation of the deformation, starting from the hectometric scale down to the metric scale. Secondly, the heterogeneity of the fracture pattern at the outcrop scale and the variable clustering of the fracture network should be integrated into modelling approaches to quantify flow properties in such lithologies and structural contexts. Such scale dependency of specific parameters controlling the connectivity and the clustering of lineaments and fractures and thus also the hydraulic and mechanical reservoir properties are demonstrated herein. It has to be considered for the multiscale characterisation of implementation sites, heat storage and deep geothermal reservoirs, or the GeoLAB tunnelling project.