Methods and uncertainty estimations of 3-D structural modelling in crystalline rocks: a case study

Exhumed basement rocks are often dissected by faults, the latter controlling physical parameters such as rock strength, porosity, or permeability. Knowledge on the threedimensional (3-D) geometry of the fault pattern and its continuation with depth is therefore of paramount importance for applied geology projects (e.g. tunnelling, nuclear waste disposal) in crystalline bedrock. The central Aar massif (Central Switzerland) serves as a study area where we investigate the 3-D geometry of the Alpine fault pattern by means of both surface (fieldwork and remote sensing) and underground ground (mapping of the Grimsel Test Site) information. The fault zone pattern consists of planar steep major faults (kilometre scale) interconnected with secondary relay faults (hectometre scale). Starting with surface data, we present a workflow for structural 3-D modelling of the primary faults based on a comparison of three extrapolation approaches based on (a) field data, (b) Delaunay triangulation, and (c) a best-fitting moment of inertia analysis. The quality of these surface-data-based 3-D models is then tested with respect to the fit of the predictions with the underground appearance of faults. All three extrapolation approaches result in a close fit (> 10 %) when compared with underground rock laboratory mapping. Subsequently, we performed a statistical interpolation based on Bayesian inference in order to validate and further constrain the uncertainty of the extrapolation approaches. This comparison indicates that fieldwork at the surface is key for accurately constraining the geometry of the fault pattern and enabling a proper extrapolation of major faults towards depth. Considerable uncertainties, however, persist with respect to smaller-sized secondary structures because of their limited spatial extensions and unknown reoccurrence intervals.


Introduction
Geological information is inherently three-dimensional (3-D) in space but often represented in 2-D (Jones et al., 2009).With increasingly available computer power, 3-D modelling or geometrical visualizations have become widespread, as they can be performed on a desktop computer (e.g.Bistacchi et al., 2008;Caumon et al., 2009;Hassen et al., 2016;Sausse et al., 2010;Stephens et al., 2015).3-D models widely serve as a basis for subsequent investigations, such as stress modelling or fluid flow modelling (e.g.Hassen et al., 2016;Stephens et al., 2015).Explicit structural modelling can further be subdivided into stochastic and deterministic methods.Deterministic approaches yield a single output for input parameters, analogous to drawing a map (e.g.Stephens et al., 2015), whereas as in stochastic approaches parameters are defined by a probabilistic density function with a component of randomness (e.g.Cherpeau and Caumon, 2015;González-Garcia and Jessell, 2016;Jørgensen et al., 2015;Koike et al., 2015).
When modelling a certain volume of Earth's intermediate deep subsurface (tens of metres to kilometres), as is often done for planning nuclear waste repositories, geothermal projects, or tunnelling work, 3-D structural modelling commonly starts from a known lithological and structural dataset from the Earth's surface or underground facilities such as tunnels or boreholes.Known information is then extrapolated towards the unknown.At the time of extrapolation, the validity cannot be proven unless additional information, such as geophysical, borehole, or excavation data, is integrated.
Previous studies report that this extrapolation represents a main uncertainty within 3-D structural modelling of known structures (e.g.Baumberger, 2015;Bistacchi et al., 2008).From environments with sparse data, for example, the topology of the fault network is known to be highly uncertain or prone to the existence of unknown faults (e.g.Cherpeau et al., 2012;Cherpeau and Caumon, 2015;Hollund et al., 2002).
More generally, for kilometre-scale models, uncertainties in accuracy related to input data (i.e.GPS location, dip-dip azimuth measurements) are small compared to the uncertainty related to the data interpolation between known locations or data extrapolation (Bond, 2015).
Uncertainties play an important role when considering decision-making based on information available from a 3-D model and have therefore been subject to extensive studies in the past (e.g.Bistacchi et al., 2008;Bond et al., 2007a;Lindsay et al., 2012;Tacher et al., 2006;Wellmann et al., 2010Wellmann et al., , 2014;;Wellmann and Regenauer-Lieb, 2012;Yamamoto et al., 2014).Since models are a function of the data used, some of the approaches tend to analyse uncertainties in the input data before modelling (e.g.Bond et al., 2007b;Jones et al., 2009).Other approaches investigate the error propagation into the models, inferring the uncertainty after modelling (Jessell et al., 2010;Lindsay et al., 2012;Viard et al., 2011;Wellmann et al., 2010).Most of these published studies were performed within sedimentary environments where parameters such as stratigraphy, layer thickness, layer orientation, and structural setting are well constrained.Uncertainty estimation and its potential reduction are less well constrained for the structural modelling of basement rocks (e.g.Svensk Kärnbränslehantering AB, 2009), which are characterized by intrusive contacts and a complex arrangement of deformation structures.
In this study, we focus on deformed basement rocks and the extrapolation of faults to depth.We follow two main goals: (i) the application of an extrapolation workflow for three different techniques for projection of surface structures to depth considering associated projection uncertainties and (ii) the design and application of a probabilistic approach to compare different extrapolation techniques in order to validate the generated models.
We focus specifically on the combination of observations in outcrops at the surface with observations in an underground facility, allowing for an extrapolation modelling approach, and propose that it is possible to link these two types of observations in a probabilistic context by taking into account uncertainties in measurements and the exact tie between observed features at the surface and in the underground facility.We investigate a local case study in a relatively simple setting in crystalline rocks.The study area is characterized by well-exposed crystalline rocks of the Aar  et al., 2017).The coordinate system of the inset map is given in WGS84, whereas the coordinates of the local geological map are given in the Swiss coordinate system (CH1903).massif in the central Swiss Alps (Fig. 1) and furthermore greatly benefits from subsurface information from the Grimsel Test Site (GTS) underground rock laboratory run by the Swiss Cooperative for Disposal of Radioactive Waste (Nagra).
This combination of good outcrop conditions at the surface and independent high-quality subsurface information allows for an extrapolation modelling approach and subsequent validation in a relatively simple and well-constrained hightopography crystalline setting.

