Multiscale lineament analysis and permeability heterogeneity of fractured crystalline basement blocks

. The multiscale analysis of lineament patterns helps defining the geometric scaling laws and the relationships between outcrop- and regional-scale structures in a fracture network. Here, we present a novel analytical and statistical 10 workflow to analyze the geometrical and spatial organization properties of the Rolvsnes granodiorite lineament (fracture) network in the crystalline basement of southwestern Norway (Bømlo Island). The network shows a scale-invariant spatial distribution described by a fractal dimension D ≈ 1.51, with lineament lengths distributed following a general scaling power - law (exponent α = 1.88). However, orientation-dependent analyses show that the identified sets vary their relative abundance and spatial organization/occupancy with scale, defining a hierarchical network. Lineament length, density, and intensity 15 distributions of each set follow power-law scaling laws characterized by their own exponents. Thus, our multiscale, orientation-dependent statistical approach can aid in the identification of the hierarchical structure of the fracture network, quantifying the spatial heterogeneity of lineament sets and their related regional- vs. local-scale relevance. These results, integrated with field petrophysical analyses of fracture lineaments, can effectively improve the detail and accuracy of permeability prediction of heterogeneously fractured media. Our results show also how the geological and geometrical 20 properties of the fracture network and analytical biases affect the results of multiscale analyses and how they must be critically assessed before extrapolating the conclusions to any other similar case study of fractured crystalline basement blocks.


Introduction
Crystalline rocks are characterized by very low intrinsic permeability, usually in the order of 10 -18 m 2 (Achtziger-Zupančič et 25 al., 2017;Brace, 1984), so that that their capability to transmit and/or store fluids is mainly related to the structural permeability associated with fracture and fault networks created by brittle deformation and the associated fluid-rock interaction (Caine et al., 1996;Caine and Tomusiak, 2003;Ceccato et al., 2021b, a;Evans et al., 1997;Pennacchioni et al., 2016;Schneeberger et al., 2018;Stober and Bucher, 2015). When studied at different scales, fracture and fault networks commonly exhibit variable geometrical and spatial characteristics, which may significantly affect the overall permeability 30 geometrical properties and the range of scales within which those laws can be applied; (iii) evaluate the difference between scaling laws retrieved from the entire network and those from individual sets. This multiscale and statistical approach tries to overcome the natural bias of lineament maps retrieved from remote sensing of natural outcrops, which are inherently incomplete due to partial exposure, resolution and analytical biases. The implications of the adoption of general scaling laws 80 on the upscaling/downscaling of fracture network properties, as well as the possible analytical biases and sources of errors in the analytical approach, are then evaluated and discussed. In this work, rather than trying to accurately quantify the scaling parameters, we focused on highlighting and analyzing the weaknesses and uncertainties that invariably accompany this sort of lineament analysis even when very robust statistical approaches and analytical methods are applied (Dichiarante et al., 2020). By integrating this information with existing field structural analyses and modeling of lineaments petrophysical 85 properties (Ceccato et al., 2021b, a;Scheiber and Viola, 2018), we provide further constraints on the multi-scale heterogeneity in magnitude, orientation, and spatial distribution of the permeability structure of the studied fractured crystalline basement.

Geological setting 90
The crystalline basement of the Island of Bømlo belongs to the Upper Allochthon units of the Caledonian orogen (Gee et al., 2008). Our lineament maps represent the fracture pattern affecting the Rolvsnes granodiorite, a pre-Scandian (466 ± 3 Ma; zircon U-Pb dating) granitoid pluton hosted in the Upper Allochthon metamorphic units   (Fig. 1). The 4 Rolvsnes granodiorite recorded a prolonged and multi-phase brittle deformation history Scheiber and Viola, 2018), only briefly summarized in the following, while the reader is referred to the cited literature for a more detailed 95 and comprehensive description of the tectonic history of the area. Overall, the whole tectonic history of the area is the expression of three main deformation episodes (Bell et al., 2014;Fossen et al., 2017Fossen et al., , 2021: (1) Caledonian convergence and continental collision from the Mid-Ordovician to the Silurian; (2) extensional tectonics related to the late-Scandian orogenic collapse during the Devonian, and (3) prolonged and multi-phase extensional tectonics related to the North Sea rifting from the Permian to the Cretaceous. During this tectonic evolution, the pre-Scandian Rolvsnes granodiorite did not record 100 penetrative ductile strain and was instead affected by pervasive brittle deformation. Each tectonic stage recorded in the granodiorite is associated with a characteristic set of fracture and fault zones that dissect Bømlo (Scheiber and Viola, 2018): (1) NNW-and WNW-striking conjugate strike-slip faults developed coevally with ENE-WSW and NE-SW-striking reverse faults during Caledonian convergence; (2) the same structures were reactivated with opposite kinematics during the early stages of late-Scandian orogenic collapse; (3) NW-and NNW-striking normal faults ascribable to the regional Permian-to-105 Jurassic rifting phase of the North Sea, which partially reactivated earlier, inherited structures. During the latest rifting stages of the North Sea, in the Early Cretaceous, new N-to NNE-striking fracture corridors and normal faults overprinted the previously formed fracture pattern. This tectonic history is reflected in the field by a sequence of three main classes of fracture and fault zones: (i) pre-Permian, ESE-WNW and NE-SW striking mineralized shear fractures and minor faults; (ii) Permo-Jurassic major normal faults, mainly NW-SE and N-S striking; (iii) Cretaceous fracture clusters striking N-S to NNE-110 SSW Scheiber and Viola, 2018).
A key structure for the detailed analysis of the timing of deformation, the geometry of the deformation structures, and the effects of deformation on the petrophysical properties of the crystalline basement of Bømlo, is the Goddo Fault Zone (GFZ, Fig. 1), (Ceccato et al., 2021b,a;Scheiber and Viola, 2018;Viola et al., 2016). The GFZ is an east-dipping normal fault that accommodated multiple slip increments during the prolonged Permian-to-Cretaceous rifting of the North Sea, recording 115 several stages of reactivation, during which a complex network of brittle structural facies developed in the fault core (sensu Tartaglia et al., 2020). Structures like the GFZ actually controlled the permeability and fluid flow evolution from rifting to current times of the crystalline basement (Ceccato et al., 2021b,a;Viola et al., 2016).
The Rolvsnes granodiorite is interpreted as the onshore analogue of the crystalline basement of the Utsira High, which, similarly to Bømlo, is affected by a complex fracture and fault network (Fredin et al., 2017;Trice et al., 2019). The Utsira 120 High crystalline basement is composed of (likely multiple) pre-Scandian igneous intrusions of similar age and composition to the Rolvsnes granodiorite (Lundmark et al., 2014;Slagstad et al., 2011). The fracture network in the Utsira High developed under tectonic conditions like those of the Bømlo crystalline basement, but with several significant differences mainly related to the structural position of the two crystalline basements within the North Sea rifting region (Bell et al., 2014). 125 5

