Establishing an integrated workflow identifying and linking surface and subsurface lineaments for mineral exploration under cover: Example from the Gawler Craton, South Australia

. Mineral exploration in areas comprising thick and complex cover represents an intrinsic challenge. Cost- and time-efficient methods that help to narrow down exploration areas are therefore of particular interest to the Australian mining industry and for mineral exploration worldwide. Based on a case study around the Tarcoola gold mine in the regolith dominated South Australian central Gawler Craton, we suggest an exploration targeting workflow based on the joint analysis of surface and subsurface lineaments. The datasets utilised in this study are a digital elevation model and radiometric data that 5 represent surface signals and total magnetic intensity and gravity attributed to subsurface signals. We compare automatically and manually mapped lineament sets derived from remotely sensed data. In order to establish an integrated concept for exploration through cover based on the best-suited lineament data, we will point out the most striking differences between the automatically and manually detected lineaments and compare the datasets that represent surficial in contrast to subsurface structures. We further show how lineaments derived from surface and subsurface datasets can be combined to obtain targeting 10 maps that help to narrow down areas for mineral exploration. We propose that target areas are represented by high lineament densities which are adjacent to regions comprising high density of lineament intersections.


Introduction
The Gawler Craton (Figure 1a) is one of the three largest Archean to Proterozoic Cratons within the Australian Continent ( Figure 1b) and the major crustal province in South Australia (Hand et al., 2007). Australian Cratons are rich in mineralisation 15 (e.g., Pilbara and Yilgarn Craton (Witt et al., 1998)), and the Gawler Craton hosts significant economic mineralisation such as Olympic Dam, Challenger, Prominent Hill, and Tarcoola. The discovery of new deposits is particularly challenging in this part of Australia due to limited surface outcrops, variable thickness and complexity of the cover, with few surface features that can be used as direct proxies for mineral exploration.
To narrow down potential areas for exploration, and to enhance the general understanding of the geology, the Geological Sur-20 vey of South Australia (GSSA) has recently acquired high-resolution magnetic, radiometric, and digital elevation data across 1 the Gawler Craton via the Gawler Craton Airborne Survey (GCAS) . Here we utilise these and previously acquired gravity datasets for the extraction of lineaments from the surface (digital elevation model (DEM) and radiometrics) and subsurface (magnetics and gravity) datasets. Lineaments are linear features that can reflect geological structures and the extraction of such features can be important for mineral exploration, as well as for the investigation of fault activity (neotectonic) 25 and water resource analysis (Vassilas et al., 2002). In images, photos, or maps, lineaments are represented by straight or slightly curved lines, linear patterns or an alignment of discontinuity patterns (Wang, 1993). A relationship between lineaments and mineralisation has been suggested for a long time and has been proven to be a useful tool for mineral exploration (O'Driscoll, 1986). We can distinguish between surface lineaments that are obtained from surface data, and geophysical lineaments that are derived from processed geophysical data. It is widely assumed that surface lineaments represent structural features which in 30 the simplest model are related to dip-slip or strike-slip faults (Florinsky, 2016). This assumes that a considerable displacement is associated with faulting in the subsurface that leads to a detectable pattern on the surface (Boucher, 1997). In contrast, geophysical lineaments represent major subsurface boundaries (e.g. Hall, 1986;Langenheim and Hildenbrand, 1997) that are not necessarily associated with faults, but rather with lithological or petrophysical contrasts. It is important to note that datasets such as digital elevation models and radiometrics represent only topographical variations and changes in chemical composition 35 of the surface geology or cover respectively.
As a case study for mineral exploration, we choose an area in the South Australian central Gawler craton around the Tarcoola mine, an Au-deposit mined for over 125 years (Daly et al., 1990). As the structurally controlled mineralisation often localises around discontinuity intersections (Wilson et al., 2018), this area represents a perfect study area for investigating the potential of surface and subsurface lineaments as a potential exploration targeting tool. The thick and complex cover (up to 500 m) overlying 40 the basement units makes it particularly challenging to identify target areas and a cost-efficient approach to exploration targeting is desirable in such a region. A great challenge is that surface impressions of non-vertical basement hosted displacement structures can be offset compared to the location of the discontinuity in the basement. Furthermore, cross-strike features are likely associated with small scale shear zones that will not be traceable in potential field data that images basement structures.
A reliable interpretation of surface and subsurface lineament sets is particularity challenging in an old crustal block such as the 45 Gawler Craton and we consider the work presented here as an first attempt to unify a lineament-based workflow for exploration targeting in such an environment.
In this study, we use the above-mentioned new GCAS datasets to identify surface and subsurface lineament features and design a workflow to automatically extract and analyse these features. We assume that elevation and radiometric data relates to surficial features, while gravity and magnetics data represent structures below the cover. This study is part of a broader 50 effort to geologically link basement architecture with surface linear features, landforms, and landscape variability in the central Gawler Craton (González-Álvarez et al., 2020(González-Álvarez et al., , 2022a. It is important to note that datasets such as digital elevation models, or radiometric data represent only the change in surface properties such as elevation or surface geology. These lineaments may or may not represent structures that extend into the subsurface. By using datasets that represent the subsurface (i.e. gravity and magnetics) lineaments extracted are directly representative of changes in the subsurface. The challenge is (1) identifying if the 55 lineaments from any dataset are geologically meaningful and (2) if lineaments from surface and subsurface datasets represent the same structure (e.g. fault, lithological boundary).
Towards that end, we further explore the generation and use of targeting maps (i.e., a map generated to highlight areas with specific features) based on surface and subsurface lineaments. Targeting maps derived from lineament analysis are often based on the density of lineaments per unit area. Density maps combining subsurface lineaments (potential field data) and 60 surface lineaments (digital elevation model and satellite imagery) were proposed as an exploration tool for groundwater (Epuh et al., 2020) and for mineral exploration (Mohammadpour et al., 2020). Lineament intersections were also used previously for the analysis of groundwater (Ilugbo and Adebiyi, 2017) and locations of intersecting structural elements were suggested to represent favourable target areas for mineral exploration (Krapf and González-Álvarez, 2018;González-Álvarez et al., 2019;Sheikhrahimi et al., 2019;González-Álvarez et al., 2020González-Álvarez et al., , 2022b. In hydrocarbon exploration, cross-strike discontinuities were 65 suggested as an exploration tool for natural gas (Wheeler, 1980). In the context of this contribution, we define targeting maps as maps which highlight areas that comprise the structural features preferable for mineral exploration.
We apply edge enhancement filtering to digital elevation data and perform manual lineament extraction on the DEM and automatic segmentation on all datasets to compare both approaches and results. We discuss the advantages and shortcomings of different edge detection filters applied to the subsurface and surface datasets and present a workflow ( Figure 1c) to help with 70 the identification of linear features that could be linked to main basement geological features.

Geological Overview
The Gawler Craton is the oldest and largest geological province in South Australia and represents one of the three major Australian Cratons (Figure 1b). The region hosts several economic iron oxide copper-gold ore deposits (IOCG) including the world-class Olympic Dam (Figure 1a). The Olympic IOCG province forms a 100 to 200 km wide north-south trending belt at 75 the eastern margin of the Gawler Craton . The largest known mineral occurrence in the area investigated is the Tarcoola mine, which hosts disseminated or veinlet-type Au-mineralisation mainly in brittle to brittle-ductile faults and shears (Hand et al., 2007). The age of the Tarcoola Au-mineralisation is considered to be 1564 Ma (Bockmann et al., 2019) and the timing of the ore-formation at the Tunkilla project (approx. 70 km S/SE of Tarcoola) has been constrained at 1590-1570 Ma (Budd and Fraser, 2004). Ongoing exploration in the central Gawler Craton targets Au, Cu-Au, Pb-Zn, 80 Fe, and Ni in the crystalline basement (Sheard et al., 2008). The expected commodities within the investigated region are mainly Au and Fe and are likely located in proximity to crustal-scale structures that provided conduits for the upwelling of deep crustal fluids. Such large scale reactivated tectonic features often form major crustal boundaries that are detectable with potential-field methods (Motta et al., 2019). It was shown that Archean gold mineralisations are often associated with such crustal-scale shear zones (Eisenlohr et al., 1989;Budd and Fraser, 2004). If a surface expression of such structures exists, 85 this may be indicative that the crustal structures remained active for a long time resulting in a high amount of deformation or strong lithological contrasts. In the framework of the central Gawler mineral-systems, the vicinity of such structures indicates potential exploration targets.
Three major orogenic events, corresponding to crustal deformation and tectono-thermal alterations, are recorded by the crystalline basement of the Gawler Craton: the Sleaford Orogeny (Paleoproterozoic, 2440 Ma), the Kimban Orogeny (Paleo-90 proterozoic, 1845-1700 Ma), and the Kararan Orogeny (Mesoproterozoic, 1650-1540 Ma) (Ferris et al., 2002;Swain et al., 2005;Kositcin, 2010;Reid et al., 2014). Three major mineral systems are related to these three deformation-magmatic events and the Tarcoola mineral field is attributed only to the Kararan system (≈ 1570 Ma) (Gum, 2019). As outlined by Hand et al. (2007), the exact timing and spatial distribution of the tectonostratigraphic sequences within the Gawler Craton remain a controversy, and reworking of existing mineralisation during reactivation of existing crustal-scale fluid conduits must be 95 considered (Gum, 2019). The last large-scale deformation in the Gawler Craton was the reactivation of shear zones between 1470 and 1450 Ma (Hand et al., 2007). After this time, only minor near-surface movements are recorded (Sheard et al., 2008).
Given the mineralisation is linked to the youngest orogeny in the area (Kararan Orogeny) and only minor tectonic activity is evident in the area after this event, we can assume that links between surficial and subsurface features point to areas of high deformation and/or neotectonic activity.

100
The bulk of the Paleo-to Mesoproterozoic rocks of the Gawler Craton enclose an Archean core in the central Gawler Craton with the oldest units being of Late Archean age (Reid et al., 2014). Internally, the Gawler Craton is subdivided into different domains based on contrasts in geophysical, structural, and geochemical characteristics (Ferris et al., 2002;Fairclough et al., 2003;Kositcin, 2010). The different rock units are often separated by crustal-scale shear zones that often coincide with the boundaries of the individual blocks. The region around Tarcoola includes four major provinces of the central Gawler Craton:

Basement and cover sequence
The key geologic features that we seek to explore are structural and lithology discontinuities in the basement along with if and how they may relate to today's landscape and current topographic relief. In the following, we briefly describe the dominant lithologies in the study area from the oldest to the youngest unit and highlight the expected variability in aeromagnetic and 110 gravity data.
The interpreted Precambrian basement geology of the Harris Greenstone Belt comprises east-northeast-trending linear magnetic highs that often correlate with broad gravity signatures and are flanked by ovoid to elongated magnetic lows and highs (Hoatson et al., 2002). Rocks of the Harris Greenstone belt are found along the Flinke and Yerda Shear Zones and at the eastern margin in the central part of Figure 2b.  (Huang et al., 2015, and references therein). The youngest structure in the region is the Mulgathing Trough in the northwest that has an inferred Permian age.

120
The magnetic signatures could help to distinguish between different rock types as suggested by Hoatson et al. (2002), who showed that Archean granitic gneisses and granites are often irregular or elongated bodies with low amplitude magnetic Proterozoic tectonic provinces. The Gawler Craton spans across nearly the entire South Australian Crustal element. c Outline of the workflow applied in this study to obtain targeting maps from surface and subsurface datasets. The first step is to categorise the datasets based on their penetration depth into surface and subsurface signals. Lineaments are then extracted with the methods outlined in sections 3.3, 3.4, and 3.5.
Note that the multi-scale edge detection ('worming') is only applied to the magnetic and gravity data, manual segmentation was only performed on the DEM, and automatic lineament mapping with the PCI Geomatica LINE module was applied to all datasets. We performed a geometric analysis on the extracted lineaments to highlight the variability between the different methods ( Figure 11). Targeting maps are derived from computed line density and intersection density maps. The first step for obtaining these maps is to merge lineaments into a single dataset if the lineaments are extracted with the same method. This is performed for surface and subsurface lineaments respectively, i.e. extracted from the DEM and the radiometrics with the LINE module of the PCI Geomatica software are merged into a single dateset representing the collated lineaments attributed to changes in surface topology and chemical composition. A summary of the merged data utilised is shown in Table 1. The density maps are then computed for datasets that comprise surface and subsurface lineaments (see Table 2).
We visualise the targeting maps as colour coded line density maps that are contoured by the intersection density. Based on a defined threshold, potential targeting areas are then highlighted in these maps. A more detailed workflow diagram is included in the appendix Figure A1.
signatures, whereas Proterozoic granites comprise both zoned and massive ovoid plutons of low and high magnetisation. For the Hitaba Suite, Schmidt and Clark (2011) already pointed out the high variability in airborne magnetic signature.
Overlying the crystalline basement of the Gawler Craton in the analysed area are Palaeozoic, Mesozoic and Cenozoic 125 sedimentary sequences that combined form significant but variably thick cover ranging in thickness from 0 to more than 600 m. However, the spatial distribution of each sedimentary sequence is poorly understood (Hou, 2004).
The terrain across most of the study area is relatively flat to moderately undulated. Prominent topographic highs are localised around dissected rocky outcrops. The surface is characterised either by aeolian sand covering variably weathered bedrock mainly in elevated parts of the landscape, or by saline playa lakes and drainage tracing topographic lows. One distinct feature The oldest preserved cover is composed of Late Carboniferous to Early Permian post-glacial sediments within the Mulgathing Trough in the northwestern corner of the study area ( Figure 2) (Nelson, 1976;Hibburt, 1995). Magnetic source depth analysis suggests that the trough is more than 600 m deep, with the maximum depths likely not detected .

135
The overall variability of the cover thickness in the study region is shown in Figure 4.

Structural Framework
An interpretation of the geological framework was undertaken around Tarcoola (Wilson et al., 2018;Pawley and Wilson, 2019) using the new GCAS aeromagnetic data (see Figure 5c) that was collected between 2017 and 2018, at 60 m ground clearance with 200 m spacing on east-west flight lines . The dominant linear structures that can be identified in the 140 aeromagnetic data are shear zones and faults. Identified faults and shears tend to be relatively discrete zones whose magnetic signature was altered by circulating fluids, or by the juxtaposition of two blocks of rock with different magnetic characters. In general, the shear zones are prominent on aeromagnetic images as they form extensive structures that often separate lithological packages with contrasting magnetic characters or are associated with changing trends of the magnetic grain (e.g. Figure 5c).
Faults often form shorter, narrow features with changing magnetic expressions that can be difficult to recognise in rocks with a 145 low magnetic response. For details on the signature of fault in aeromagnetic data we refer to Grauch and Hudson (2007). The area in Figure 2 can be divided into a series of structural domains with distinct fault patterns.
Few faults are identified in the Christie Domain, although this could be due to the bland magnetically low character of the region. One exception is the >80 km-long, northwest-trending Mulgathing Trough (Figure 2b), which can be recognised as it is filled with Permian glacial sediments that affect its magnetic response.

150
The Wilgena Domain contains northwest-trending faults that are particularly prominent as narrow non-magnetic zones in the magnetic Hiltaba Granite plutons ( Figure 2b). Some faults are relatively straight to curvilinear and can be traced for more than 80 km (e.g. the Lake Labyrinth Fault), whereas others are shorter and form anastomosing and bifurcating structures. The northwest-trending faults typically show apparent dextral offset, and usually cut the major shear zones. An exception to this trend is the north-northeast-trending Tarcoola Fault that appears to propagate from the southern Finke Shear Zone.

155
The faults form several patterns in the Nuyts Domain ( Figure 2). To the northwest, the faults are northwest-trending with dextral offset. The eastern Nuyts Domain is characterised by northeast-trending sinistral faults, aligned sub-parallel to the Koonibba Shear Zone located to the southwest of Figure 2 (see detailed map of (González-Álvarez et al., 2020)). In the central Nuyts Domain, a pluton from the Hiltaba Suite has a long, straight northwest-trending margin that is bound by the Kooniba The NW-and NE-trending faults across a large portion of the study area cut the ≈ 1585 Ma Hiltaba Suite plutons. Therefore, the faults formed at or after ≈ 1585 Ma, likely during Kararan Orogeny that occurred either at 1570 Ma (Hand et al., 2007) or 1600-1570 Ma (Reid et al., 2017, and references therein). There is evidence that pre-existing structures (i.e. Gulgunnia, Muckanippie, and Coorabie shear zones (see Figure 2) were reactivated during the Kararan Orogeny (Direen et al., 2005;Reid 165 and Dutch, 2015). In conclusion, most of the large-scale fault zones could have provided pathways for mineralising fluids.
Based on the observable offsets of pluton contacts a N-S directed shortening can be assumed during the formation of the Au-deposit mentioned above.

Methodology
In the following section, we introduce the datasets used for lineament extraction, then describe the lineament analysis we 170 employ to compare lineaments extracted by different techniques, followed by a description of the three lineament extraction techniques used.

3.1 Datasets
Lineaments extracted from a subset (region 9A (Childara)) of the Gawler Craton Airborne Survey (GCAS), the world's largest high-resolution airborne geophysical and terrain imaging program at 200 m line spacing, were analysed by Foss et al. (2019).

175
The data that is released under the Creative Commons Attribution 4.0 International Licence includes total magnetic intensity (TMI), radiometrics (RAD), and digital elevation model (DEM).
We utilised the DEM derived from laser altimeter subtracted from differential GPS heights (Figure 5a), radiometric data (total dose) processed using the Noise Adjusted Singular Value Decomposition (NASVD) (Hovgaard and Grasty, 1997) (Figure 5b), total magnetic intensity reduced to pole (Figure 5c), and gravimetric data gridded to 100 m with a station spacing between 180 50 m and 50,000 m ( Figure 5d).
The data presented is freely available through the South Australian Resources Information Geoserver (SARIG).

Lineament analysis
Lineament analysis allows for obtaining unbiased metrics for comparison of the data in terms of their dominant strike directions. For each dataset the principal orientations are automatically obtained. Assuming that the strikes represents a multimodal 185 distribution, we first calculate the kernel density estimation (kde) using normal distributed kernels. The kde is a non-parametric estimator that basically smoothes each data point into a small density function based on the underlying kernel and bandwidth.
The obtained kde function then represents the sum of all these subfunctions. By this we obtain an smoothed estimation of the distributions probability density function. This analysis was performed using the FracG software (Kelka and Westerlund, 2021). For a general overview of how kde works we refer to Chen (2017). 190 For simplification, we assume that the individual principal strike directions follow Gaussian distributions. In this case, we obtain a best-fit model allowing for up to ten principal strike directions (Gaussians) fitted to the obtained kde. The parameters of the Gaussians are obtained via maximum likelihood fitting. Goodness-of-fit testing is performed iteratively based on a modified Akaike Information Criteria (Akaike, 1998). The amplitudes of the Gaussians are normalised and are proportional to the number of lineaments that belong to this distribution. 195 Figure 6 shows the principal orientations of the structural interpretation ( Figure 2a). Raw kernel density estimates are plotted as a dotted line to which the Gaussian model is fit to. Clearly, two perpendicular directions dominate the data with a subordinate set of roughly N-S striking lineaments.

Manual lineament extraction
Lineaments are pattern breaks within data that the human eye can depict as a straight or somewhat-curved feature in an 200 image (Boucher, 1997). This is dependent on the person's visual ability as well as technical experience and hence mapping the presence and location of surface lineaments can vary significantly between individuals. By applying different types of preprocessing (e.g. edge detection filtering, hill shading etc.) different features in the raster image can be enhanced thus leading to different lineament sets segmented from the same data set. Direct observation-based surface lineament mapping has been  Figure 7 shows the manually mapped lineaments in the DEM (a) and in the raster representing the mean gradient component (b). The mean gradient component was calculated as the arithmetic mean of the horizontal and vertical gradient components obtained through Sobel convolution filtering. This edge detection filter can utilise different kernels for enhancing 210 edges of different orientation (e.g. horizontal or vertical) (Shrivakshan and Chandrasekar, 2012). The mean gradient component we calculated as a preprocessing step prior to the manual segmentation contains enhanced edges in the horizontal and vertical orientation. The dominant orientation in both datasets is around 106 and 110 degrees, respectively. In both datasets, three Gaussians provide the best fit, whereas the two subordinate directions of both data sets differ significantly.

215
A multi-scale edge detection technique has been applied to the potential field data which produces edge features called "multiscale edges" (or colloquially "worms"). This technique (Hornby et al., 1999;Holden et al., 2000) relies on a wavelet transform based on the Green's function of vertical gravity or reduced-to-pole (RTP) total magnetic intensity. A low-resolution multi-scale edge mapping of the whole Gawler Craton was performed by Heath et al. (2009 applied a higher resolution mapping using the more recent GCAS Region 9A magnetic field data and updated gravity coverage. An fundamental part of 220 this automatic edge detection technique is the upward continuation that acts similar to a filter suppressing shallow sources (see Hornby et al. (1999) and Foss et al. (2019)). In this contribution, we selected different heights of upward continuation for the gravity and magnetic data respectively in order to derive edge maps that comprise similar detail of the subsurface.
The reasoning behind our choice is on one hand that gravity and magnetic field decay at different rates away from the source and on the other hand the datasets we utilised have different resolutions. In order to derive lineament maps that comprise a 225 comparable level of detail from the gravity and magnetic data we choose an upward continuation of 930 m for the gravity data and an upward continuation of 2070 m for the magnetic data. Note that these only represent a single example of the geophysical lineament ensembles computed presented in Foss et al. (2019). In the framework of this study, where the potential field data was utilised to derive signal from the deeper subsurface, the two chosen datasets represent a good example for this automated lineament mapping technique.

230
The potential field data was processed using upward continuation (UC) to generate edge features that can be considered representative of different depths. Upward continuation suppresses high frequencies in the data, and increases the weighting of signals from deeper physical property contrasts. Calculation of edge enhancement transforms at different upward continua- tion heights, producing a series of edge mappings ('multi-scale edges' or 'worms'). Edges derived from the gravity data map subsurface density contrasts, and those from the magnetic field data map subsurface magnetisation contrasts. The edges (in 235 particular the shallow edges) depend considerably on data distribution which is very regular for the magnetic field data, but highly irregular for the gravity data. In areas of sparse gravity coverage, it is not possible to map detail in the shallow gravity multi-scale edges. In compensation, the gravity field better expresses contributions from deeper property contrasts than the magnetic field. However, the principle value of having multi-scale edges derived from both gravity and magnetic field data is that they map quite separate physical properties, even though both properties depend on lithology. In some cases the contact 240 between two lithologies is both a density and magnetisation contrast, and the two multi-scale edge vectors are strongly correlated, but in other cases a lithology contact may cause only a significant density contrast or only a significant magnetisation contrast, giving rise to edge vectors in only one of the fields. The combination of the two sets of edge vectors is therefore much more informative than either one alone. By their nature potential fields measured above a physical property interface are automatically smoother than the trace of that interface. They cannot include abrupt changes of trends and at higher upward 245 continuations the potential field and multi-scale edge expression of any straight-line property contrast becomes progressively more curved. There are therefore compromises in matching naturally curved multi-scale edges with corresponding straight lineaments extracted from the same dataset. The principal orientation of the lineaments is roughly E-W for the gravity and magnetic data (Figure 8). The lineaments derived from the magnetic data ( Figure 8a) exhibit three main directions. In contrast, the strike of the gravity lineaments (Figure 8b) are uniform with only one clear principal orientation.

Automatic lineament extraction
Lineaments have also been extracted from the DEM, radiometrics, gravity, and magnetics using PCI Geomatica's LINE function (Geomatics, 2005). This technique relies on properties inherent to the image (e.g. pixel intensity) making use of Canny edge detection (Shrivakshan and Chandrasekar, 2012) as the basis.
The gradient of an image is computed and pixels that are not a local maximum are suppressed. Edge strength threshold of 255 pixels produces a binary image that, after applying a thinning algorithm, results in skeleton curves that represent edges. If a curve meets a minimum length criterion, the curve is approximated by line segments within an error threshold. Lineaments are the result of linking line segments if they have similar orientation. PCI Geomatics defines a lineament as a "straight or somewhat-curved feature" (Geomatics, 2005) and the parameters control the extent to which edges detected in the image may result in a line feature. Originally intended to be used on radar images, the 260 technique has been widely used in various remote sensing applications (e.g., Pandey and Sharma, 2019). Edges identified relate to significant changes within a given image and the resulting lineaments are highly dependent on user-specified parameters that control length and segment linkage. For a detailed description of the parameters used in this study we refer to González-Álvarez et al. (2020) and Pawley et al. (2021).
The lineaments were automatically extracted for datasets representative of surface ( Figure 9) and subsurface ( Figure 10  In contrast to the uniform distribution of lineament orientation obtained for the surface layers, the automatically detected subsurface lineaments exhibit two nearly equal principal directions for the total magnetic intensity ( Figure 10a) and a more uniform distribution obtained for the gravity data (Figure 10b). The dominant direction in each subsurface dataset differs significantly and are oriented nearly perpendicular to each other.
270 Figure 10. Lineaments automatically extracted with PCI Geomatica LINE module from the total magnetic intensity (reduced to pole) (a) and the gridded gravity data (b). To the right of the maps the rose diagrams show the orientation distribution (bin size of 10 degrees) and the models fitted to the probability density of the empirical data.

Comparison of lineament datasets
The lineament datasets from the three techniques applied differ the most in terms of their length distribution (Figure 11a).
In particular, the automatically generated surface lineaments exhibit a narrow distribution with the highest density limited to regions around the mean and median. The length distributions of the lineaments automatically extracted from the subsurface data show a wider range slightly skewed towards smaller values. It is worth noting that the length tolerance of lineaments is an 275 input parameter and the bias is well represented in the resulting lineament datasets. The automatic gradient extraction ('worms') yields more distributed length of lineaments where smaller lineaments dominate the data. Manually extracted lineaments show a similar length distribution independent of whether the interpretation was performed on the processed or unprocessed elevation data. The structural interpretation that was performed mainly on the total magnetic intensity data exhibits the widest distribution. We note that the different techniques will yield variable results and different information can therefore be obtained 280 from the same dataset by applying for instance the worming (Figure 8) or automatic segmentation ( Figure 10). We do not seek to compare the geological or physical information inherent to each dataset but rather perform a statistical analysis to point out the most striking geometrical differences (e.g., length and orientation). An in-depth evaluation of which extraction methods yields more reliable or realistic geological information is beyond the scope of this study. The principal strikes exhibit a common E-W trend in the surface and subsurface datasets (Figure 11a, b). The dominant 285 orientations of the surface lineaments (Figure 11b) scatter around the E-W plane with only one subordinate orientation that is somewhat oriented perpendicular (the line labelled 4: automatically extracted from radiometric data). The orientations of the subsurface lineaments are more diverse but also scatter mainly around the E-W plane. Subordinate orientations are more common in the subsurface lineament datasets and most pronounced in the worms (Figure 8a, b) and in the lineaments automatically extracted from the total magnetic intensity (Figure 10a). The latter exhibits a bimodal distribution that is comparable to the 290 manual structural interpretation ( Figure 6). Overall, the automatically extracted surface lineaments tend to yield uniform distributions for length and orientation whereas the automatically extracted subsurface lineaments exhibit wider length distributions and non-uniform strike directions. The automatic gradient extraction produces lineament sets with wider length distributions and orientations dominated by an E-W trend with subordinate orientations that are nearly perpendicular to the principal strike.
Apart from the manual structural interpretation, this method is the only one presented in this study that produces strongly 295 curved lineaments. Lineaments that are manually extracted show a comparable length distribution independent of whether the data was processed by edge detection filtering. The locations and orientations are influenced by the processing and so is the number of features (Figure 7a, b). In summary, a dominant orientation trend is observable in the automatically generated sur-face and subsurface lineaments that is roughly east-west. Only the manually extracted surface lineament set shown in Figure 7b and the structural interpretation ( Figure 6) exhibit considerable divergence from this orientation. The main difference between 300 the different lineament sets is demonstrated by their length distributions (Figure 10).

Lineament density maps as an exploration tool
Here we present a workflow for exploration targeting based on lineament density and intersection density per unit area that utilises remotely sensed surface and subsurface data. Lineament datasets that are obtained with the same method and correspond to signals either both from the surface or both from the subsurface are merged (see Tables 1 and 2). The resulting dataset 305 comprises the collated lineaments associated with topographical or chemical changes in the surface. We performed the same for the potential field data (gravity and magnetics) to obtain a collection of all geophysical lineaments detected by the respective extraction method. The aim of this analysis is to identify areas of maximum line and intersection densities in the combined surface and subsurface signals. The merged datasets with their name used in this section are shown in Table 1. The density maps are calculated for several combinations of surface and subsurface layers (Figure 1c, Tables 1 and 2, and Figure A1). The 310 two types of density maps for deriving potential exploration targets are: -Density maps of lineaments per area (P20) -Density maps of lineament intersections per area (I20) P20 maps represent the number of lineaments per unit area derived for rectangular sampling windows of size 2 km by 2 km.
The pixel resolution of the derived raster file is set to the search window size. I20 maps are derived by converting the lineament 315 data into a graph representation where intersections between lineaments are vertices (Sanderson et al., 2019). The number of intersections is derived using a pixel size of 2 km by 2 km and circular sampling windows with a radius of 2.5 km. The maps are then up-sampled via bilinear interpolation to a cell size of 500 m by 500 m. To obtain the density maps we used the open-source software FracG (Kelka and Westerlund, 2021).
The targeting maps are derived by overlaying the P20 maps with contours of the I20 maps (Figures 12, 13, and 14). We used 320 combined lineament data that represent surface or subsurface signals (see Table 1). In cases where we utilised the structural interpretation ( Figure 2) as the subsurface datasets the density maps comprise three individual lineament datasets (Figures 12   and 13 a) whereas the other targeting maps are comprised of four lineament sets (Figures 13 d and 14). The reason behind this is that the remotely sensed data comprise two datasets for the surface (DEM and radiometrics) and two datasets representing subsurface (TMI and gravity) respectively. Furthermore, the signals might detect features at different scales whereby the 325 sensitivity for structures at a particular scale is a function of penetration depth and spatial resolution. Especially, the gravity and magnetic data used in this study have very different resolution ( Figure 5) and will therefore yield lineaments maps with different detail. Datasets used for obtaining the respective targeting map and extraction methods used for obtaining the utilised lineaments are summarised in Table 2. We further classified the targeting maps into "manual", "automatic", or "semiautomatic"  Subsurface auto PCI Line module TMI RTP (Figure 10c) Gravity (Figure 10d) indicating whether the underlying lineament sets were derived with purely manual segmentation, represent a combination of 330 manually and automatically extraction, or are obtained solely by automatic segmentation (Table 2).
Considering not only the overall lineament density but also the intersection density allows us to further constrain potential targeting areas. By obtaining intersection densities, cross-strike features can be identified that are thought to represent zones of enhanced permeability (Wheeler, 1980;Southworth, 1985). Areas of enhanced structural complexity or numbers of crossstrike discontinuities could therefore represent zones for preferential upwelling of mineralising fluids. We suggest that the 335 adjacent areas that comprise an overall high density of lineaments represent preferential exploration targets. For identifying these mineral potential zones, we set a threshold of 9 intersections per 500 m by 500 m pixel size and then visually identified the areas of overall high densities in the vicinity of these specific points as favourable targeting areas. The threshold is kept constant across datasets in this study to ensure a better comparability but would need to be adjusted depending on the underlying data for more reliable targeting.
340 Figure 12. a Targeting map derived from "Surface manual" combined with the structural interpretation ( Figure 2a). b Lineament density map (P20). c Intersection density map (I20). Potential targeting areas are indicated by red ellipses. Figure 13. a Targeting map derived from "Surface auto" and the structural interpretation. b Lineament density map (P20). c Intersection density map (I20). d Targeting map derived from "surface manual" and the "subsurface worm"). e Lineament density map (P20). f Intersection density map (I20). Potential targeting areas that are indicated by red ellipses.

Discussion
This is the first study to utilise the newly acquired high-resolution datasets of South Australia (GCAS Region 9A (Childara)) to investigate the applicability of lineament mapping/extraction from a variety of data sets as a potential exploration tool. As both manual and automatic lineament extraction methods are subject to bias, we first compared the results obtained from both approaches.

345
During manual interpretation bias arises from subjectivity that is introduced by different interpreters (Raghavan et al., 1993) and can also be caused by scale or variable processing techniques for edge enhancement such as illumination azimuth (Scheiber 19 Figure 14. a Targeting map derived from "surface auto" and the "subsurface worms". b Lineament density map (P20). c Intersection density map (I20). d Targeting map derived from "surface auto" and the "subsurface auto"). e Lineament density map (P20). f Intersection density map (I20). Potential targeting areas that are located at the margins of the domain indicated by red ellipses. et al., 2015; Masoud and Koike, 2017). We do not seek to explore these aspects, but do note that interpretation bias may be playing a role in creation of the manual lineament datasets used. Automatic mapping methods can also be subjected to bias related to the type of applied edge detection filter and underlying segmentation algorithm. We applied the lineament extraction 350 algorithm from the commercial software PCI Geomatica, as one of the main topics of this study is comparing different conventional methods with automatic detection methods. Some input parameters of the PCI Geomaticas LINE module will certainly introduce a bias towards certain length distributions of the extracted lineaments. The two most important parameters are the threshold to length and the threshold to angular difference that defines the maximum difference in orientations for uniting Table 2. Summary of the targeting datasets (Figures 12, 13, and 14) and name convention used in the combined targeting map (Figure 15).

Targeting datasets
Name ( two segments (see Table 7 in Kelka and Martinez (2020) for details on the parameters). We note that recently several new 355 approaches were suggested for automatic lineament mapping (e.g., Zhang et al., 2006;Hashim et al., 2013;Mohammadpour et al., 2020;Xu et al., 2020) and utilising one of these might yield results different to the ones presented here.
We will first discuss the similarities and differences of the surface and subsurface lineaments sets obtained with the different methods pointing out the individual strength or weaknesses. Lineaments obtained by manual mapping in this study scatter around a particular length ( Figure 11). The length distribution of automatically extracted surface lineaments is even narrower, 360 pointing towards a bias in lineament detection potentially related to the parameter combination applied for automatic mapping (see section 3.5). One other major difference is that the manual mappings yield three main directions for each dataset whereas the automatic detection exhibits uniform distributions. The principal orientation of automatically and manually extracted lineaments are similar for the radiometrics and DEM, respectively, but differ by a maximum of 18 • between the methods when compared to the dominant directions obtained from the manually derived lineaments.

365
The automatic detection method yields data that is more representative for local scale geomorphological features visible by the strong influence of the sand ridges in the southwestern part. In contrast, a human interpreter tends to identify the general trends in the data where results can be biased by the preprocessing of the data such as edge enhancement filtering. In summary, a general trend of superficial features extracted automatically and manually scatter around a E-W to NNW-SSE direction and therefore this orientation is likely a characteristic regional feature of the cover across the investigated area. Whether 370 automatic lineament extraction is superior to manual mapping cannot be stated with certainty based on the data presented in this study. Choosing one method over the other might depend on the type and the thickness of cover. Generally human investigators identify more regional trends, which are probably hard to detect automatically. The reason behind this is that a human interpreter will consciously detect the trends of lineaments and will merge them even if there are large gaps between the identified edges.

375
Geophysical lineaments are obtained from the gravity and magnetic data via automatic gradient extraction. The lineaments obtained from each geophysical dataset differ significantly in terms of length distribution and principal orientation ( Figure 8).
This could be attributed to different depths the upward continuation represent but is more likely to reflect the difference in resolution of the data and the physical properties each dataset is sensitive to. The gravity data yields lineaments that are attributed to major lithological boundaries pronounced by density contrasts. In the study region the prominent boundaries are 380 the margins of the domains that often coincide with large crustal-scale shears and the Mulgathing Trough in the northwest (see Figures 2a, b). While the major crustal scale elements are traced by the lineaments extracted from the gravity data, the lineaments obtained from the magnetic data seem to reveal a more detailed picture of the subsurface structural framework.
In addition to domain boundaries the magnetic lineaments also outline large intrusive bodies. In line with the the structural interpretation, three main directions are detectable for the magnetic "worms" (Figures 2 and 6) but only two for the gravity 385 "worms" (Figure 8). While the edge vector's orientation and length distributions differ significantly, a reasonable correlation considering their locations is observable (Figure 8). In line with Foss et al. (2019), this suggests that the mapping of gravity and magnetic contrasts with worms allows for correlating the magnetic and gravity field anomalies. It must be pointed out that the mapped edges only act as approximate markers of the horizontal contrasts in density or magnetisation . We note that utilising different values for the upward continuation will yield lineament maps comprising different detail and the 390 datasets shown in this study represent only a single example from the range of upward continuations used in Foss et al. (2019).
We note that using different values will alter the results and represents a source of uncertainty. However, for the purpose of this study the two upward continuation values were chosen to allow a more reliable comparison between the different datasets.
The presence of magnetic remanence may alter the field anomaly, rendering the reduction-to-pole data we used less useful.
However, for our purposes of extracting lineaments from multiple datasets, the uncertainty in the degree of magnetic remanence 395 is of less concern compared to the uncertainty associated with the different automated and manual techniques in extracting lineaments. We note that the magnetic and gravity datasets are of different resolution and in particular the resolution of the gravity dataset is non-uniform. As the upward continuation acts similar to a low-pass filter the difference in resolution becomes negligible allowing for a more reliable comparison of the potential filed data used in this study.
Automatic lineament mapping performed with PCI Geomatica (Geomatics, 2005) yields a picture less consistent with the 400 structural interpretation of the subsurface framework. While the automatically extracted gravity lineaments still outline some crustal-scale boundaries (especially evident for the graben structure in the northwest) the information associated with the automatically extracted geophysical lineaments is inferior compared to the information that can be obtained by automatic gradient extraction and the latter method should be favoured for the mapping of geophysical lineaments. The "worms" trace the subsurface in greater detail and profound physical meaning can easily be attributed to the location of the lineaments as they 405 are associated with strong lithology contrasts within the basement units. In contrast, the automatically extracted lineaments pick up a rough impression of the structural framework of the subsurface with only major elements detectable in the data, such as major shear zones and the Permian graben. We conclude that automatic gradient extraction is the superior technique for extracting geophysical lineaments from high resolution magnetic data and from gravity data with variable resolution.
In section 5 we tested an approach that integrates surface and subsurface lineaments in a simple framework for exploration 410 targeting that is based on identifying areas of high lineament density and high intersection density. The justification for this approach is that welling of hydrothermal fluids is often associated with structurally complex zones that comprise a high intersection density (e.g. Dimmen et al., 2017). Areas comprising an overall high density of discontinuities that are adjacent to such zones of high structural complexity represent preferential exploration targets for hydrothermal mineralisation. By combining surface and subsurface datasets we not only account for intersections in subsurface datasets but also for intersection of 415 subsurface and surface lineaments. Such cross-strike discontinuities are an additional indicator for structurally complex zones and are taken into account in our workflow.  Figure 15). Here it is important to note that the gold deposits in the central Gawler Craton exhibit some similar characteristics considering the mineralisation style, as the gold is dominantly hosted in sulphide-poor structurally controlled quartz veins that seem to be spatially related to the Hiltabas Suite (Daly, 1993). We do not directly identify the deposit exploited at the Tarcoola mine but an area to the northeast that is situated along the margin of the Tarcoola formation (known host-rock (Pawley and Uncertainty in manual lineament mapping is directly related to the person's experience and the scale they are intending on mapping. The manual extraction of lineaments in this study focused on the regional linear trends (lineament greater than 1 km).
Addressing uncertainty for the automatic lineament mapping is difficult and directly related to the resolution of the underlying datasets. In case of the automatic gradient extraction, the upward continuation can pose another source of uncertainty related 435 to the loss of detail that increases with higher upward continuations. For an example of how uncertainty in lineament mapping can be assesses statistically, we refer to Pawley et al. (2021). For a reliable interpretation of the obtained lineament maps the geological history of the respective area must be considered. This means that differently oriented lineament sets could correspond to different tectonic and fluid flow events. In areas comprising multiphase deformation some extracted lineaments might therefore not be of relevance for the targeted mineral system and, for instance, directional constraints on the utilised 440 lineaments would need to be applied. In the case study presented here, the youngest orogenic event (Kararan Orogeny) is thought to be linked to the mineralisation Bockmann et al., 2019) and to the reactivation of pre-existing structures (Direen et al., 2005;Reid and Dutch, 2015). In this special case there is no need for applying strict constrains on the extracted lineament sets but this might be necessary for other regions.
We note that relatively little known mineralisation coincides with the targets identified by the lineament analysis ( Figure 15) 445 and further research is needed to validate the reliability of the presented workflow. Geological knowledge of the area might help to reduce the number of false positives obtained by lineament-based exploration targeting.

Conclusion
In this study we pointed out the differences between subsurface and surface lineaments in the Gawler Craton in South Australia mapped/extracted with different methods and from a variety of remotely sensed and geophysical data. We determined the 450 principal orientations of each dataset by automatically deriving a best-fit Gaussian model of the data. Overall an E-W direction dominates in surface and subsurface datasets that likely represent the structural grain of the area. Surface lineaments manually mapped are clearly subjective and can be biased due to preprocessing of the data. Compared to automatic extraction the main difference seems to be the scale on which the extracted lineaments play a role; the manual interpretation picks up regional scale trends whereas automatically extracted features represent smaller scale, locally relevant structures. We found that the 455 automatic gradient extraction ("worming") is superior to automatic lineament extraction performed with PCI Geomatica as the worms detect more details that are related to lithological and structural contrast. In this study we showed that automatic gradient extraction yields geophysical lineaments associated with a profound geological meaning compared to automatically detected edges. In terms of the surface lineaments we conclude that the automatic extraction represents a method that picks up more local scale features and is likely well suited for well exposed areas. However, in areas that comprise a thick, reworked 460 cover manual mapping of lineaments yield a more regional scale picture and seems to represent a reliable method.
An integrated workflow that utilises surface and subsurface lineaments should include density per unit area and intersection density per unit area. We found that a combination of geophysical lineaments derived by automatic gradient extraction combined with either manually or automatically mapped surface lineaments represents the most promising combination of data for exploration targeting. The gravity and magnetic worms will coincide with present lithological boundaries or major structural 465 features. To clearly state whether edges or lineaments observable in the surface data are correlated with crustal scale features such as shear zones requires further research. For efficiently combining intersection density and line density for targeting, spatial clustering algorithms might yield more reliable results compared to the simple approach presented in this study. The crucial parameters will be setting an appropriate threshold for intersection and line density for determining target areas.
Code and data availability. Datasets and the code for automatic lineament analysis are freely available:
Author contributions. Ulrich Kelka wrote the manuscript with input from all authors, analysed the data, and developed the computational framework. Cercia Martinez wrote parts of the manuscript, analysed data, performed the automatic extraction of lineaments, and revised the manuscript. Carmen Krapf devised the project, wrote parts of the manuscript, performed the manual lineaments extraction, and revised 475 the manuscript. Stefan Westerlund developed the computational framework, helped with data visualisation, and revised the manuscript.
Ignacio Gonzalez-Alvarez devised the project, helped with data interpretation, and revised the manuscript. Mark Pawley wrote parts of the manuscript, performed the structural interpretation, helped with data interpretation, and revised the manuscript. Clive Foss performed the automatic gradient extraction for the geophysical datasets, helped with interpretation, and revised the manuscript We thank the four anonymous reviewers for their thorough revision which helped to improve this manuscript.