Geological setting
The study site is located in the Haslital in the Central Alps (Switzerland; Fig. 1) within the Aar massif, an external crystalline massif in the Alps representing exhumed basement rocks of the former European continental margin and thus belonging to the palaeogeographic Helvetic domain of the Alps (e.g.Mercolli and Oberhänsli, 1988;Pfiffner, 2009;von Raumer et al., 2009).
Meta-basic dykes, formerly called lamprophyres (Oberhänsli, 1986), intrude into the granitoid bedrock without altering the granitoid, indicating only slightly younger intrusion ages of former basic dykes with respect to the calcalkaline granitoids.
Several authors have described the deformation related to Alpine orogeny in the vicinity of the study area (e.g.Baumberger, 2015;Challandes et al., 2008;Choukroune and Gapais, 1983;Goncalves et al., 2012;Keusen et al., 1989;Marquer et al., 1985;Rolland et al., 2009;Steck, 1968;Wehrens et al., 2016Wehrens et al., , 2017)).Ductile deformation is expressed by a pervasive foliation and localized high-strain zones (shear zones).The exact geometry of the 3-D shear zone network, which occurs at a variety of scales ranging from several kilometres down to millimetres, is complex (Choukroune and Gapais, 1983).It is, however, possible to extract a pattern of kilometre-long major shear zones interconnected by hectometre-long subordinate bridging structures.The major shear zones tend to be quasi-planar (Baumberger, 2015;Wehrens et al., 2017) and we therefore assume a considerably simplified shear zone pattern with quasiplanar to planar geometries of the major shear zones grouped according to their strike orientation (Fig. 2).
The kinematic framework of Alpine deformation is controversial.Several kinematic models have been proposed for the shear zone network genesis in the study area, including single-phase (Choukroune and Gapais, 1983) and multistage evolution models (Herwegh et al., 2017;Rolland et al., 2009;Steck, 1968;Wehrens et al., 2016Wehrens et al., , 2017)).This study aims to reconstruct the present day 3-D geometry.Although the kinematic evolution is beyond the scope of this study, the generated models have been validated for kinematic inconsistency with respect to the known tectonic framework (e.g.models with unrealistic dip values have been removed; dip < 60 • or north verging).The different orientations of the structures are therefore used without kinematic implications.The major orientation of structures within the area are NE-SW (group A), E-W (group B), and NW-SE (group C) trending (Schneeberger et al., 2016;Wehrens et al., 2017).
The pervasive foliation and the highly localized shear zones form mechanical anisotropies, which favour subse- quent brittle localization (Belgrano et al., 2016;Kralik et al., 1992).Deformation in the brittle regime is expressed by fracturing and cataclasis, often resulting in fault gouges (Bense et al., 2014;Wehrens et al., 2016).The spatial distribution of fractures and their reactivation in the form of fault gouge development is heterogeneous (Bossart and Mazurek, 1991;Mazurek, 2000).
Although the shear zones experienced a severe ductile deformation history, most of them were reactivated in a brittle manner during the exhumation history (Wehrens et al., 2017).Subsequently, we therefore use the term fault as a summary term for high T ductile shear zones, low T ductile shear zones, and their reactivation by brittle shearing leading to cohesive (protocataclasite, cataclasite) or non-cohesive (fault gouge) fault rocks.
Present day seismic activity (Pfiffner and Deichmann, 2014) indicates ongoing recent tectonic activity in the deep subsurface of the Aar massif.
Glaciation and glacial retreat contributed to the latest history of the area (Wirsig et al., 2016).Basal erosion and the latest young (17.7 ka, Wirsig et al., 2016) retreat ages produced excellent outcrop conditions, as most outcrops are glacially polished and above the treeline, exposing bare bedrock.
Owing to deglaciation, exfoliation jointing occurred (Ziegler et al., 2013).Given the restricted near-surface occurrence of these exfoliation joints and their small dimensions, we exclude these deformation features from further consideration in this study.3 Methods