Methods
In this study, the word "lineament" is used to refer to any linear feature of the topography as detected on a digital representation of the surface. The topography of Bømlo is cut through by deep linear grooves resulting from the penetration at depth of erosional processes exploiting fracture and fault zones of the crystalline basement (Scheiber and Viola, 2018). Therefore, the mapped lineaments represent fracture and fault zones identified in the field. In the followings, the terms 130 fracture(s) and fault zone(s) refer only to the geological structures observed in the field, which indeed assemble the remotely sensed lineaments. The lineament maps (Fig. 2a) used for the presented multiscale analyses have been generated in ArcGIS 10.8 by manually picking the same digital terrain model (DTM) of selected areas of the Island of Bømlo at different scales of observation. DTM's from high-resolution (1 m/pxl) airborne Light Detection and Ranging (LiDAR, kindly provided by Norges Geologiske Undersøkelse) surveys ( Fig. 1) have been used for the manual picking of lineaments at the 1:5,000; 135 1:25,000 and 1:100,000 scales. The details of LiDAR data acquisition and DTMs elaboration can be found in Scheiber et al. (2015). In addition, the dataset of lineaments interpreted from LiDAR DTM at the 1:5,000 scale was integrated with the interpretation of aerial orthophotos from the Bing Maps database (https://www.bing.com/maps). Bing aerial imagery was also adopted to distinguish between natural and man-made linear structures and to check for artefacts and potential misinterpretation of linear features on LiDAR-derived DTMs in the absence of systematic ground truthing. The 1:100 140 outcrop-scale lineament picking was performed on digital orthophotos of a key GFZ outcrop (Figs. 1, 3a) as obtained from the elaboration of the imagery collected via UAV-drone surveys through Structure-from-Motion (SfM) algorithms. Details on this acquisition and its elaboration methods can be found in Ceccato et al. (2021a). Topographic lineaments were traced as single, linear segments (not polylines) interpreting their topographic expression on DTMs. The obtained lineament maps are included in the dataset related to the present paper available at https://data.mendeley.com/datasets/3ymhkpmr9s/1. This 145 interpretation technique introduces two major analytical biases on the obtained lineament maps: 1) the interpreted length may only partially represent the entire lineament (which may be covered by deposits or be differently expressed in the topography, thus being not visible in its entire length, e.g., Cao and Lei, 2018; 2) as a consequence, abutting relationships, intersections between lineaments and lineament network topology and connectivity remain highly speculative and susceptible to subjective biases (Andrews et al., 2019). The orientation of mapped lineaments, expressed as azimuth angle 150 from the geographic north, was calculated in ArcGIS 10.8 using Easy Calculate 10 (https://www.ianko.com/free/free_arcgis.htm) and the Orientation Analysis Tools (https://is.muni.cz/www/lenka.koc/prvnistrana.html). Rose diagrams plotting lineament azimuths were produced with the MARD 1.0 software (Munro and Blenkinsop, 2012). Lineament density P 20 (m -2 ) and intensity P 21 (m/m 2 ) (Dershowitz and Herda, 1992) were calculated as the ratio between the total number of lineaments and total length of lineaments, respectively, over the total area of the land exposure in each 155 lineament map.

Fractal dimension -Box-counting method
The fractal dimension of each lineament map at different scales was computed with the box-counting method (Bonnet et al., 2001;Gillespie et al., 1993) by using the freely available function boxcount.m in MATLAB R2019b (http://www.fast.upsud.fr/~moisy/ml/boxcount/html/demo.html). The box-counting method consists in subdividing the analyzed image in 160 progressively smaller square boxes of side b and counting how many of them contain a segment of the analyzed lineament network. Plotting the number of boxes N b containing at least one lineament against the side length b on a log-log diagram should yield a straight curve, whose slope defines a power-law function with D as the fractal exponent (Bonnet et al., 2001).
The fractal dimension obtained from the box-counting method quantifies the scaling properties of the spatial occupancy of the lineament network (Bonnet et al., 2001). We have analyzed (i) the entire lineament maps at each scale of observation, 165 and, additionally, (ii) selected sub-areas of the lineament maps completely exposed on the land surface (that is, without any sea cover). Entire lineament maps inherently present discontinuous exposures of lineaments, which crop out on different islands of the Bømlo archipelago and are "separated" by sea branches and bays. Therefore, we have analyzed by boxcounting also the spatial occupancy (and the related fractal dimension D) of the exposed land surface of the Bømlo Island archipelago as shown in Fig. 1 in order to check if the fractal characteristics of the exposed land surface affect the fractal 170 parameters of the fracture network. The selected areas, instead, represent continuous exposures of lineaments over several m 2 or km 2 (depending on the scale of observation). For each lineament map we have then analyzed selected areas of continuous exposure (Fig. Supplement S1). For the lineament map of the GFZ outcrop we have analyzed eight sub-areas of 25 (four areas) and 100 (four areas) m 2 each. For the 1:5,000 lineament map, we have analyzed six sub-areas of 1 km 2 each. For the 1:25,000 lineament map, we have analyzed ten sub-areas of 1 (six areas) and 4 (four areas) km 2 each. For the 1:100,000 175 lineament map, we have analyzed four sub-areas of 4 (three areas) and 9 (1 area) km 2 each. Selected areas and lineament maps, and the complete results of the box-counting analyses are reported in the Fig. Supplement S1 and in the Online Repository (https://data.mendeley.com/datasets/3ymhkpmr9s/1).

Cumulative length distribution analyses
Length data of lineaments have been organized as cumulative distributions and plotted in log-log diagrams of the length L of 180 lineaments on the X-axis versus N (l>L) , the cumulative number of lineaments with length l > L (Fig. 2b-c). The cumulative length distributions were then normalized by the area of the land surface reported on each map over which the lineaments were picked. Single-scale distributions report the lineament lengths at a specific scale of observation. Multiscale distributions report the sum of the cumulative length distributions observed at different scales. We have analyzed the multiscale/single-scale normalized distribution functions of (i) all lineaments included in each lineament map at different 185 scales and (ii) each lineament orientation set at different scales.

Fitting of multiscale cumulative length distributions
The single-scale cumulative length distributions have been merged to form a single multiscale distribution. The mathematical functions fitting multiscale cumulative distributions are commonly retrieved by manual fitting of the distributions by assessing the slope of the tangent to the observed distributions (e.g., Bertrand et al., 2015;Bossennec et al., 190 2021;Castaing et al., 1996;Le Garzic et al., 2011). Here we apply a more quantitative method by adopting least square regression (LSR) in Microsoft Excel to the multiscale cumulative distributions. Manual fitting has been performed as well (results are reported and compared to LSR fitting in Fig. Supplement S2).

Fitting of single-scale cumulative length distributions
Single-scale cumulative length distributions have been analyzed by means of the Maximum Likelihood Estimation (MLE) 195 and Kolmogorov-Smirnov (KS) statistical tests to retrieve the best fitting mathematical function (Dichiarante et al., 2020;Kolyukhin and Torabi, 2013;Rizzo et al., 2017). The mathematical functions considered were negative exponential, powerlaw and log-normal (Fig. 2b). The advantage of adopting MLE-KS statistical tests derives from the possibility to also retrieve the function parameters (namely the exponent λ for the exponential, the exponent α for the power-law and the mean µ and standard deviation σ for the log-normal functions) in addition to the mathematical function best approximating the 200 observed cumulative length distributions. A dedicated MATLAB script implementing the freely available functions provided in the latest version of FracPaQ Rizzo et al., 2017) was used to this purpose. The results of the MLE-KS tests are reported in "checkerboard" diagrams, following the method proposed by Dichiarante et al. (2020) (Fig. 2d). Such diagrams allow to image the results of the MLE-KS tests on the selected portions of the cumulative distribution, i.e., the best fitting mathematical function for varying subdomains of the cumulative distribution. A subdomain is defined as a segment of 205 the cumulative distribution curve bounded by a lower and upper cut value ( Fig. 2c-d). The upper cut (UC) value represents the distance, expressed in terms of percentage of the total number of elements contained in the cumulative distribution, from the shortest observed length. The lower cut (LC) value represents the distance, in terms of percentage of the total number of elements contained in the cumulative distribution, from the longest observed length. On the checkerboard diagrams, the LC values are plotted versus the UC values (Fig. 2d). Each point of the checkboard represents a specific percentage range of the 210 total cumulative distribution between the upper and lower cut limits over which the MLE-KS tests were run. The plotted symbol represents the mathematical function among those considered (power-law, exponential, log-normal) for which the MLE-KS tests yielded the highest fitting score, whereas the symbol is color-coded according to the retrieved value of the fitting scores (namely the H-percentage HP and the P-percentage PP parameters, see Rizzo et al., 2017;Dichiarante et al., 2020). This analytical approach allows for the determination of the mathematical function that best fits the truncated 215 cumulative distribution and for the evaluation of the effects of truncation and censoring biases ( Fig. 2c-d). We report in the main text and tables the range of UC values for which each function fits best: the upper cut values quantify the "truncation" 8 of the cumulative distribution at short lengths, and this has been demonstrated to deeply affect the results of MLE-KS tests (Dichiarante et al., 2020).

Spatial distribution analysis 220
The spatial distribution of lineaments has been quantified following the approach by Sanderson and Peacock (2019). We analyzed the spacing between lineaments collected along virtual scanlines computed with the NetworkGT toolkit in QGis 3.12.2 (Nyberg et al., 2018) (e.g., Fig. 2e). The lineaments were classified and grouped into orientation sets. A grid of equally spaced virtual scanlines (100 m spacing for the 1:5,000 scale and 500 m spacing for the 1:25,000 and 1:100,000 scales) oriented perpendicular to the selected lineament set orientation was drawn upon the imported lineament map with 225 NetworkGT (e.g., Fig. 2e). Intersections between each virtual scanline and map lineaments were also obtained in NetworkGT. For each scanline, we analyzed the statistics (mean -µ, standard deviation -σ S and minimum and maximum values) for several parameters (Fig. 2e): (i) spacing (S) between lineaments; (ii) Coefficient of Variation (CoV) of the spacing, defined as the ratio between the standard deviation of spacing along a scanline and its average (CoV = σ S /µ) (Gillespie et al., 2001); (iii) coefficient of heterogeneity (V f ) and its statistical significance (V*) according to Sanderson and 230 Peacock (2019). The CoV of spacing is commonly adopted to assess the spatial organization (clustering vs. uniform distribution) of lineaments along scanlines (e.g. Gillespie et al., 2001). CoV values > 1 are usually related to clustered lineaments; CoV = 1 should represent a (negative) exponential-random distribution of spacing intervals, and CoV < 1 is usually related to log-normal (uniform) spacing distributions (Gillespie et al., 2001;McCaffrey et al., 2020;Odling et al., 1999). The spacing heterogeneity, i.e., the deviation of the spacing distribution along a scanline from a uniform distribution, 235 is quantified by the Vf and V* coefficients computed with the Kuiper method . The coefficient of heterogeneity V f quantifies the deviation from a theoretical uniform distribution of the observed spacing distribution along a given scanline expressed as the sum of the moduli of the positive and negative deviations (see Sanderson and Peacock, 2019;Fig. 2e). The coefficient V* quantifies the statistical significance of the heterogeneity factor V f : where N i represents the number of lineaments intersected by the scanline. Stephens (1970) demonstrated that for V* > 1.75, 2.0 and 2.3, the null hypothesis of uniformity can be rejected at the 95%, 99% and 99.9% levels, respectively. Thus, the coefficient V* can be used to quantify the probability that a certain spacing distribution is uniform or not.
We present the results of this analyses as CoV-V* diagrams (Fig. 2f) in which we plot the statistical distribution of the values of CoV vs. V* as box-and-whiskers. In doing so, we can qualitatively evaluate if, statistically, a set of lineaments has a 245 random or organized spatial distribution . With this method, four main spatial organization types can be distinguished ( Fig. 2f): (i) uniform distribution, characterized by CoV << 1 and V* < 1.75; (ii) random distribution, characterized by CoV ≈ 1, V* < 1.75; (iii) corridor/clustered distribution, characterized by CoV > 1 and V* > 9 1.75-2.00; (iv) fractal distribution, characterized by CoV >> 1 and V* >> 1.75-2.00. Scanlines with more than 10, 5 and 3 lineament intersections were considered on maps at 1:5,000, 1:25,000, and 1:100,000, respectively. 250 Given the limited number of intersections recorded by each virtual scanline in our maps (N i << 30), more advanced and upto-date analyses of the spacing variability (Marrett et al., 2018;Sanderson and Peacock, 2019;Bistacchi et al., 2020) were not possible. The Vf-V* method proposed by Sanderson and Peacock (2019) yields statistically meaningful results for datasets containing a number of intersected lineaments N i > 6. Our datasets unfortunately do not always satisfy this requirement. To perform the analysis on a statistically meaningful number of samples (total number of scanlines, N SL >10), 255 we had to reduce the minimum number of spacing data each scanline had to contain to be included in the dataset. Therefore, our analyses report the results from scanlines showing a minimum number of intersections N i >10, >5, and >3 for the analyses performed at the 1:5,000, 1:25,000 and 1:100,000 scale, respectively.
All parameters and the related abbreviations are reported in Table 1.

Lineament maps description
The manual picking of topographic lineaments on different digital representations of the selected areas of Bømlo led to the production of maps at different scales ( Fig. 3, Supplement S1-S3-S4). The orthophotos retrieved from UAV-drone surveys and the related lineament map (Fig. 3a, Supplement S3a) helped to characterize the main outcrop of the GFZ along the eastern shoreline of Goddo Island (Ceccato et al., 2021b,a;Viola et al., 2016). The investigated areas extend for 2127 m 2 265 over which we picked 930 lineaments. Lineament mapping on LiDAR DTM and aerial imagery at the 1:5,000 scale ( Fig. 3b, Supplement S3b) was performed on the best exposed areas along the coastline of the Goddo Island and nearby smaller islands. The resulting lineament map covers more than 17 km 2 and includes 3,835 lineaments. Furthermore, we generated additional lineament maps from the interpretation of the LiDAR DTM at the 1:25,000 and 1:100,000 scales over the same area (83 km 2 ; Fig. 3c-d, Supplement S3c-d). The 1:25,000 lineament map contains 894 lineaments, whereas the 1:100,000 270 map contains 249 lineaments.

Fractal dimension
The fractal dimension of the lineament maps at all scales was evaluated by applying the box-counting method (Bonnet et al., 2001;Gillespie et al., 1993). The number of filled boxes N b decreases with increasing box size b following a power-law relationship (Fig. 4). The power-law exponents (the fractal exponents) retrieved from the box counting analyses of the entire 275 lineament maps at different scales range between 1.45 and 1.61 (Fig. 4). On average, the lineament network exposed on the discontinuous outcrops of the Bømlo islands archipelago is characterized by a fractal dimension D = 1.51± 0.14 (2σ S ). The analysis of the spatial occupancy of the exposed land surface of the Bømlo islands archipelago yielded a fractal dimension D = 1.72 (Fig. 4, see supplementary dataset in the Online Repository). Figure 4 shows the box-and-whiskers diagram of the fractal exponent D for each scale of observation as obtained from the 280 analyses of selected areas. The fractal dimensions D obtained from the analyses of both entire lineament maps (cross marker in Fig. 4) and selected areas (box-and-whiskers in Fig. 4) are always smaller than the fractal dimension D of the exposed land surface on the Bømlo islands archipelago (D = 1.72; gray solid line in Fig 4). Similarly, the fractal dimension of the whole lineament maps at 1:5,000, 1:25,000 and 1:100,000 scales always over-estimate the fractal dimensions D obtained from the analyses of selected areas. This is not true for the GFZ lineament map at the 1:100 scale, in which the average 285 fractal dimension of the considered sub-areas is systematically larger than that of the entire lineament map (D = 1.46). The

Lineament Orientation
The comparison of the rose diagrams at different scales of observation allows to define some dominant trends ( Fig. 5a-b).
The five main orientation sets are (Fig. 5a, Supplement S4, Table 2 . 5a). Therefore, we have subdivided Set 1 into Set 1a, including NNE-SSW-striking lineaments, and Set 1b, including N-S-to NNW-SSE-striking lineaments. This subdivision will be adopted for discussing the spatial organization of the lineaments.

Cumulative length distributions
The results of LSR fitting and MLE-KS tests are reported and summarized in Fig. 6 and Table 3; the checkerboards diagrams 305 are reported in the supplementary material (Fig. Supplement S5).
The results of the MLE-KS tests suggest that a log-normal function best approximates the entire single-scale distribution in all considered cases (Fig. Supplement S5; Table 3). Variably truncated distributions are best approximated by either negative exponential or power-law functions (Table 3). In particular, the truncated length probability distributions for both single sets and the entire lineament network mapped at 1:100 are best represented by negative exponential functions, with λ ranging 310 between 0.65 and 1.25. Truncated distributions retrieved from lineament maps at 1:5,000 are best fitted, in most cases, by power-law functions with a minimum exponent α of 2.2. Truncated length distributions for lineaments mapped at 1:25,000 and 1:100,000 scale are well approximated by negative exponential functions, with an average λ of 0.004 and 0.0017, respectively (Table 3). Figure 6a reports the cumulative length distributions for the entire set of lineament maps normalized to the area of 315 investigation at each scale of observation. The multiscale normalized cumulative distributions obey a general power-law relationship valid over five orders of magnitude (1 m to 10,000 m). The power-law exponent α is 1.88 (Fig. 6a). Figure 6b reports the multiscale cumulative length distributions for each lineament set normalized for the area of investigation. Also in this case, multiscale distributions obey a general power law scaling with a characteristic exponent α for each set ranging between 1.62 and 2.12 (Fig. 6b, Supplement S2). 320

Lineament density and intensity
As also suggested by the cross-scale variation of the relative proportions of the orientation sets ( Fig. 5a-b), also the normalized density P 20 (m -2 ) and intensity P 21 (m/m 2 ) of each lineament set vary across scales. The variations of density and intensity are both described by a power-law relationship in log-log diagrams plotting the scale on the X-axis (e.g., 10 5 = 1:100,000) and the density P 20 or intensity P 21 on the Y-axis ( Fig. 7a-b) (e.g., Castaing et al., 1996). The variation trend for 325 the total lineament density P 20 of each map at different scales is characterized by power-law exponents β = 1.77 (Fig. 7a).
Sets 1, 2, and 3 display β values larger than the average value; Sets 4 and 5 display β values smaller than the average values.
Similarly, the variation trend for P 21 is characterized by a power-law exponent δ = 0.86 (Fig. 7b); Sets 1, 2, and 3 show δ values larger than the average value. Sets 4 and 5 display δ values smaller than the average value.

Spacing and organization at different scales 330
The CoV-V* diagrams highlight a similar trend for all analyzed lineament sets with increasing scale of observation (from 1:100,000 to 1:5,000). At the 1:5,000 scale (Fig. 8a), Set 1a and 1b lineaments are characterized by CoV ≤ 1 and V* ≤ 1.75, suggesting a random-to-uniform spatial distribution. At smaller scales, CoV for Set 1 exhibit a tendency towards random-touniform distribution (Fig. 8b-c). Set 2 lineaments display CoV on average >1 at the 1:5,000 scale, and V* >1.75 for a significant number of data (>40% of the total number of data), suggesting a clustered spatial distribution. At smaller scales, 335 both CoV and V* values generally decrease, although some of the analyzed scanlines still display CoV >1 and V* >1.75-2.00.
Set 3 lineaments are too scattered and sparce to allow for a meaningful analysis of their spatial arrangement and, therefore, they are not reported in Fig. 8. Set 4 lineaments mapped at the 1:5,000 scale on average show CoV values >1, but V* is rarely >1.75. At smaller scales, both CoV and V* decrease progressively. CoV and V* for Set 5 lineaments are generally <1 and <1.75, respectively. The most significative variation across scales in spatial distribution occurs for Set 2 and Set 4, both of 340 which exhibit a tendency towards clustering at the large scale (1:5,000; CoV > 1, V* > 1.75; Figs. 8a-9d), whereas they exhibit a tendency toward a random-to-uniform distribution at the smaller scales (1:25,000 and 1:100,000; CoV < 1; V* < 1.75; Figs. 8b-c, 9e). None of the lineament sets show a tendency to develop fractal behavior with a power-law spacing distribution.

Discussion 345
In the following, we firstly assess the geometrical characteristics, scaling laws and exponent values obtained for the Rolvsnes granodiorite fracture network. Then, we evaluate the possible biases affecting the analyzed datasets and their effect on the quantification of the fracture network geometrical properties. In addition, we discuss the implications of applying the scaling relationships to the quantification of fracturing and reservoir permeability at different scales by integrating our results with field geological data. 350

Characterization of geometric properties of the Rolvsnes granodiorite lineament network
The fractal dimensions D retrieved from the analysis of the 2D entire lineament maps cluster around 1.5 (Fig. 4), similar to what is commonly reported from other case studies on lineament pattern fractal dimensions (Bonnet et al., 2001;Hirata, 1989). Also, the normalized cumulative distribution of lineament lengths effectively defines a single scaling law, which can be best described by a power-law relationship with an exponent α = 1.88 ( Fig. 6a; Table 5). The general scaling law obtained 355 for the overall lineament network is very similar to that derived from many other case studies of fracture networks affecting both crystalline basements and (meta)sedimentary rocks, with an average power-law exponent close to α = 2 (cf. Bertrand et al., 2015;Bonnet et al., 2001;Bossennec et al., 2021;Chabani et al., 2021;Le Garzic et al., 2011;McCaffrey et al., 2020;Odling, 1997;Torabi and Berg, 2011). The power-law scaling relationship defined by the lineament density P 20 values is characterized by a power-law exponent β = 1.77, similar to the value of 1.7 commonly observed in many other fault 360 networks (Castaing et al., 1996;Bonnet et al., 2001, and references therein).

Scale-invariant lineament network
A similar fractal dimension D and power-law scaling relationship are commonly used as evidence for the occurrence of a fracture network whose geometrical properties (size of fractures, i.e., length, and spatial correlation and organization) are scale-invariant (Bonnet et al., 2001). This suggests that, at a first approximation, the documented lineament pattern in the 365 Rolvsnes granodiorite is self-similar at any scale of observation. However, the Rolvsnes granodiorite case appears to be more complex than it would seem at a first glance. The detailed analyses of lineament maps, fractal dimension D and the geometrical properties of orientation sets revealed similar patterns of variation through the scales of observation.

Lineament Types within a hierarchical fracture network
For each scale of observation, the box-counting analyses of sub-areas of the lineament maps yielded similar D values (2σ 370 always < 0.1), suggesting that the lineament network is rather homogeneous in space in terms of spatial occupancy.

13
Nonetheless, these results highlight a general decrease of the fractal dimension D from the largest (1:100) to the smallest (1:100,000) scale of observations (box-and-whiskers of Fig. 4). Even if there is a small variability in the fractal dimension D values, none of the measure fractal dimensions, either from small sub-areas or entire lineament maps, are similar to the fractal dimension D = 1.77 retrieved from the areal occupancy of the exposed land surface of the Bømlo Island archipelago. The relationship between fractal dimensions of exposed land surface and fracture maps remains to be understood and would deserve further analyses, which however go beyond the scope of the present paper. 390 The Rolvsnes granodiorite lineament network is composed of five main orientation sets with variable relative abundance, density, and intensity across scales (Figs. 5 and 7, Table 2). The observed variations of density and intensity are predictable and can be described by a general power-law function, the exponent of which is characteristic of each orientation set ( Fig. 7; Table 5). Even though the single-scale cumulative length distribution for each orientation set can be best approximated by scaling laws other than power-law (Table 3), the multiscale cumulative length distribution is best approximated by a power-395 law scaling relationship ( Fig. 6b; Table 5). Again, each orientation set is characterized by its own power-law exponent ( Fig.   6b; Table 5), which differs slightly from that computed for the entire lineament network. Some lineament sets display similar trends of variation of the relative abundance and intensity, such that they can be grouped into two main set types (Figs. 6b, 9; Table 5): (1) Type A includes Sets 1, 2 and 3, characterized by comparable P20 and P 21 variation trends across scales (β ≈ 1.90; δ ≈ 0.95) and length distributions (α ≥ 2); (2) Type B includes Sets 4 and 5, 400 characterized by comparable P 20 and P 21 variation trends across scales (β ≈ 1.60; δ ≈ 0.70) and length distributions (α < 1.7).
This classification into Type A and B lineament sets is not directly reflected in the CoV-V* diagrams (Fig. 8a-b), the latter rather suggesting a scale-dependent organization of spacing for each lineament set. Therefore, the observed density/intensity variation and spatial organization trends indicate the occurrence of a hierarchical (scale-dependent) organization of lineament sets within a network presenting overall scale-invariant geometrical properties 405 14 (e.g., Le Garzic et al., 2011). In this hierarchy, Type B lineaments represent the higher-order structures, controlling the geometrical properties of the network at the regional scale (Fig. 9d-e). Type A lineaments represent lower-order structures and control the geometrical properties of the network at the local-to-outcrop scale (Fig. 9d-e).
At the smallest investigated scale, the homogeneously spaced, WNW-to-NW-striking Type B lineaments (Fig. 9d-e) dominate the network. These lineaments are characteristic of -and predominant over -the whole of onshore western and 410 southwestern Norway (Gabrielsen et al., 2002;Gabrielsen and Braathen, 2014;Tartaglia et al., 2022), as well as offshore (Preiss and Adam, 2021). The power-law exponent for Type B lineaments (α < 1.7) suggests that long lineaments represent a substantial part of the overall lineament population. These observations also suggest that Type B lineaments probably result from the homogeneous distribution of deformation at the regional-scale, while still representing localized zones accommodating significant deformation at the outcrop-scale, when compared to Type A structures (Ackermann et al., 2001). 415 Therefore, these lineaments probably represent major fractures and normal fault zones formed and repeatedly reactivated during the prolonged brittle tectonic history of the Rolvsnes granodiorite (Ceccato et al., 2021a;Preiss and Adam, 2021;Scheiber and Viola, 2018;Viola et al., 2016). The schematic representation of lineaments in Fig. 9d-e highlights an heterogeneous distribution of Type B lineaments, which is not captured by the statistical analysis of spacing heterogeneity.
Indeed, the Rolvsnes granodiorite can be subdivided into several domains of the lineament maps where either Set 4 or Set 5 420 lineaments are predominant at the regional scale ("Set 4-5 domain" -grey and dashed transparent areas in Fig. 9d-e). These domain-type distribution of regional lineaments were already reported by Scheiber and Viola (2018). At the largest analyzed scale, the lineament network is mainly dominated by random-to-clustered, NNW-SSE to NE-SW-striking Type A lineament sets (Fig. 9d-e). Accordingly, the general power-law exponent (α ≥ 2) suggests that, among Type A lineaments, short lineaments represent a significant part of that population, probably resulting from an incipient stage of distributed faulting 425 and deformation accommodation (e.g., Ackermann et al., 2001).

Analysis of reliability, biases and limitations of the scaling laws
The scaling laws described here have several limitations in their applicability related to: (i) occurrence of different orientation sets; (ii) network heterogeneity at different scales; (iii) analytical biases. 430

Reliability and biases behind orientation set definitions
The lineament sets defined by this study are grouped in classes based on their azimuthal orientation. As such, they may or may not share a genetic relationship. However, field analyses (Scheiber et al., 2015Scheiber and Viola, 2018) have demonstrated the occurrence of systematic sets of fractures, genetically related in terms of chronology, tectonic phase, and orientation, which can be identified from remote sensing techniques. Thus, we assume that the identified orientation sets 435 effectively represent groups of genetically correlated fractures (see below, Section 5.3).

15
The lineament orientation reported by the 1:100 rose diagram differs significantly from the orientations of all other diagrams (Fig. 5). Even though the number of lineaments interpreted from UAV-drone imagery is statistically significant (N Lin = 930), the N-S-trending outcrop exposure, its 3D topography, and the location of the exposed area along a major fault zone (Ceccato et al., 2021b,a) are such that it is necessary to question whether the obtained results are truly representative of the 440 larger-scale lineament network. The observed variations of density, intensity, and relative abundance of orientation sets across scales could be affected by several analytical and interpretative orientation-dependent and classification biases.
First of all, under-sampling of specific lineament orientations during manual interpretation may be due to (Scheiber et al., 2015): (i) interpretative biases of the operator; (ii) changes in resolution of the digital representation of the terrain (DTMs and orthophotos) with the changing scale of observation; (iii) constant direction of the light source adopted for the LiDAR 445 DTM hill shading (from NW in our study). The change in resolution would affect equally each orientation set, thus maintaining a constant relative abundance across scales. Likely, constant direction of the light source may affect the detection of lineaments at specific orientations, but systematic effects have not been identified by previous studies (Scheiber et al., 2015). Also, the topography, surface extension, perimetral shape and exposure of outcrops with rough topography (such as the GFZ outcrop analyzed here at 1:100 scale) might affect the exposure and detection of specific lineaments and 450 thus orientation sets -some sets may be more visible than others. In particular, gently dipping fractures or fractures parallel to the surface of the outcrop might be underrepresented. Moreover, rose diagrams only report the number of lineaments, and do not consider their spatial persistence (length), such that fractures related to the main GFZ and expected to be dominant at the local scale might be represented by very few, but longer lineaments. Thus, the small number of lineaments can be diluted and obscured by the large number of short lineaments related to background fracturing. 455 This notwithstanding, by considering the relative abundance of orientation sets retrieved only from the 1:5,000; 1:25,000 and 1:100,000 lineament maps, we can still observe the relative decrease of Sets 1, 2, 3 along with the increase of Set 5 lineaments at progressively smaller scales of observation (Fig. 4b). This suggests that the observed variation trends reflect an effective variation in relative abundance of lineaments across scales and represent a real characteristic of the lineament network. 460

Network vs. orientation set scaling laws
The studied lineament network exhibits some general power-law relationships describing the multiscale behavior of both lineament length distribution (α = 1.88), density P20 (β = 1.77) and intensity P 21 (δ = 0.86) ( Table 5). These general powerlaw scaling laws may effectively be adopted to retrieve fracture network properties (geometrical properties and permeability) at any scale of observation. However, the adoption of a general scaling law for the geometrical properties without taking into 465 consideration the peculiarity of each orientation set building up the network, may lead to an erroneous extrapolation of the analyzed properties. For example, lineament sets exhibit different power-law exponents for density P 20 and intensity P 21 distributions, which, in our case study are systematically smaller for Type A sets, and larger for Type B sets than the exponent of the network taken as a whole (Fig. 7). Adopting power-law exponents larger than the actual exponent of the individual lineament set would lead to an overestimation of the network properties at larger scales, and vice-versa. In the 470 case of the Rolvsnes granodiorite lineament network, this overestimation/underestimation could be significant and reach one order of magnitude in terms of intensity and density (Fig. 7).

Spatial organization and scaling limitations
Field investigations (Ceccato et al., 2021b,a;Scheiber and Viola, 2018) have revealed the highly heterogeneous distribution of fractures at the outcrop scale. Most of the identified fracture sets at the outcrop occur with either a clustered spatial 475 organization or a variable intensity over short distances (50-100 m; Ceccato et al., 2021a). This clearly represents a limitation to the acritical extrapolation of the general power-law determined in this study, and thus the lower bound for the application of the proposed power-law scaling (Bonnet et al., 2001). Similarly, the "zonal" spatial distribution of Set 4 and 5 lineaments identified at small scale of observation (see the identified "Set 4 domain" and "Set 5 domain" reported in Fig. 9de) needs to be accounted for when evaluating the upper limit of applicability of the general scaling laws defined here. The 480 outcrop-scale spatial heterogeneity and the overestimation/underestimation effects due to applying a general power-law scaling become relevant when considering the role that different fracture sets may have in the definition of the net permeability of a fractured crystalline basement, as highlighted by field studies (Ceccato et al., 2021b,a;Gabrielsen and Braathen, 2014;Torabi et al., 2018).

Analytical biases and statistical methods to overcome them 485
In the light of the discussion above, it is therefore important to analyze the mapped distributions with appropriate methods and to always consider the regional lineament network as composed of different orientation sets, each of which is characterized by its own geometrical and scaling properties. Additionally, we need to consider all potential sources of error and estimate the analytical/human biases that potentially affect our analyses.
It is worth noting that, manual fitting of our multiscale distribution slightly overestimates the power-law exponents retrieved 490 from LSR methods (Fig. Supplement S2). This might be since LSR methods consider the entire distribution, including the portions of the distributions affected by censoring and truncation, which are inherently excluded by manual fitting.
Most published length datasets and the related fitting results have been analyzed qualitatively (e.g., manual multiscale fitting or LSR of truncated distributions). Indeed, the review process of similar published datasets (for example, particle size distributions) with updated statistical methods has revealed a fundamental mis-fitting of mathematical distributions, and 495 therefore a fundamental misunderstanding of the mathematical function describing the distribution of the investigated geological parameter (Phillips and Williams, 2021). MLE-KS tests have already demonstrated their strength in the analysis of fault attribute distributions (Dichiarante et al., 2020;Kolyukhin and Torabi, 2013). In our case, the results of length cumulative distributions fitting with MLE-KS tests differ significantly from the multiscale LSR power-law relations reported in Fig. 6. Even when the distribution is best approximated by a power-law function (e.g., most lineament sets mapped at 500 1:5,000; Table 3), the values of the power-law exponents retrieved from MLE-KS tests (α >2.2) differ from those obtained from the LSR fitting of multiscale distributions (α < 2.12).
Even though we might apply very robust statistical methods for the evaluation of the significance of the results of lineament map analyses, we still struggle to accurately quantify the scaling parameters or to define general, statistically robust scaling laws. Nevertheless, the first step to improve our accuracy in the quantification of scaling laws and scaling parameters is to 505 consider and discuss the most important biases and analytical errors that might affect the analysis results and lead to deviation from the "natural" scaling laws.
The observed deviations in curve fitting results (both that of single-scale distributions from the power-law, and that of the single-scale power-law exponents from the multiscale ones) are commonly observed in almost all natural fracture networks.
Power-law fitting is usually retrieved from short segments in the central portions of a "truncated" cumulative distribution 515 (e.g. Dichiarante et al., 2020). Truncation and censoring biases may affect large portions (even >50% of data) of the cumulative distribution. This would mean that most (>>50% of data) of the analyzed dataset is biased, and thus of little use to any kind of statistically significant analysis, such that the related fitting results are not statistically meaningful. The multiscale distribution analysis can reveal the overarching scaling law of the fracture network and of each lineament set, but the actual values of the exponents of the fitting laws need to be carefully evaluated by also considering the statistical 520 significance of the analyzed dataset.
Segmentation of long lineaments into shorter segments may be due to several causes, both introduced into the dataset by analytical/interpretative biases, and intrinsically related to the network topology, fracture chronology and geological fracture formation processes. Segmentation may result from partial exposure and cover of the fracture network, and it may increase the power-law scaling exponent, without affecting the type of scaling-law function (Cao and Lei, 2018). Segmentation may 525 be related to the progressive growth stages of fault/joint patterns evolving with increasing accommodated deformation and faulting maturity from a network composed of completely isolated short fractures to a network formed by a few long, single lineaments, through fracture interaction and interconnection (Ackermann et al., 2001;Michas et al., 2015;Scholz, 2002).
This has been demonstrated to affect both the shape of the mathematical function describing the length distribution (exponential vs. power-law), as well as the power-law exponent at a specific scale of observation (Schultz et al., 2013). 530 However, this may explain the difference in scaling relationships observed during the evolution of a fracture network through time but not at different scales of observation. In addition, the subjective choice of tracing single segments composing a longer lineament as separate fractures rather than tracing a single, continuous, long lineament, may likely affect the cumulative length distributions of the fracture network. Tracing single segments would increase the number of short segments compared to longer segments, at constant P 21 intensity, increasing the total number of traced lineaments, and thus 535 increasing the power-law exponent of the distribution (Xu et al., 2006). This segmentation bias may justify the fact that power-law exponents of the LSR multiscale length distributions of each fracture set (Fig. 5b) are systematically smaller than those obtained from single-scale MLE-KS tests at 1:5,000 scale. Whether or not this sampling bias may effectively affect the mathematical shape of the cumulative distribution would deserve further investigations, which goes beyond the scope of the present paper. 540 Thus, the most plausible option is that the fracture network may effectively present different scaling properties at different scales of observations (Kruhl, 2013). Indeed, fault and fracture networks may exhibit a hierarchical organization, which inherently implies scale-dependent geometrical properties and spatial distribution of lineaments (Castaing et al., 1996;Le Garzic et al., 2011). In fact, this is also consistent with the observed variation of relative abundances of orientation sets across scales: each lineament set contributes differently to the overall fracture network geometrical characters and thus the 545 variation of the relative abundance may also lead to variations in geometrical properties (spatial organization and length distributions) at different scales (e.g., Le Garzic et al., 2011).

Integration of remote sensing and field observations
The fracture pattern of the Rolvsnes granodiorite includes three main classes of fractures and fault zones Scheiber and Viola, 2018): (i) pre-Permian, ESE-WNW and NE-SW striking mineralized shear fractures and minor 550 faults; (ii) Permo-Jurassic major normal faults, mainly striking NW-SE and N-S; (iii) Cretaceous fracture clusters striking N-S to NNE-SSW.
Smaller normal faults and mineralized shear fractures described by Scheiber and Viola (2018) are subparallel to the lineament Sets 2-3-4 defined here. They form the background fracture pattern of the Rolvsnes granodiorite (Ceccato et al., 2021a). In particular, Set 4 lineaments are subparallel to the ESE-WNW orientation of the (relatively) oldest generation of 555 fractures identified in the field (Scheiber and Viola, 2018). Indeed, the regional distribution of Set 4 and the geometrical characteristics discussed above stress their importance as regional structures accommodating significant deformation through the brittle deformation history of the Rolvsnes granodiorite (Scheiber and Viola, 2018). Conversely, NE-SW minor faults and shear fractures, subparallel to Sets 2 and 3 lineaments, accommodated only limited deformation, which may be compatible with their local-scale distribution and geometrical characteristics (e.g., Ackermann et al., 2001). Permo-Jurassic 560 normal fault zones are oriented NW-SE to N-S, similarly to our Set 5 and Set 1(b) lineaments. The geometrical characters and spatial distribution of Set 5 lineaments suggest their role as important zones of deformation accommodation at the regional scale. N-S to NNE-SSW fracture clusters are comparable to Set 1(a) orientations; they are indeed local-scale structural features and are inferred to have accommodated limited deformation during Cretaceous rifting of the North Sea (Scheiber and Viola, 2018). 565 Summarizing, although our lineament classification in Type A and B is certainly an oversimplification of the actual complexity of the natural fracture network, it provides valuable information as to the geometrical characteristics of faults and fractures and their regional/local importance as zones of deformation accommodation.