Extrapolation workflow
In order to represent the 3-D geometry of faults, we developed a workflow based on a combination of remote sensing and fieldwork (Fig. 3).As a first step, we generated a lineament map using remote sensing data.We use the term lineament as defined by Gabrielsen andBraathen (2014) andO'Leary et al. (1976): a lineament is a mappable linear or curvilinear feature identified by remote sensing, possibly representing the intersection between a planar to subplanar structural anisotropy and the Earth's surface.Lineament mapping followed the methodology presented by Baumberger (2015).Aerial photographs (swisstopo) and a digital elevation model (DEM; swisstopo) with resolutions of 0.5 and 2 m, respectively, served as a basis.
Using the DEM, hillshade images (i.e.greyscale relief images) with distinct illumination angles (0-360 • illumination azimuth with 45 • steps constant at a 30 • altitude angle) were calculated, resulting in eight hillshade images and illuminating different parts of the investigation area.On a pixel-based map, the possible strike angle of a line depends on the number of pixels of the raster matrix in which the line is enclosed (Heilbronner and Barrett, 2014).Our approach requires an angular resolution < 10 • , and thus a minimum length of 10 pixels for a specific lineament was necessary to fulfil this criterion.Hence, shorter lineaments (< 5 m) were discarded.Lineaments were manually digitized and are composed of a minimum of two endpoints and potentially several points in between.
The strike of lineaments was defined as the angle measured clockwise from north.Two different approaches to analyse the strike of lineaments were compared: (i) single strike values from endpoint to endpoint and (ii) strike values for in-dividual segments between a lineament's nodes.In both approaches, a weight is added to the strike proportional to the length of the lineament.
In addition to the aforementioned remote sensing approach, conventional structural surface mapping over an area of 13 km 2 was performed.Spatially restricted outcrop observations at the surface were extrapolated along strike using the lineament map; thus combining fieldwork and remote sensing allowed us to obtain a structural surface map.Ductile deformation was mapped by differentiating pervasive background strain and localized high-strain zones (shear zones).At the surface, mapping of brittle deformation focused on the occurrence of fault gouges.In addition, mapping in the GTS underground facility was performed similarly to surface mapping on a decametre scale and in more detail regarding brittle structures (Schneeberger et al., 2016).
Structural modelling was performed using Move ™ software (Midland Valley) on two distinct scales: a local scale (decametre) for the GTS and a regional scale (kilometre) for the entire study area.Underground 3-D structural modelling was performed on the basis of underground mapping and drill core data, which resulted in fault traces and orientations.This information provided the basis for the 3-D reconstruction of fault planes.Regional 3-D structural modelling was performed following published workflows using the surface fault map as a basis (e.g.Baumberger, 2015;Bistacchi et al., 2008;Kaufmann and Martin, 2009;Zanchi et al., 2009).Surface faults were extrapolated to depth by assigning a dip value to individual surface traces, where a trace is the intersection between the Earth's surface and a fault.Three different extrapolation approaches were applied: (i) extrapolation along measured dip and dip azimuth (fieldwork-based approach).Data from outcrops were considered within an orthogonal distance of < 20 m to inferred fault traces and a strike differing less than 20 • compared to the fault's mean strike as defined by remote sensing.The fault's mean strike was calculated via linear regression through all points defining its trace.(ii) Delaunay triangulation is a meshing algorithm that produces a triangulation for several points such that for a given point cloud, no point of the point cloud is inside the circumcircle of any triangle connecting three points of the point cloud (Delaunay, 1934).It results in a 3-D surface interpolating the selected points.Based on this surface, the entire fault trace can be extrapolated.Noise can arise because of rugosities of the fault planes, uncertainties in tracing the fault intersection at the surface, and too-low vertical variations in topography.In the case of near planar faults, the noise is reduced in the case of high variations in altitude between valleys and mountain peaks and by preferring projections of long fault traces over those of short fault segments.(iii) The ribbon tool is a Move ™ internal interpolation algorithm based on a three-points approach in which three points form a triangle and the orientation is averaged over a defined number of triangles (Midland Valley).The maximum dip orientation of each average triangle is represented as a stick at the location of the starting point.The combination of all sticks along a trace results in a plane for the given trace.More details on the method used by the ribbon tool are given in Fernandez (2005) and Baumberger (2015).
For each approach, the surface fault trace was extrapolated to depth using the obtained specific orientation.Subsequently, the intersection line between the extrapolated plane and a horizontal plane at GTS elevation (approx.1730 m a.s.l.) was calculated.Then, the resulting intersection lines were compared with the underground structural map in order to find the "best-fitting" underground structure to the obtained intersection line.The degree of fit between the intersection line at the surface and the trace of the underground structure was estimated using the orthogonal distance (distance misfit) starting from the intersection with the main gallery and the angular difference (angle misfit) between the two linear features (Fig. 4).Only structures within the same orientation group (groups A, B, C) were compared.
Furthermore, the degree of fit was compared between the different extrapolation approaches and thus for every surface fault.Considering all approaches, a best-fitting underground fault was assigned based on the aforementioned criteria.This assignment served as basis for the following structural modelling step in which every surface fault was linearly interpolated with the assigned best-fitting underground fault, yielding a "best-estimate" model.