Constraints on the multiscale permeability structure of crystalline basement
As constrained by field studies, each lineament/fracture set may contribute differently to the bulk permeability of a fractured 570 crystalline basement block (Ceccato et al., 2021b,a;Gabrielsen and Braathen, 2014;Torabi et al., 2018). The effects of the background fracturing of the Rolvsnes granodiorite on permeability is secondary (Ceccato et al., 2021a), it being mainly composed of minor sealed faults and mineralized fractures belonging to Set 2-3-4 lineaments. It is the regional scale structures, like our Type B lineaments (Set 5; e.g., the GFZ) that effectively control permeability, fluid flow, and reservoir compartmentalization at the regional scale (Ceccato et al., 2021a;Holdsworth et al., 2019). Results from in-situ 575 petrophysical analyses and discrete fracture network modelling of fault zone permeability have shown that these structures behave as mixed conduit-barrier for fluid flow, and they are characterized by a strongly anisotropic permeability tensor (Caine et al., 1996;Ceccato et al., 2021a). Fluid flow is promoted parallel to the main fault plane, especially parallel to the (sub-horizontal) intersection directions of the dominant fracture sets within the fault damage zone, whereas the anisotropic permeability of the fault core brittle structural facies buffers cross-fault fluid flow (Ceccato et al., 2021b,a;Tartaglia et al., 580 2020). Conversely, fracture clusters, comparable to our Set 1(a) lineaments, may represent effective fluid pathways at the outcrop scale, acting as preferential conduits for vertical fluid flow within the basement (Ceccato et al., 2021b, a;Torabi et al., 2018;Place et al., 2016;Souque et al., 2019).
In summary, by integrating field and remote sensing data we can improve the conceptual models and their dimensioning in an attempt to describe the anisotropic permeability structure of a fractured crystalline basement at different scales (e.g., Fig.  585 11 of Ceccato et al., 2021b). Our results constrain the heterogeneous structure of a fractured basement block in terms of orientation and spatial distribution of permeability. The permeability of the fractured basement at the regional-scale is characterized by the occurrence of rhombohedral-shaped compartments (the fault-bounded polygonal domains of Ceccato et al., 2021b) that are homogeneously distributed and defined by the higher-hierarchical order Type B lineaments. Their extension is determined by the spacing of Type B lineaments, ranging in the order of 500-1,000 m at the regional scale 590 (Table 4). Fluid flow is promoted along the major fault planes and parallel to the sub-horizontal intersections of fracture sets dominant within the fault damage zones (Ceccato et al., 2021a). Within these rhombohedral compartments, permeability is heterogeneously distributed at the 50-100 m scale, following the random-to-uniform spacing distribution of lowerhierarchical order Type A (Set 1b) lineaments (Table 4). At the outcrop scale, these N-S lineaments are represented by fracture clusters, which promote vertical fluid flow (Ceccato et al., 2021b,a). 595 Accordingly, any underestimation/overestimation of spatial distribution and density of the lineaments may deeply affect the accuracy of hydrological and petrophysical models of fractured basement blocks at the outcrop and at the sub-seismic resolution scale (Bertrand et al., 2015;Le Garzic et al., 2011).

Conclusions 600
The fractured crystalline basement of the Rolvsnes granodiorite on Bømlo is characterized by the occurrence of a fractal fracture network controlled by a general power-law scaling law for the distribution of fracture lengths. However, detailed orientation-dependent analyses have revealed that this first-approximation scale-invariant lineament network is composed of lineament sets, which individually exhibit a scale-dependent hierarchical spatial distribution, and parameter variation trends with the scale of observation. Different trends of intensity/density variation across scales for each orientation set have been 605 detected, as well as different scaling laws for length distribution of each orientation set. These results, integrated with field observations, suggest that the documented lineament network results from the summation of different geological structures (e.g., faults vs. joints, major fault zones vs. incipient minor faults) organized in a hierarchical manner and characterized by different geometrical and scale-dependent properties.
The hierarchical lineament network affecting the Rolvsnes granodiorite controls the anisotropy and directionality of the 610 permeability structure of the basement at different scales. At the regional scale, the crystalline basement is characterized by a rhombohedral pattern of basement compartments bounded by regional fault zones impermeable to cross-fault fluid flow.
Within these compartments, the permeability structure is controlled by local-scale fracture clusters, promoting subvertical N-S striking fluid flow.
Our study allows us to draw some general conclusions about the methods for characterization of fracture network and their 615 analysis: -Firstly, the presented multiscale analytical workflow may represent a valid option for the quantification of large, inherently incomplete (due to analytical and subjective biases) lineament datasets. The lineament maps retrieved from digital terrain and surface models of the Rolvsnes granodiorite offer very large datasets, which are inherently incomplete due to partial exposure and/or incomplete sampling of lineament due to partial exposure, resolution or 620 human biases. Thus, a statistical approach such as that proposed in this paper is highly recommended when aiming to retrieve relevant information from datasets that, for several reasons, are only partially representative for the entire fracture network.
-Detailed orientation-dependent, multiscale analyses of the lineament network can provide the different scaling laws and geometrical properties for each constituent lineament set, which can be adopted to improve the detail and tune 625 the accuracy of permeability models of fractured crystalline basements considering outcrop-scale structural features.

21
-The integration of multiscale length distribution analyses, multiscale intensity/density estimations and multiscale description of spatial organization provides useful information for the classification of topographic lineaments as different geological structures (e.g., fracture clusters vs fault zones) with specific hierarchy and control on the permeability of the fractured basement. 630  Table 1. Summary table of the parameters and related nomenclature adopted in this paper. Table 2. Table presenting Supplementary Fig. S1, the whole dataset is available at the Online Repository).