Bayesian inference
For a better description of the system taking into consideration the inherent uncertainty in the extrapolation methods In this study scenario, we assign a parameter to each surface fault at the tunnel level.We represented the uncertainty about the exact value with a Gaussian distribution and a constant standard deviation of 40 m in the horizontal axis (Fig. 5) based on the dip uncertainty of 10 • based on the dip variation in multiple orientation measurements along a single fault (Fig. 5).As a mean value, we assign the best-estimate model from the previous interpolation.Interpolated planes were grouped according to their orientation into three separate groups (identified by A, B, C in the following; Fig. 2).
Within each orientation group, we expect faults to be mostly parallel with limited intersections based on field observations.To capture this idea, we assigned a penalty factor that reduces the log-likelihood of a parameter set for an increasing number of intersections (by 0.05 per intersection, to be precise).The number of intersections per iteration was calculated using the Bentley-Ottmann algorithm (Supplement; Balaban, 1995;Bentley and Ottmann, 1979).
The described Bayesian inference cannot be performed directly due to the complexity of multiple parameters in several groups and the non-linearities due to the fault intersections.We therefore apply a computational sampling method based on an adaptive Metropolis MCMC approach (Haario et al., 2001) implemented in the probabilistic programming package PyMC 2 (Patil et al., 2010) and previously successfully used in a geological context (de la Varga and Wellmann, 2016).
Final posteriors were discretized to match the locations of measured faults in the underground tunnel by a simple nearest location classifier (Fig. 5).Therefore, the final result of the inference is a discretized distribution of each of the parameters.In order to compare the 3-D models obtained by the three extrapolations approaches, we then use the maximum a posteriori value, i.e. the highest probability value of the posterior distributions.4 Results

Lineament map
In total, 5277 lineaments with a spatially heterogeneous distribution and lengths ranging from 5 to 1941 m were mapped (Fig. 6a).Lineaments are generally more concentrated along topographic highs and lows.Within certain areas (areas i and ii in Fig. 6a), the lineament's strike tends to be parallel to the dip azimuth of the slopes, yielding uniform orientations.In contrast, in domains with relatively low topographic variations, a variety of strike orientations become discernible (area iii in Fig. 6a).
Looking at the bulk data, lineaments show a major NE-SW and a minor NW-SE trend (Fig. 6b).Long lineaments (> 400 m) are mainly oriented NE-SW (Fig. 6c), whereas short lineaments show a considerable variation in strike (Fig. 6d).The two methods for estimating lineament orientations, endpoint to endpoint (Fig. 6b) and individual segments (Fig. 6e), yield similar results.

Field observations and data
Data obtained from fieldwork combined with a compilation of several published maps (Baumberger, 2015;Keusen et al., 1989;Vouillomaz, 2009;Wehrens et al., 2017;Wicki, 2011) yielded a surface fault map (Fig. 7, see also Schneeberger et al., 2016).Based on their orientation, we discriminated different groups of faults (Fig. 7): group A are mainly steep SE-dipping faults.Their average orientation (dip azimuth to dip) is 149/74.Group A faults mostly show steeply plunging stretching lineations resulting from ductile shearing.Group A can be correlated with faults formed during the Handegg phase (22-17 Ma) as defined by Wehrens et al. (2017), while group B and group C would correspond to faults formed during the Oberaar phase (14-12 Ma; Wehrens et al., 2017).Group B are mainly steep S-dipping (mean orientation: 178/72) faults.Lastly, group C are SW-dipping faults coeval with group B with an average orientation of 196/72.Group C faults are subparallel to meta-basic dykes and often co-occur spatially with the latter.Groups B and C mostly show oblique to horizontal stretching lineations.For multiple orientation measurements along individual faults, the standard deviation of the mean dip azimuth was below 15 • and the mean dip below 10 • .Generally, the GrGrdominated southern area shows an increased number of faults (Figs. 7 and 8).Detailed underground mapping resulted in a lithological (Fig. 8a) and a structural map of the GTS (Fig. 8b).
Meta-basic dykes occur as three distinct swarms, two located within the CAGr domain (Fig. 8a).The northern two swarms strike NW-SE, whereas the southern swarm strikes E-W; however, it is less clearly marked.Numerous dykes are overprinted by an Alpine foliation, which is sometimes oblique to the dyke boundary.Furthermore, dykes are often overprinted by localized ductile and brittle deformation expressed by shear zones and fault gouges.
Faults occur along three NE-SW trending swarms, two E-W trending swarms, and two NW-SE trending swarms, leading to a heterogeneous strain distribution along the underground facility (Fig. 8b).
The NE-SW trending swarms correspond to group A faults with an average spacing of ca.16 m.In total, 31 group A faults were mapped underground.They can be further subdivided into 17 moderately to steeply dipping faults (between 45 and 75 • ) and 14 sub-vertically dipping ones (> 80 • ).The E-W trending swarms correspond to faults with orientations that are similar to group B. In total, 12 of these E-W striking faults were mapped.
The NW-SE trending fault swarms are localized mainly along dykes (Fig. 8) and represent group C structures.In total, 25 NW-SE striking faults occur within the GTS.
Faults in the CAGr (northern part) seem to preferentially localize along pre-existing anisotropies, i.e. hightemperature brittle fractures (biotite coating) or meta-basic dykes, and thus form discrete faults (centimetre sized) with marked contacts to the host rock.In contrast, faults in the GrGr-dominated southern part form strain gradients over larger distances (metres).This observation is in agreement with the findings of Wehrens et al. (2017).

3-D structural modelling
The GTS model size is 600 × 250 × 100 m, whereas the regional model size was 4×3 km with a projection depth reach-  ing the underground facility for all faults.The projection depth was defined arbitrarily but no larger than half of the fault trace's length.All 3-D models are provided in the Supplement.

GTS model
The combination of the above-presented underground map with measured surface orientation data resulted in a 3-D geometric visualization of meta-basic dykes and faults mapped underground.Swarms of meta-basic dykes tend to join towards less numerous dykes with depth.Based on geomet-rical considerations, we infer the occurrence of three major dykes from which all others either fan out or form relay structures between the major dykes.Based on the field observation that the major faults and relay structures dip steeply sub-vertically towards the south, we discriminated 8 major group A faults and 23 relay structures.Major group A faults occur within each NE-SW trending swarm discriminated on the map view.Group B deformation structures can be further subdivided into six major and seven relay faults.Group C deformation structures can be subdivided into 6 major and 32 relay deformation structures, some of which are very short (14 m).

Regional model
The surface fault map (Fig. 7) served as a basis for the generation of the three different kilometre-scale 3-D models (see above).All three modelling approaches yielded the 3-D geometrical visualization of the surface fault pattern.They all share the same fault traces at the model surface.As mentioned above, projection specific dip values were used for each of the models.However, not all surface faults were extrapolated with each approach.Of the 21 possible surface faults, 10 were extrapolated with the fieldwork-based approach, 11 using the Delaunay triangulation, and 13 with the ribbon tool method.Missing projections can be due to a lack of outcrop description or the absence of sufficient topographic relief for remote-sensing-based approaches.
By combining all three approaches, at least 1 (but up to 3) degrees of fit with underground faults were calculated for each surface fault.Based on the different degrees of fit, a best-fitting underground structure was assigned to each surface fault.By linearly interpolating the two traces, we obtained a model which we called best-estimate model.In total, 11 group A faults reach the GTS.From the total 11, 7 have a dip < 80 • , which would correspond to the major structures defined in the above-presented GTS-scale model, whereas the 4 steeper faults correspond to relay structures.Moreover, two group B and eight group C faults connect the surface with the GTS.The combination of all faults yields an average spacing of 25.4 m, and faults appear to converge with depth.

Bayesian inference
For each model that is obtained when each surface point (intersection between surface fault and 2-D section along the GTS) is interpolated with a specific underground point, the number of intersections was calculated and the likelihood of the model compiled based on the number of intersections.In total 10 000 models were calculated and for each a probability for a certain interpolation of a specific surface point with an underground point was obtained (Fig. 9).For certain surface points, a clear maximum a posteriori value was  found (Fig. 9a-d); however, for other surface points no underground point could be assigned unambiguously (Fig. 9e).
Based on the maximum a posteriori value, a 3-D structural model was obtained by linearly interpolating each surface point to the underground point with the maximum a posteriori value.We call this model the "maximum a posteriori" model (Fig. 10).Note that the maximum a posteriori model only adds information to the initial model through consideration of a likelihood, i.e. the assumption that crossing faults at a large scale are unlikely.Note that the smaller-scaled relay structures are not considered in this approach.
This maximum a posteriori interpolation model served as a basis for comparing different employed extrapolation approaches.The comparison did not yield a clear "best" extrapolation approach; however, it seems that fieldwork-based approach results in the most accurate extrapolation.

Lineament map
A comparison of the remote-sensing-based lineament map and field data showed that in intact granitic rocks, purely ductile shear zones without later brittle overprinting not detected by remote sensing.Brittle deformation generating fractures, cataclasites, or even fault gouges responsible for mechanical weakening is necessary to form morphologically detectable structures (Fig. 11a; Baumberger, 2015).Moreover, the orientations of the slopes play an important role, as faults striking in the down-dip directions of slopes are prone to the most effective erosion processes driven by gravity.Different orientations observed on the lineament map (Fig. 6a, areas i and ii) for the eastern and western flank of the Hasli valley are interpreted to result from such preferential erosion.In contrast, the surface area (iii) in Fig. 6a is nearly horizon- tal, thus reflecting a homogeneously eroded pattern of intersection for lineaments.The dependence on erosion for the formation of morphological incisions leads to the observed heterogeneous lineament density distribution as ridges and valleys show higher lineament densities.Endpoint-to-endpoint strike and the strikes of individual segments of lineaments are very similar (Fig. 6b and e), indicating only small variation in the strike of the lineaments themselves.Therefore, underlying structures should be quasilinear to linear in 2-D and planar in 3-D.We also observe that the longest lineaments are NE-SW striking and that the variability shown in Fig. 6e is mostly due to varying strike orientations of very short lineaments (< 20 m).In addition to the NE-SW striking maximum, few long lineaments strike NW-SE.Both major orientations are similar to those reported from field observations (Figs. 7 and 8) and correlate with previous studies (Rolland et al., 2009;Steck, 1968;Wehrens et al., 2017), indicating that lineament maps are suitable to obtain the general trend of steep faults in wellexposed crystalline terrain.Much care is needed, however, when further interpreting lineament maps, as the geologic meaning of the lineament is ambiguous and lineament maps are strongly operator dependent (e.g.Scheiber et al., 2015).

Field observations and data
Differences between the surface map and the underground map are relatively small.The spacing of faults at the surface is lower, but general orientations are comparable (Figs.7 and  8) and the two mappings are thus discussed jointly.
Faults commonly show little variation in orientation along strike as evidenced by consistent orientations of the dip and dip azimuth of multiple outcrop descriptions along the same fault.In conjunction with the small variability in strike for lineaments, this is clear evidence for the planarity of largescale faults.At the surface, the 2-D length of faults is between 229 and 5591 m (mean 2199±1603 m).Therefore, extrapolation of surface faults to depths similar to the depths of the underground faults, which have an overburden between 420 and 520 m, is well in the projection depth range assuming a circular shape for the plane as a minimum estimation for their lateral extent.
Localization processes seem to differ between the two host rocks (CAGr and GrGr; Wehrens et al., 2017).The higher amount of biotite in the GrGr could influence the rock's rheology towards more ductile behaviour.In contrast, the relatively higher amount of quartz and K-feldspars renders the CAGr more brittle than GrGr at similar pressure P and temperature T conditions and thus enforces brittle fracturing and possible subsequent ductile shear zone widening, as observed in other crystalline rocks (Guermani and Pennacchioni, 1998;Mancktelow and Pennacchioni, 2005;Wehrens et al., 2016Wehrens et al., , 2017)).Hence any mechanical anisotropy, such as along pre-existing structures in the form of magmatic shear zones, meta-basic dykes, or aplitic dykes, served in the CAGr as sites for strain localization when suitably oriented with respect to the stress field.

3-D structural modelling
Our 3-D structural models were generated as a contribution to a project monitoring several parameters, such as microseismicity and in situ stress conditions, on the kilometre scale (large-scale monitoring; Nagra).Therefore, 3-D structural models were required mostly for visualization purposes.A deterministic explicit modelling workflow was required, as often is used in applied projects.It is, however, clear that for model updating, an implicit modelling approach would result in faster data handling.The deterministic approach was chosen because we attempted to obtain a geometrically satisfying product within the simplest geological setting possible.Furthermore, we were interested in the actual geometry of the faults dissecting granitoid rock bodies.Lastly, the uncommonly well-constrained setting of our study site (high-resolution underground data) was used to test and potentially validate extrapolation techniques for common application.Therefore, the underground data were only as validation and not as a constraint during interpolation.

Three different approaches to obtain
extrapolation 3-D structural models (kilometre-scale models) Uncertainty related to the assignment of specific dip values to lineament traces (Baumberger, 2015;Bistacchi et al., 2008) led to the comparison of three different approaches.
Validation attempts by comparison with underground mapping are purely geometrical and were based on two criteria, namely angle and distance misfit (Fig. 4).All three extrapolation approaches yielded similar results and no significant differences were observed.Moreover, in order to allow for a thorough comparison between the different extrapolation approaches solely based on the angle and distance misfit, the underground faults would need to be homogeneously distributed, which is not the case (Fig. 8).
The validation procedure could be refined using fault core thickness.However, fault thickness varies substantially along strike (e.g.Torabi and Berg, 2011) and thus is not a clear distinction criterion.
In addition to the average dip, the maximum and minimum dips could be used, which would yield a projection cone similar to the uncertainty visualization suggested by Baumberger (2015).Applying this approach to a restricted area, such as the underground rock laboratory investigated in this study, resulted in total coverage and no possible distinction between different faults.However, for a final representation of the uncertainty related to the dip value on a regional scale (kilometre scale), the approach of visualizing projection cones would suit.

GTS (decametre-scale model) compared with kilometre-scale "best-estimate" model
As a result of differences in outcrop conditions, the number of observed faults is significantly higher in the underground laboratory compared to the surface (Fig. 12).Underground, nearly 100 % of polished outcrop is accessible along the tunnel walls (Fig. 11b), whereas at the surface faults are often covered with vegetation, even in relatively vegetation-poor domains.Furthermore, we observe convergence of surface faults with depth in our best-estimate model, which could be a modelling artefact.The N-S extent of surface area is larger compared to the GTS area, leaving a northern and a southern surface part underneath for which no underground data exist (Fig. 10).Hence faults in these domains are forced by the model set-up to be connected to the underground, leading to artificial fault orientations in these two cases (N and S rims of Fig. 13a).For that reason, only faults in the central part of the best-estimate model will be further considered.model was considered.Assuming no intersections within a large-scale fault set is simplistic, but from field observations it seems plausible (Fig. 7) as a first approach for faults belonging to a specific orientation group (group A, B, C).The maximum a posteriori model is based on a N-S cross section along the GTS, and this orientation implies that faults would only cross if their dip varies strongly.Such a strong variation in the dip value is improbable based on measured dip values (Fig. 7).Therefore, this assumption seems feasible.This simplistic representation of nature enabled us to obtain a probability for all possible interpolations between a specific surface point and all underground points of the corresponding orientation group.As previously mentioned, the margins of the interpolation space show boundary artefacts, and thus the following surface points at the model margin were not further considered: J03, J06, J20, and J21 (Fig. 8).As expected, probability densities are skewed towards the area of lesser fault density (Fig. 9).At this point it is important to remember that probability density is given as area, and therefore we cannot directly compare the discretized posteriors since they are a function of the distance between nearby faults.

"Maximum a posteriori" model
We compared the initial three extrapolation techniques based on the maximum a posteriori model (Fig. 13).When comparing group A faults, fieldwork-based extrapolation closely fit the maximum a posteriori interpolation, which indicates either that the fieldwork-based model yields the best results or validates the Bayesian inference approach depending on whether the reference state is the statistical interpolation or the measured field data.Generally, dips of the maximum a posteriori models are slightly steeper than measured dips during fieldwork (Fig. 14a).However, the dip differences between the fieldwork-based extrapolation 3-D structural model and the maximum a posteriori interpolation model are small (Fig. 14b).Also, the dip differences between the ribbon-tool-based 3-D structural model and the maximum a posteriori interpolation model are small, but dips obtained via the ribbon tool are systematically steeper, which does not correspond to the measured dips (Fig. 14a).The extrapolation 3-D structural model obtained via Delaunay triangulation is less close to the maximum a posteriori interpolation model and obtained dips vary substantially.
The comparison for the group B and C faults is less clear.Fieldwork-based and ribbon tool extrapolations are close to the maximum a posteriori model (Fig. 13).Therefore, we conclude that fieldwork is still necessary for 3-D structural modelling in crystalline environments and that the ribbon tool (Move ™ ) offers numerous options to tune the obtained plane; however, this tuning requires a profound conceptual background model.

Possible model refinements
Presented surface models include only major faults (Fig. 15).However, for further applications, such as groundwater flow modelling or slip tendency analysis, not only major faults are of interest but also their relay structures.Based on the orientation information gained from the regional kilometre-scale models and on the intersection pattern observed during lineament mapping, it is possible to infer a near surface 3-D model not only with the major fault but also the relay structure.Furthermore, the increased level of detail in the GTS model (decametre scale) forms a similar model in the underground.The unknown space between the two models would require probabilistic modelling with several key parameters, for example fault spacing, fault orientations, apertures, or cross-cutting relationships.

Conclusions
The exceptional opportunity for a surface and underground data comparison over 3-D structural modelling approaches led us to the following conclusions: Lineament maps enable the identification of major faults but are highly sensitive to preferential erosion.Structural surface mapping allowed for a discrimination of three orientation groups of faults.
Comparison based on geometrical criteria (distance and angle misfit) of the three approaches to extrapolate to depth surface traces yielded comparable results for all extrapolation approaches.
Interpolation of surface data with underground data based on a Bayesian inference problem showed that the fieldworkbased approach is the most accurate extrapolation technique.However, this could also validate the interpolation approach.
We conclude, similarly to Zanchi et al. (2009), that for 3-D structural modelling a high-topography area within crystalline bedrock, classical fieldwork as an information source and as a basis for a conceptual background model on which interpolations or extrapolations performed within 3-D structural modelling can be examined for validity.In terms of general fault networks, our approach can be applied to (i) pervasive regional fault or fracture patterns.Currently it will fail in the case of (ii) discrete large-scale faults (e.g.strike-slip faults) consisting of one fault core and an associated damage zone.In such cases, more elaborate probabilistic models have to be generated in future, including 3-D variations in terms of spacing and orientation of secondary faults and splay faults.(i) Even with limited or missing underground information, our approach can be used to predict a surface-based 2-D model including a probability evaluation (e.g.variable dip angles) with depth.If available, this evaluation can be tested with individual depth points such as drill core information.Additionally, an expansion towards 3-D would require probability attributes for dip azimuths.

Figure 1 .
Figure1.Geological map of the study area (modified afterBerger et al., 2017).The coordinate system of the inset map is given in WGS84, whereas the coordinates of the local geological map are given in the Swiss coordinate system (CH1903).

Figure 2 .
Figure 2. Schematic bloc diagram showing geometrical relationships between faults of different orientation groups (modified after Wehrens et al., 2017).

Figure 3 .
Figure3.Employed modelling workflow to generate a 3-D structural model of the area based on a surface lineament map.As a major step, the workflow also considers the uncertainty related to the connection between mapped faults at the surface and underground.

Figure 4 .
Figure 4. Schematic drawing of hypothetical example for validation of 3-D models based on angular and distance misfit (map view).The contours of the GTS are shown in grey with a mapped fault trace and a fault trace resulting from projection of the fault plane from the surface.

Figure 5 .
Figure 5. Schematic cross section illustrating the statistical modelling methodology for one example.(a) Reference state is defined by aforementioned workflow.A search window of 80 m is assigned to the reference underground point.(b) The code picks a possible underground point within the 80 m search window based on a normal Gaussian distribution and calculates the connection.(c) For every surface point one underground point is picked and connected by a projection line (hypothetical fault plane).By connecting each surface point with one underground point, a connection pattern is formed.Ten thousand runs were performed, each connecting different surface and underground points.The 10 000 connection patterns are compared with each other and evaluated by a penalty criterion.This penalty criterion addresses the number of crossings and rate solutions with a lowest and highest number of fault plane crossings, yielding a probability for connecting a specific surface point to a certain underground point.

Figure 6 .
Figure 6.(a) Lineament map of the study area with underground rock laboratory.Topography contours are based on swissALTI3D (reproduced by permission from swisstopo; BA17063).(b)-(d) Length-weighted rose diagrams showing the endpoint-to-endpoint strike of all lineaments (b), lineaments longer than 400 m (c), and lineaments shorter than 400 m (d).(e) Length-weighted rose diagram showing the orientation of each segment of all lineaments. Fig.11a

Figure 7 .
Figure 7. Surface fault map with faults grouped by strike orientation (group A, B, C).Hillshade image underlying the map is based on swissALTI3D (reproduced by permission from swisstopo; BA17063).Fault exposure lines are dashed over uncertain areas and labelled in cases for which a connection to GTS exists.Lower hemisphere equal area projection with plane poles grouped according to strike.The map is based on the Swiss coordinate system.

Figure 8 .
Figure 8.(a) Petrographic underground map.(b) Structural mapping (1 : 1000) of the underground rock laboratory (GTS) with faults grouped according to their strike.Indicated labels correspond to surface fault labelling and represent the maximum a posteriori interpolation.

Figure 9 .
Figure 9. Probability distributions of five selected examples: panels (a) to (d) show the highest probabilities achieved, whereas panel (e) shows an example without clear maximum probability.On top, the positions of the underground deformations zones are indicated and grouped according to their strike.Additionally, the maximum a posteriori interpolation is highlighted with an arrow.

Figure 10 .
Figure 10.Cross section showing maximum a posteriori connections between surface and underground faults.Faults are grouped and coloured according to their strike.Underground faults are represented by short ticks; the less transparent ones have a connection to the surface.

Figure 11
Figure 11.(a) Mountainside with incisions and exfoliation joints.(b) Detailed picture of underground outcrop showing outcrop conditions and key structural features: a ductile shear zone (SZ) and a fault gouge (FG).

Figure 12 .
Figure 12.Histogram showing the number of faults grouped per strike at the Earth's surface and underground (GTS).

AFigure 13 .
Figure 13.Comparison of maximum a posteriori interpolation with three extrapolation approaches used to assign dip to fault exposure line.Figure subdivided into (a) group A (NE-SW), (b) group B (E-W), and group C (NW-SE).Group B and group C are displayed jointly as group B, which contains only two faults.

Figure 14
Figure 14.(a) Box plot showing dip value for different extrapolation approaches and for maximum a posteriori (MAP) interpolation.(b) Box plots for dip comparison between different extrapolation approaches and the maximum a posteriori (MAP) interpolation.

Figure 15 .
Figure 15.Representation in 3-D of the maximum a posteriori model of fault geometry with three different angles of view.N is indicated by the black triangle.The black tunnel is 717 m long.