Ambient seismic noise analysis of LARGE-N data for mineral exploration in the Central Erzgebirge, Germany
Ambient seismic noise tomography is a novel, low-impact method for investigating the earth's structure. Although most passive seismic studies focus on structures on crustal scales, there are only a few examples of this technique being applied in a mineral exploration context. In this study, we performed an ambient seismic experiment to ascertain the relationship between the shallow shear wave velocity and mineralized zones in the Erzgebirge in Germany, one of the most important metal provinces in Europe. Late Variscan mineralized greisen and veins occurring in the Geyer-Ehrenfriedersdorf mining district of the Central Erzgebirge were mined from medieval times until the end of the 19th century. These occurrences represent a significant resource for commodities of high economic importance, such as tin, tungsten, zinc, indium, bismuth and lithium. Based on ambient noise data from a dense “LARGE-N” network comprising 400 low-power, short-period seismic stations, we applied an innovative tomographic inversion technique based on Bayesian statistics (transdimensional, hierarchical Monte Carlo search with Markov Chains using a Metropolis/Hastings sampler) to derive a three-dimensional shear wave velocity model. An auxiliary 3D airborne time-domain electromagnetic dataset is used to provide additional insight into the subsurface architecture of the area. The velocity model shows distinct anomalies down to approximately 500 m depth that correspond to known geological features of the study area, such as (a) gneiss intercalations in the mica schist-dominated host rock, imaged by a SW–NE striking low-velocity zone with a moderately steep northerly dip, and (b) a NW-trending strike-slip fault, imaged as a subvertical linear zone cross-cutting and offsetting this low-velocity domain. Similar to the velocity data, the electromagnetic data exhibit north-dipping (high-conductivity) structures in the mica schists, corresponding to the strike and dip of the predominant metamorphic fabric. An unsupervised classification performed on the bivariate 3D dataset yielded nine spatially coherent classes, one of which shows a high correspondence with drilled greisen occurrences in the roof zone of a granite pluton. The relatively high mean shear velocity and resistivity values of this class could be explained by changes in density and composition during greisen formation, as observed in other areas of the Erzgebirge. Our study demonstrates the great potential of the cost-efficient and low-impact ambient noise technology for mineral exploration, especially when combined with other independent geophysical datasets.
Passive seismic methods using the ambient seismic noise field have become very popular over the last few years for imaging the earth's subsurface structure. In particular, the discovery of deriving the Green's functions between many seismic stations by data cross correlation turned out to be a game changer in seismology (Campillo and Paul, 2003). These Green's functions correspond to a “virtual source” located at any given seismic station being recorded at all other stations, thus allowing the velocity structure (and, in turn, the geological structure) to be inferred without earthquake signals or controlled sources necessary for traditional imaging methods (Bensen et al., 2007). The use of ambient noise fields (from wind, traffic, sea waves etc.) instead of signals from controlled sources (vibroseis, explosions) has several advantages: data acquisition is more environmentally friendly, logistically less demanding and comparatively cheap. Passive seismic can be applied to a wide range of spatial scales and can yield similar near-surface spatial resolutions to refraction seismic imaging. Recent advances in electronic components made the development of low-power recording systems possible, which are used for deployments of large/dense, temporary networks (“LARGE-N”, formed by hundreds to thousands of nodes). Although ambient seismic noise analysis has successfully been used in geodynamic investigations (e.g. Green et al., 2020), shallow site characterization (e.g. Hannemann et al., 2014) and geothermal site characterization (Martins et al. 2020), applications in mineral exploration are relatively sparse (e.g. Hollis et al., 2018; Da Col et al., 2020; Xie et al., 2021; Colombero et al., 2021).
In this paper we describe a passive seismic survey that has been conducted to study technical and methodological aspects of passive seismic methods as well as their usefulness for the application in an exploration context.
A worldwide growing demand in raw materials, largely in response to the widespread deployment of green energy technologies, electrification of transport and digitalization (e.g. Sovacool et al., 2020), has led to renewed exploration activity in the Erzgebirge metallogenic belt (e.g. Dittrich et al., 2020; Burisch et al., 2021), which is well-endowed with high-tech metals such as Sn, W, Zn, In, Ge, Bi and Li (Štemprok and Seltmann, 1994). The study area in the Geyer-Ehrenfriedersdorf mining district (Central Erzgebirge) was selected given the local occurrence of greisen and vein deposits with significant potential for future ore discoveries. Because of the physicochemical processes associated with ore formation, ore bodies may differ in their petrophysical properties from the surrounding rocks. Hence, through observation and interpretation of variations in geophysical properties, areas of ore formation can be delineated (Parasnis, 1986; Telford et al., 1991), complementing data from costly field and drilling-based exploration work. The geophysical method used for the delineation of mineral deposits needs to be adapted to the geological setting and type of mineralization. In this case, tin and lithium-bearing greisen and vein deposits similar to those occurring in the Geyer-Ehrenfriedersdorf area have been documented to differ from the surrounding host rocks in terms of their seismic velocities and electrical conductivities (Müller-Huber and Börner, 2017). However, the region is densely populated and contains environmentally sensitive areas, placing additional constraints on the chosen mineral exploration techniques. Ambient seismic and airborne electromagnetics are considered environmentally sensitive and socially acceptable exploration technologies that are well suited in this context. The study is carried out within the framework of the EU-funded INFACT project (2017–2021), which sought to establish infrastructures and protocols for the trialling of exploration technology. Being located within an INFACT reference site, this study benefits from a large portfolio of exploration data for validation and integration purposes. The results of this study, apart from increasing the robustness of the benchmark targets of the INFACT reference site, will contribute to improving targeting of polymetallic greisen-type mineralization.
2.1 Study area, petrophysical constraints and data acquisition
In July/August 2020 we carried out a passive seismic experiment at the INFACT reference site in Geyer-Ehrenfriedersdorf, Erzgebirge, SE Germany. The study area measures 1.0×1.7 km, has an elevation of 570–710 m a.s.l. and is located in a dense forest between the villages of Ehrenfriedersdorf to the SE and Geyer to the SW (Fig. 1). Underneath a thin (on average approximately 3.5 m, locally thickening to 10 m) veneer of Quaternary sediments, the Geyer granite pluton of late Variscan age (325–318 Ma, Förster and Romer, 2010) underlies the area, emplaced into a suite of Proterozoic to Early Palaeozoic metamorphic rocks ranging from high-pressure two-mica schists with muscovite gneiss intercalations in the SE to phyllites in the NW. The Greifensteine, a locally renowned rock formation just outside the study area to the NW (Fig. 1a), forms the exposed top of this granite, and, based on drilling evidence, it is known to drop off steeply due SE to a depth of approximately 400 m below the surface in the SE part of the study area (Fig. 1b). The Geyer granite is a medium-grained monzogranite that belongs to a suite of high-F, high-P2O5 Li-mica granites in the Central Erzgebirge (Förster et al., 1999). It is locally associated with porphyritic microgranite dykes. The dominant Late Palaeozoic metamorphic fabric trends SW–NE and generally dips NW, with the dip angle increasing from east to west (Hösel et al., 1994). Post-orogenic tectonics led to the formation of a series of steep, NW–SE and NE–SW striking faults, including the prominent Greifenbach fault, which is associated with variable lateral offsets of up to 100 m (Fig. 1a).
The Erzgebirge is endowed in world-class magmatic-hydrothermal ore deposits hosted by a variety of vein and metasomatic structures that have been extensively mined for silver, tin, bismuth and cobalt since the Bronze Age (Tolksdorf et al., 2019). The most important mineralized structures in the study region are closely associated with the emplacement of the Geyer granite, and include stockwork-like greisen bodies and veins occurring in the endo- and exo-contact of the granite, respectively (Hösel et al., 1994). Greisen are metasomatic granitic rocks composed of quartz, mica and topaz, and are a common host of tin, tungsten and lithium ores in the Erzgebirge (Štemprok and Seltmann, 1994; Dittrich et al., 2020). According to Jung and Seifert (1996) the (Sn-W) mineralization occurred in three stages resulting in a zonal distribution with respect to the distance from the granite: (i) a first stage of greisenization is associated with the formation of lithium micas and topaz in the endo-contact of the pluton as 2–10 thick peripheral zones and apophyses, (ii) a second stage of greisenization led to the formation of a Sn-W-As paragenesis with cassiterite, wolframite and arsenopyrite as a vein-shaped exo-contact mineralization of 0.1 to 1 m thickness (i.e. Sn-greisen veins in Fig. 1a), and (iii) a late-stage hydrothermal stage involved fluorite and sulfide precipitation. The Sn-W-As succession contains a variety of commodities such as Sn and W and other metals of economic importance such as In, Cu, Co, Bi, Sb and Au. Although largely eroded or mined out near the surface, exploration drilling in the study area from about 1960 to 1990 has documented the remaining mineralized greisen occurrences at depth, with reserves of stringer and greisen ores in the Ehrenfriedersdorf mining district being calculated at 17 100 t (Hösel et al., 1994). Although the endo-contact greisen mineralization is considered impoverished in terms of Sn-W content (Hösel et al., 1994; Jung and Seifert, 1996), it is known to contain lithium, which, being the dominant battery technology for electric vehicles, is currently in high demand worldwide. Li-Sn-W greisen deposits in the nearby Zinnwald/Cínovec area are currently being developed into an operating mine (Dittrich et al., 2020).
Direct information on surface wave velocity is not available for the rocks of the study area, but petrophysical information from a variety of sources that relate to this parameter can be used to support the geological interpretation. For instance, bulk density is a parameter that critically influences the seismic velocity in rocks (e.g. Salah et al., 2018). Based on a report by Klemps and Lindner (1985), the two-mica schist in the Geyer-Ehrenfriedersdorf area has a median density of 2.76 g cm−3 (hornfels variety up to 2.80 g cm−3), muscovite gneiss a median of 2.65 g cm−3 and granite a median of 2.62 g cm−3 (Fig. 2a). Greisen in the study area exhibit nominally higher densities (median = 2.69 g cm−3) than the granite. Mineralogical changes induced by greisen alteration of granitic rocks are known to cause an increase in grain density owing to the metasomatic replacement of less dense minerals such as albite or K-feldspars by denser minerals such as mica and topaz during greisenization. However, greisenization also induces a porosity increase owing to the molar volume decrease associated with mineral replacement reactions (Launay et al., 2019). Density values given in Klemps and Lindner (1985) represent water-saturated density and not bulk density, i.e. they only partially include porosity (from sealed-off or occluded pores), so the bulk density of granite and greisen might not be different if greisen had a higher porosity than granite. Borehole neutron gamma logs from drill holes in the study area (see representative log in Fig. 2b), which are sensitive to the amount of hydrogen atoms in a formation and are commonly used as a proxy for porosity, show no discernible difference between greisen and their host granite (according to a report by Sonntag, 1977, summarizing borehole geophysical data from a variety of drill holes in the area). This indicates that the porosity that initially developed in the greisen facies may have been reduced by subsequent mineral precipitation (Launay et al., 2019). Greisen occurrences of similar age and composition from the nearby Altenberg area in the Erzgebirge also exhibit low porosity values (median = 2.63 %; Müller-Huber and Börner, 2017), even lower than their granite host rocks (median = 6.7 %). In terms of electrical resistivity, the primary lithologies two-mica schist, muscovite gneiss and granite cannot be differentiated according to a summary report for borehole geophysical data from a drilling campaign just north of the study area (Tita, 1978). Variations in resistivity are instead attributed to differences in the degree of fracturing, alteration and mineralization. There is no reference to greisen in this report, but in the nearby Altenberg area, the greisen characteristically exhibits higher specific electrical resistivities (∼ 3250 versus ∼ 1485 Ωm) than their granite host rocks (Müller-Huber and Börner, 2017).
A total of 400 seismic stations were deployed to record ambient seismic noise continuously for 10 d. The instruments were installed in an area of approximately 1700×1100 m on a dense, regular grid with ∼70 m spacing between the grid nodes. We used DIGOS/Omnirecs digital CUBE recorders equipped with 1- and 3-component seismic sensors (4.5 Hz Eigenfrequency). Accurate timing of the recordings was realized through cycled GPS, recording at 400 samples per second.
2.2 Data processing and inversion
Ambient noise techniques have been applied to recover the Green's functions between all stations (80 000) by data cross-correlation (Campillo and Paul, 2003; Wapenaar, 2004; Curtis et al., 2006). As mentioned before, these Green's functions correspond to virtual sources at a given station recorded at all other stations. To derive Green's functions, we applied the data processing to the vertical component noise recordings described in Bensen et al. (2007): it consists mainly of splitting the data records into daily segments, one-bit normalization, spectral whitening and cross-correlation of the traces. The resulting daily cross-correlations are finally stacked. As the vertical components of all receivers are identical, we did not remove instrument responses. The Green's functions are typically dominated by dispersive surface waves (Rayleigh waves).
Before calculating the Green's functions we analysed the characteristics of the ambient noise, i.e. the frequency content and its directional properties. We checked the directional properties of the noise by plain-wave beam forming in the time domain (see for example Schweitzer et al., 2012; Harjes and Henger, 1973). For the analysis we used a subset of 22 stations and calculated the beam power of 10-min time windows after bandpass filtering (1.2–4.5, 4.5–10, 10–20 Hz). Six consecutive beam-forming results were stacked to yield hourly beam-forming plots. Figure 3 shows examples of these beam-forming plots (normalized beam power) together with a spectrogram (power spectral density; single station) for an arbitrary day (31 July 2020). Noise levels vary significantly with time. Low-frequency noise is often coming from north-north-westerly directions (apparent velocities 2.5–4 km s−1); however, other directions are also represented. In the higher-frequency bands the noise is also coming from a broad range of directions and slownesses; however, some rather constant sources can also be identified (e.g. signal around 10 Hz in the spectrogram of Fig. 3 coming from ∼ NW with an apparent velocity of around 6.6 km s−1).
The virtual shot gathers constructed from the vertical component recordings are shown in Fig. 4a for all station pairs and Fig. 4b for a single virtual shot gather. The latter gather shows a significant asymmetry of the causal and acausal phases.
In the next step the Rayleigh dispersion curves were determined in an interactive way for a large number of cross-correlations (FTAN method from Dziewonski et al., 1969, with band = 0.25 and β = 3.15, GUI tool implementation by Ryberg et al., 2021a). To avoid near-source effects we limited our dispersion curve analysis to traces with offsets >1 km and further decimated the data set to ∼ 5000 station pairs, by requiring a minimum spatial receiver station separation. We derived dispersion curves at 30 frequencies ranging from 1.2 to 20 Hz in steps equally spaced with log(frequency). We started picking for every trace when an energetic arrival was present at frequencies between 2 and 5 Hz, from there extending the picking to lower and higher frequencies until a coherent (continuous) dispersion curve branch disappeared. We used the instantaneous frequency to avoid systematic dispersion curve shifts (caused by non-uniform frequency distribution of energy in the trace) for the following inversion step. Figure 5 shows two typical examples of a dispersion curve for an individual station pair. Despite the fact that a number of traces suffer from poor signal-to-noise ratio, because all stations had been deployed in a densely wooded region, ∼ 45 000 travel time picks (frequency dependent travel times of group velocity arrivals) could be determined. This subset represents ∼ 7 % of the dispersion curves of the entire data set (Fig. 4a, 400×400 stations2 = 80 000 station pairs). Bootstrapping inversion tests using an even further decimated data set resulted in very similar final 3D models, so our 7 % subset appears to be quite representative. Finally, we used this travel time (or dispersion curve) data set to invert for a three-dimensional shear-wave velocity model below the seismic array.
Instead of using a classical linearized inversion with a Tikhonov-type regularization (including damping, smoothing, a starting model etc.) we applied a tomographic inversion technique based on a Bayesian statistical method (transdimensional, hierarchical Monte Carlo search with Markov chains using a Metropolis/Hastings sampler, Bodin et al., 2012) to derive the 3D distribution of shear wave velocity (and its uncertainty), liberating us from the potentially subjective choice of starting velocity models, damping and smoothing parameters etc. that are necessary when using classical inversion techniques. We follow the fully 3D technique by Zhang et al. (2018) to invert for a 3D velocity model, but replace the involved travel time calculations by a fast finite-difference based eikonal solver (Podvin and Lecomte, 1991). The latter, given its high computational speed, is essential when using the Markov chain Monte Carlo inversion approach as the forward problem (i.e. calculating the travel times for a given model) has to be solved orders of magnitude more often than classical inversions.
The 3D seismic velocity field is discretized by a set of polyhedral Voronoi cells, with a shear wave velocity value assigned to every cell. The model misfit, i.e. the differences between the travel time observed and predicted by the model, is calculated using the following steps: the 3D model is horizontally gridded (50-m spacing), at every location a vertical 1D shear velocity model is extracted from the 3D velocity model, and then the group velocity dispersion curve is calculated (Herrman and Ammon, 2004; Herrmann, 2013) for the 1D model, assuming Vp Vs ratios and density values according to Brocher (2005). When repeated for all locations of the horizontal grid, 2D velocity maps for all frequencies are constructed. For all 2D velocity maps, travel times (with respect to the virtual sources and receivers) are calculated using a fast eikonal solver (Podvin and Lecomte, 1991) and eventually the misfit is determined.
The Markov chains are started with a model consisting of a random number of Voronoi cells with random cell centres and shear wave values. By randomly changing the model (changing the velocity and/or position of a cell, adding or removing a cell etc., details in Zhang et al., 2018) we construct a Markov chain of consecutive models. After a burn-in phase, the (well-fitting) models are decimated (i.e. only every 200th model is selected from all 15 000 well-fitting models) gathered for further derivation of a reference solution. For even more extensive model space exploration, models from 1000 separately evolving chains are investigated, which is easily achieved on computer clusters. The final results, after the burn-in phase (several thousands of well-fitting models) are combined to analyse their statistical properties, i.e. averages and uncertainty.
The inversion is transdimensional because the number of Voronoi cells is not fixed, and it is hierarchical because we invert for data noise, i.e. all remaining data misfit, which cannot be explained by the model. This approach has the advantage that the derived models (distribution of velocities and their uncertainties) are almost completely data driven, i.e. no prior choice of damping, smoothing, starting model etc. is needed.
Figure 6 shows the evolution of 1000 Markov chains, both for misfit (panel b) and model complexity (panel a, number of Voronoi cells). The advantage of using a hierarchical version of the Markov chain Monte Carlo method resulted in the derivation of the (remaining) data noise for every signal frequency, i.e. the part of the data that could not be explained by the models. This noise is caused by mispicks of the dispersion curves, anisotropy (not considered in the inversion), noisy Green's functions etc. This noise is typically increased at the upper and lower parts of the frequencies with a lower number of data points, thus potentially with poor signal-to-noise ratios, as can be seen in Fig. 7.
2.3 Electromagnetic data
For geological interpretation, the passive seismic data are supplemented by resistivity data based on an airborne electromagnetic survey using the Geotech versatile time domain (VTEM™ ET = Versatile Time-Domain Electromagnetic - Early Time; Eadie et al., 2017) system. The VTEM™ ET transmitter-receiver loop was in concentric-coplanar and z-direction oriented configuration and was towed at a mean distance of 30 m below a helicopter. The survey was flown in a south to north, 40∘ NW direction, with a traverse line spacing of 50 m.
Fifty-four time measurement gates were used for the final data processing in the range from 5.79 to 9481 µs. The full waveform EM(electromagnetic)-specific data processing operations included half cycle stacking (performed at the time of acquisition), system response correction, as well as parasitic and drift removal. In addition, a three-stage digital filtering process was applied to reject major sferic events and the signal-to-noise ratio was further improved by the application of a low-pass linear digital filter.
The VTEM ET data were reprocessed to eliminate coupling with mapped infrastructures, to increase , and assign noise. The data were inverted with spatial constraints in a layered earth approach (SCI). Where present, induced polarization effects were taken into account, using the dispersive resistivity model of Cole and Cole (see Viezzoli et al., 2017). The result is a 3D distribution of Cole–Cole parameters, the most relevant being direct current resistivity and chargeability.
2.4 Data integration
In order to facilitate a comparison between the VTEM resistivity data and the seismic velocities, the inverted EM data, represented by a layer thickness with a corresponding resistivity value, was discretized along each 1D resistivity-depth array to a resolution of 2 m and cropped to the depth of investigation. Subsequently, Bentley Systems Leapfrog Geo (version 5.1.1) was used to create a numerical model for the VTEM point data at 10 m resolution using spheroidal interpolation with spatial constraints including a 25-m distance buffer and topography. Similarly, the seismic velocity data with an original point spacing of 10 m were interpolated using a linear interpolant using the same spatial constraints as the VTEM data. Finally, both parameters were mapped into one and the same block model with -m cells.
To systematically evaluate the spatial distribution and geological significance of specific value ranges within this integrated block model, an unsupervised classification was performed on the bivariate 3D dataset (Fig. 8). The vector of variables at each grid cell was extended with its neighbourhood to account for the spatial structure of the data (spatial clustering, see Talebi et al., 2020). Two different neighbourhoods were considered: one containing the cross of 6 grid cells immediately in contact with the target grid cell, and one containing the 26 cells of a cube surrounding the target. Additionally, both spectral clustering (Ng et al., 2002; von Luxburg, 2007) and conventional k-means were applied to both neighbourhood settings. Finally, a nine group classification based on k-means with 26 immediate neighbours was chosen because it provided a well-balanced distribution of data into clusters with the largest spatial coherence. k-means was chosen as a segmentation algorithm rather than a density-based clustering approach not to exclude value ranges that are only represented by a few points.
3.1 Ambient noise inversion results
Instead of running a single Markov chain, we simultaneously investigated 1000 chains in parallel. All of them had a different, randomly chosen starting model. The evolution of the models along every chain (Fig. 6) shows that the data misfit quickly decreases, and reaches stable values after ∼ 3000 models, whereas the model complexity (= number of Voronoi cells to describe the individual model) still increases. The latter starts to stabilize from model number ∼ 5000. From model number 10 000 we assumed that the burn-in phase was finished, and stable systematic model sampling in the model space was reached. All models beyond 10 000 have been collected and analysed by calculating average S-wave velocities and their uncertainty at every location in the 3D model. The velocity model is generally characterized by strong velocity anomalies (Fig. 9a–c, g–h). It is a rather smooth model with gradual velocity variations, as it is typical of tomographic inversions (i.e. not showing first-order discontinuities). The velocity uncertainty (Fig. 9d–f), as expected, is quite low within shallow depth ranges, and quickly reaches higher values (red) at greater depths, a typical behaviour when inverting surface wave data. This uncertainty (velocity resolution) has to be taken into account when interpreting the recovered S-wave velocities.
The shear wave velocity distribution shows strong lateral contrasts and is well resolved down to approximately 500 m below the surface. The following velocity anomalies can be delineated:
SA1: The seismic velocity model from the surface down to a depth of approximately 400 m is marked by relatively high shear velocities in the NW portion, and lower shear velocities in the SE portion of the area, giving rise to a prominent NW–SE spatial trend. The differences between these regions are expressed most strongly at the northern end of the study area, where adjacent high- and low-velocity zones form a pronounced SW–NE striking, vertical discontinuity (anomaly SA1, Fig. 9b). This discontinuity aligns well with the dominant tectonic grain of the area and the SW–NE striking, steep geological contact between two-mica schist hornfels within the contact metamorphic aureole of the granite (high velocities) in the NW and mica schist and muscovite gneiss (lower velocities) due SE (see Fig. 1).
SA2: In the southern part of the study area, the prominent NE–SW spatial trend described above (SA1) is crosscut (and offset?) by a NW-striking linear zone, extending from the surface to a depth of about 100 m (anomaly SA2, Fig. 9a, h). This zone coincides with a known, regionally important fault zone (the Greifenbach fault, see Fig. 1a). It also marks the trend of the Greifenbach stream and associated Quaternary deposits.
SA3: Two lobe-shaped, NW-dipping low-velocity zones occur in the two-mica schist at depths down to approximately 400 m (anomaly SA3, Fig. 9c, g). These lobes are orientated parallel to the main metamorphic fabric and are also imaged by resistivity data (anomaly EA2).
Interestingly, the seismic data only broadly image the granite–mica schist interface, and we observe both high- and low-velocity zones cutting obliquely across this lithological boundary. Thus, the observed velocity distribution is likely related to factors that are superimposed on the primary lithology, such as weathering, alteration and/or mineralization.
3.2 Electromagnetic data
The VTEM airborne electromagnetic partially overlaps with the area of the passive seismic survey (Fig. 10a), and provides valuable complementary information. In depth, the resistivity model is constrained by the depth of investigation, which varies from approximately 50 to 600 m, depending on geology, EM system attitude, noise in the data and inversion stability. Within the study area, i.e. the bounds of the seismic array, we observe the following resistivity anomalies (Fig. 10):
EA1: A conductive top layer that shows a good correlation with the mapped extent of Quaternary sediments and cover thickness extracted from drill cores (stippled lines in Fig. 10b–d).
EA2: NNE–SSW striking high-conductivity zones in the eastern part of the study area, corresponding to the location of muscovite gneiss intercalations in the two-mica schists. The 20–30∘ NW dip of these anomalies follows the approximate trend of the main foliation (Fig. 1). EA2 is also imaged by the passive seismic data (anomaly SA3).
EA3: Four to five distinct high-conductivity zones, circular in horizontal section (Fig. 10a), and plume-shaped in cross-section (Fig. 10c). These anomalies are related to anthropogenic features, i.e. game fences to protect trees in forest cultures, confirmed by field validation. The presence or absence of an anomaly at the fence is due to the fence characteristics, e.g. density and quality of the fence stranding (metal weaving), whether grounded to the earth, orientation, size, metal quality etc.
3.3 Integration and geological significance
The identified anomalies in velocity and resistivity space provide valuable insights into the subsurface structure of the study area. To investigate anomalies in the bivariate space and facilitate geological interpretation, particularly in terms of mineralization, a subdivision of value ranges of the combined 3D dataset using a spatially constrained K-means clustering approach is applied. The resulting nine classes are spatially relatively coherent (Figs. 11c, 12c) and integrate both the shear velocity and the resistivity structure. For example, the classification clearly reproduces the gently NW-dipping low-resistivity zones (Fig. 12, rows 1 and 4) as well as the ± vertical contacts between zones of high and low shear velocities (Fig. 12, rows 1, 3, 6 and 7). The classes exhibit distinct ranges of shear velocities and resistivities (Figs. 8b, 13a), some of which can be associated with identified anomalies in velocity and/or resistivity space. Cluster 1 has some of the highest resistivities and seems to be spatially associated with two-mica schists on the SW side of the Greifenbach fault. Cluster 2 has the lowest values of both resistivity and shear velocity and corresponds to the seismic anomaly SA3 and resistivity anomaly EA2. It exhibits a prominent NW dip and, together with the less conductive cluster 6, most likely corresponds to the muscovite gneiss intercalations in the two-mica schists. Cluster 3 corresponds with the seismic anomaly SA2 and resistivity anomaly EA1. It is mostly confined to the surface and thickening along the Greifenbach fault, and is therefore interpreted to be associated with faulting-related fracturing and Quaternary deposits. Cluster 4 exhibits the highest shear velocities of all clusters and occurs mostly in the apical part of the granite on the NE side of the Greifenbach fault. Cluster 5 spatially coincides with the mapped two-mica schist within the contact metamorphic aureole of the granite. Cluster 7 exhibits intermediate resistivity and shear velocity values, and, occurring mostly in the endo-contact part of the granite, may represent the normal, non-mineralized granite together with cluster 4. Cluster 8 is characterized by intermediate shear velocity and high resistivities, and is located at the granite interface, where it spatially coincides with most of the drilled greisen occurrences in the considered volume (Figs. 11c, 12c, 13). Cluster 9 is characterized by relatively low velocities and resistivities and spatially coincides with the anthropogenic plume-like conductive EM anomalies EA3 (Figs. 10c, 11b, 12b).
A number of clusters (1, 3, 4, 5 and 8) coincide with drilled occurrences of greisenized rocks (Fig. 13b). The cluster spatially most closely associated with the drilled greisen occurrences is cluster 8, which has a width of up to 750 m, a length of 1350 m and a thickness of 200 m (Figs. 11, 12). Its lateral dimension is constrained to a certain degree by the model boundaries, the chosen classification parameters and the inherent smoothness of the 3D inversion models. The mean shear velocities and resistivities of this cluster (average = 2.4 km s−1, 3.8 Ωm) are relatively high compared with those of the classes representing the surrounding non-greisen rocks (classes 2, 4, 6 and 7, average = 2.1 km s−1, 2.9 Ωm). This is in line with the petrophysical constraints outlined above. Although speculative given the described uncertainties, our data thus hint at the presence of a large greisenized zone at the periphery of the Geyer granite at depths of 50–250 m underneath the surface. Such a body is currently not documented, but greisen bodies of similar dimensions to those of cluster 8 are known from the Ehrenfriedersdorf area at the Sauberg and Vierung prospects (Brosig et al., 2020). In light of the fact that the Erzgebirge is considered underexplored with regard to its Li mineral resources (Dittrich et al., 2020), the stage-1 greisen zones in the study area, which are considered to be barren with respect to Sn-W (Hösel et al., 1994; Jung and Seifert, 1996), may simply not have received much attention.
The geophysical exploration of greisen ores in the Geyer-Ehrenfriedersdorf area is a challenge because of their occurrence in a populated area, their petrophysical indistinctiveness and their limited size. As demonstrated in this study, some of these challenges can be overcome by using a combination of different low-impact geophysical methods complemented by geological and petrophysical data. Our analysis demonstrates the potential of 3D passive seismic data from dense networks for the geological characterization of the subsurface in a crystalline rock environment. Passive seismic data are particularly well suited to application in populated areas, because of their low environmental impact and the fact that they are unperturbed by anthropogenic electromagnetic noise, which is a problem for magnetic and electromagnetic prospecting (e.g. Morris et al., 2021). Nonetheless, as shown in this study, airborne time-domain electromagnetic data, which were used to identify the cover–basement interface and conductive layers in the metamorphic host rocks, can provide useful indirect constraints as complementary information for the delineation of greisen mineralization.
We conducted an ambient noise seismic study in the Geyer area to gain insight into the 3D structure of the subsurface and evaluate the region's potential for greisen-hosted, polymetallic mineralization. Based on a novel 3D inversion technique that has been applied to a passive seismic dataset, supplemented by and integrated with 3D resistivity data from airborne time-domain electromagnetics we are able to:
identify distinct shear velocity anomalies down to a depth of approximately 500 m that most likely correspond to primary geological features such as the transition from mica-schist hornfels to muscovite gneiss dominated host rocks and near-surface fracturing related to strike-slip faulting;
use resistivity anomalies to delineate the cover–basement interface as well as conductive zones parallel to the main fabric of the metamorphic host rocks;
employ spatial clustering of the combined seismic-EM dataset and drill core validation to identify parameter ranges corresponding to a potential greisen mineralization.
The study showcases the benefits of a combined use of low-impact, geophysical data collection and processing technologies for mineral resource exploration. Crystalline settings (as in the Geyer-Ehrenfriedersdorf case) are particularly challenging for high-resolution controlled-source seismic. The resolution of the ambient noise velocity models is lower than those expected for reflection seismic imaging in a sedimentary setting but passive surveys can provide very valuable constraints at reasonable costs for an enhanced understanding of subsurface architecture and the distribution of ore bodies, particularly if integrated with other geophysical parameters.
The study area turned out to be a well-suited test area because of the large amount of information on the subsurface owing to the existence of numerous archive data that allowed the benchmarking of the results and a geologically consistent interpretation.
The passive seismic data set was gathered with a large network of regularly spaced seismic sensors, but only a small part of the entire data set was used to derive a 3D velocity model. By using other subsets of the entire data set, optimized experiment designs, i.e. number and geometric distribution of seismic sensors, could be developed and evaluated for future studies.
Seismological data are archived at the GEOFON/EIDA archive (https://doi.org/10.14470/KS7575632014, Haberland et al., 2021). The velocity-resistivity 3D block model is available under Ryberg et al. (2021b, https://doi.org/10.14278/rodare.1222) and viewable under https://bit.ly/3B2aEvD (last access: 7 March 2022).
TR designed the field methodology of passive seismic investigations, developed the inversion software, performed the passive seismic inversions, participated in the field work, and contributed to the original article preparation. MK contributed to concept development, participated in field work, conducted the integrative, geological analysis and substantially contributed to writing the paper. CH contributed to the inversion methodology, participated in the field work, data curation and contributed to the original article preparation. RTD coded and performed the spatial clustering analysis and contributed to writing the paper. AV performed QC, processing and inversion of the electromagnetic data, and contributed to the preparation of the article. RG assumed management and coordination responsibilities for this project, helped to design the experiments, acquired part of the funding that led to this publication and contributed to editing of the paper.
The contact author has declared that neither they nor their co-authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This article is part of the special issue “State of the art in mineral exploration”. It is a result of the EGU General Assembly 2020, 3–8 May 2020.
The inversion of seismic data was carried out on the high-performance computing cluster at GFZ. We gratefully acknowledge Geotech Ltd. for the acquisition and pre-processing of the EM data. Sam Thiele and Sean Walker are gratefully acknowledged for their help with data processing. Thanks also to Leïla Ajjabou, Yuleika Madriz, Erik Herrmann, Hernan Flores, Thomas Lüttke, Juan-Felipe Bustos, Ariane Siebert, Falk Brethauer, Matthias Bosdorf and Michael Weber for their help in the passive seismic field campaign. Moritz Kirsch acknowledges Tobias Duteloff (Sächsisches Landesamt für Umwelt, Landwirtschaft und Geologie) for his kind assistance and provision of archive material. Furthermore, we thank the city of Ehrenfriedersdorf for their support and hospitality during field work.
We would like to thank Nicholas T. Arndt and the anonymous reviewer for providing comments that helped to improve the article.
Funding was provided by the GFZ and by the European Union's Horizon 2020 research and innovation programme under grant agreement no. 776487. Instruments for the seismic network were provided by the Geophysical
Instrument Pool Potsdam (GIPP, GFZ), grant GIPP202010.
The article processing charges for this open-access publication were covered by the Helmholtz Centre Potsdam – GFZ German Research Centre for Geosciences.
This paper was edited by Juan Alcalde and reviewed by Nicholas T. Arndt and one anonymous referee.
Bensen, G. D., Ritzwoller, M. H., Barmin, M. P., Levshin, A. L., Lin, L., Moschetti, M. P., Shapiro, N. M., and Yang, Y.: Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int., 169, 1239–1260, https://doi.org/10.1111/j.1365-246X.2007.03374.x, 2007.
Bodin, T., Sambridge, M., Rawlinson, N., and Arroucau, P.: Transdimensional tomography with unknown data noise, Geophys. J. Int., 189, 1536–1556, https://doi.org/10.1111/j.1365-246X.2012.05414.x, 2012.
Brocher, T. M.: Empirical relations between elastic wavespeeds and density in the Earth's crust, B. Seismol. Soc. Am., 95, 2081–2092, https://doi.org/10.1785/0120050077, 2005.
Brosig, A., Barth, A., Knobloch, A., and Dickmayer, E.: Rohstoffprognosen für Zinn, Wolfram, Fluss- und Schwerspat im Mittelerzgebirge, Bergbau in Sachsen, Band 19, Bergbaumonografie, Sächsisches Landesamt für Umwelt, Landwirtschaft und Geologie (LfULG), 2020.
Burisch, M., Frenzel, M., Seibel, H., Gruber, A., Oelze, M., Pfänder, J. A., Sanchez-Garrido, C., and Gutzmer, J.: Li-Co–Ni-Mn-(REE) veins of the Western Erzgebirge, Geminimum layer thicknessrmany – a potential source of battery raw materials, Miner. Deposita, 56, 1223–1238, https://doi.org/10.1007/s00126-021-01061-4, 2021.
Campillo, M. and Paul, A.: Long-range correlations in the diffuse seismic coda, Science, 299, 547–549, 2003.
Colombero, C., Papadopoulou, M., Kauti, T., Skyttä, P., Koivisto, E., Savolainen, M., and Socco, L. V.: Surface-wave tomography for mineral exploration: a successful combination of passive and active data (Siilinjärvi phosphorus mine, Finland), Solid Earth Discuss. [preprint], https://doi.org/10.5194/se-2021-121, in review, 2021.
Curtis, A., Gerstoft, P., Sato, H., Snieder, R., and Wapenaar, K.: Seismic interferometry – Turning noise into signal, Leading Edge, 25, 1082–1092, 2006.
Da Col, F., Papadopoulou, M., Koivisto, E., Sito, L., Savolainen, M., and Socco, L. V.: Application of surface-wave tomography to mineral exploration: A case study from Siilinjärvi, Finland, Geophys. Prospect., 68, 254–269, https://doi.org/10.1111/1365-2478.12903, 2020.
Dittrich, T., Helbig, M., Kühn, K., Bock, W.-D., and Müller, A.: The Zinnwald Lithium Project: Transferring legacy exploration data into new mineral resources, European Geologist Journal, 49, 12–18, 2020.
Dziewonski, A., Bloch, S., and Landisman, M.: A technique for the analysis of transient seismic signals, B. Seismol. Soc. Am., 59, 427–444, 1969.
Eadie, T., Legault, J. M., and Plastow, G.: Introducing VTEM-ET: An Improved Helicopter Time-Domain EM System for High Resolution Applications, 1–5, Extended Abstracts – 15th SAGA Biennial Conference & Exhibition 2017, https://doi.org/10.1071/ASEG2018abW9_3H, 2017.
Förster, H. J. and Romer, R. L.: Carboniferous magmatism, in: Pre-Mesozoic Geology of Saxo-Thuringia – From the Cadomian Active Margin to the Variscan Orogen, edited by: Linnemann, U. and Romer, R. L., Schweizerbart, Stuttgart, 287–308, ISBN 978-3-51065-2-594, 2010.
Förster, H.-J., Davis, J. C., and Tischendorf, G.: Multivariate analyses of Erzgebirge granite and rhyolite composition: implications for classification of granites and their genetic relations, Comput. Geosci., 25, 533–546, https://doi.org/10.1016/s0098-3004(98)00171-x, 1999.
Green, R. G., Sens-Schönfelder, C., Shapiro, N., Koulakov, I., Tilmann, F., Dreiling, J., Luehr, B., Jakovlev, A., Abkadyrov, I., Droznin, D., and Gordeev, E.: Magmatic and sedimentary structure beneath the Klyuchevskoy volcanic group, Kamchatka, from ambient noise tomography, J. Geophys. Res.-Sol. Ea., 125, e2019JB018900, https://doi.org/10.1029/2019JB018900, 2020.
Haberland, C., Ryberg, T., and Kirsch, M.: The GEYER seismic network, GFZ Data Services, Other/Seismic Network [data set], https://doi.org/10.14470/KS7575632014, 2021.
Hannemann, K., Papazachos, C., Ohrnberger, M., Savvaidis, A., Anthymidis, M., and Lontsi, A. M.: Three-dimensional shallow structure from high-frequency ambient noise tomography: New results for the Mygdonia basin-Euroseistest area, northern Greece, J. Geophys. Res.-Sol. Ea., 119, 4979–4999, https://doi.org/10.1002/2013JB010914, 2014.
Harjes, H.-P. and Henger, M.: Array-Seismologie, Zeitschrift für Geophysik, 39, 865–905, 1973.
Herrmann, R. B.: Computer programs in seismology: an evolving tool for instruction and research, Seismol. Res. Lett., 84, 1081–1088, 2013.
Herrmann, R. B. and Ammon, C. J.: Computer Programs in Seismology: Surface Waves, Receiver Functions, and Crustal Structure, Version 3.30, Dept. of Earth and Atmospheric Sciences, Saint Louis University, available at: https://www.eas.slu.edu/eqc/eqccps.html (last access: 7 March 2022), 2004.
Hollis, D., McBride, J., Good, D., Arndt, N., Brenguier, F., and Olivier, G.: Use of ambient-noise surface-wave tomography in mineral resource exploration and evaluation, SEG Technical Program Expanded Abstracts 2018, https://doi.org/10.1190/segam2018-2998476.1 2018.
Hösel, G., Hoth, K., Jung, D., Leonhardt, D., Mann, M., Meyer, H., and Tägl, U.: Das Zinnerz-Lagerstättengebiet Ehrenfriedersdorf/Erzgebirge, Bergbaumonographie, 1, 189 pp., 1994.
Jung, D. and Seifert, T.: On the metallogeny of the late Hercynian tin deposit “Röhrenbohrer field”/Greifenstein area, Sn-W district Ehrenfriedersdorf-Geyer, Erzgebirge, Germany, Freiberger Forschungshefte C, 467, 131–150, 1996.
Klemps, P.-J. and Lindner, K.: Dokumentationsbericht Ergebnisse – Petrophysik Ehrenfriedersdorf, VEB Geophysik Leipzig, Leipzig, 1985, 25 pp., 1985.
Kirsch, M. and Steffen, M.: Erstellung eines geologischen 3D-Modells für das Pilotprojekt ROHSA 3.1 (Rohstoffe Sachsen), unpublished scientific report, Helmholtz-Institute Freiberg for Resource Technology, 84 pp., 2017.
Launay, G., Sizaret, S., Guillou-Frottier, L., Fauguerolles, C., Champallier, R., and Gloaguen, E.: Dynamic Permeability Related to Greisenization Reactions in Sn-W Ore Deposits: Quantitative Petrophysical and Experimental Evidence, Geofluids, 2019, 1–23, https://doi.org/10.1155/2019/5976545, 2019.
Martins, J. E., Weemstra, C., Ruigrok, E., Verdel, A., Jousset, P., and Hersir, G. P.: 3D S-wave velocity imaging of Reykjanes Peninsula high-enthalpy geothermal fields with ambient-noise tomography, J. Volcanol. Geoth. Res., 391, 106685, https://doi.org/10.1016/j.jvolgeores.2019.106685, 2020.
Morris, W. A., Ugalde, H., Kirsch, M., Gloaguen, R., Ibs-von Seht, M., Siemon, B., and Meyer, U.: A multi-parameter approach for recognition of anthropogenic noise in aeromagnetic data collected over populated areas: Erzgebirge, Germany, Geophys. Prospect., https://doi.org/10.1111/1365-2478.13179, online first, 2021.
Müller-Huber, E. and Börner, F.: Multi-parameter petrophysical characterization of Variscan greisen rocks from the Southern Bohemian Batholith (Austria) and the Eastern Erzgebirge Volcano-Plutonic Complex (Germany), Aust. J Earth Sc., 110, 78–100, https://doi.org/10.17738/ajes.2017.0006, 2017.
Ng, A., Jordan, M., and Weiss, Y.: On spectral clustering: analysis and an algorithm, in: Advances in Neural Information Processing Systems, edited by: Dietterich, T., Becker, S., and Ghahramani, Z., 14, 849–856, MIT Press, 2002.
Parasnis, D. S.: Principles of applied geophysics, 4th Edn., Chapman and Hall, New York, 402 pp., ISBN 978-94-009-5814-2, 1986.
Podvin, P. and Lecomte, I.: Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools, Geophys. J. Int., 105, 271–284, 1991.
Ryberg, T., Haberland, C., Green, R. G., and Sens-Schönfelder, C.: A Fast GUI-Based Tool for Group-Velocity Analysis of Surface Waves, Seismol. Res. Lett., 92, 2640–2646, https://doi.org/10.1785/0220200425, 2021a.
Ryberg, T., Kirsch, M., Haberland, C., Tolosana Delgado, R., Viezzoli, A., and Gloaguen, R.: Block model of passive seismic shear velocity and airborne electromagnetic resistivity in the Geyer area, Erzgebirge, Germany (Version 1.0), Rodare [data set], https://doi.org/10.14278/rodare.1222, 2021b.
Salah, M. K., Alqudah, M., El-Aal, A. K. A., and Barnes, C.: Effects of porosity and composition on seismic wave velocities and elastic moduli of lower cretaceous rocks, central Lebanon, Acta Geophys., 66, 867–894, 2018.
Schweitzer, J., Fyen, J., Mykkeltveit, S., Gibbons, S. J., Pirli, M., Kühn, D., and Kværna, T.: Seismic Arrays, in: New Manual of Seismological Observatory Practice 2 (NMSOP-2), edited by: Bormann, P., Deutsches GeoForschungsZentrum GFZ, Potsdam, 1–80, https://doi.org/10.2312/GFZ.NMSOP-2_ch9, 2012.
Sovacool, B. K., Ali, S., Bazilian, M., Radley, B., Nemery, B., Okatz, J., and Mulvaney, D. R.: Sustainable Minerals and Metals for a Low-Carbon Future, Science, 367, 30–33, https://doi.org/10.2139/ssrn.3513313, 2020.
Štemprok, M. and Seltmann, R.: The metallogeny of the Erzgebirge (Krušné hory), in: Metallogeny of collisional orogens, edited by: Seltmann, R., Kämpf, H., and Möller, P., Czech Geol. Survey, Praha, 61–69, ISBN 80-7075-152-5, 1994.
Talebi, H., Peeters, L. J. M., Mueller, U., Tolosana-Delgado, R., and van den Boogaart, K. G.: Towards Geostatistical Learning for the Geosciences: A Case Study in Improving the Spatial Awareness of Spectral Clustering, Math. Geosci., 52, 1035–1048, https://doi.org/10.1007/s11004-020-09867-0, 2020.
Telford, W., Geldart, L., Sheriff, R., and Keys, D.: Applied Geophysics, 2nd Edn., Cambridge University Press, New York, 792 pp., 1991.
Tita, J.: Zusammenfassender Bericht über die in den Bohrungen des Teilobjektes Röhrenbohrer-S, Ehrenfriedersdorf durchgeführten geophysikalischen Bohrlochmessungen, VEB Geophysik, DB BLM, Interpretation Gotha, GF 23-763/78, unpublished report, 1978.
Tolksdorf, J. F., Schröder, F., Petr, L., and Herbig, C.: Evidence for Bronze Age and Medieval tin placer mining in the Erzgebirge mountains, Saxony (Germany), Geoarchaeology, 35, 198–216, https://doi.org/10.1002/gea.21763, 2019.
Viezzoli, A., Kaminski, V., and Fiandaca, G.: Modeling induced polarization effects in helicopter time domain electromagnetic data: synthetic case studies, Geophysics, 82, E31–50, 2017.
von Luxburg, U.: A tutorial on spectral clustering, Stat. Comput., 17, 395–416, 2007.
Wapenaar, K.: Retrieving the elastodynamic Green's function of an arbitrary inhomogeneous medium by cross correlation, Phys. Rev. Lett., 93, 254301, https://doi.org/10.1103/PhysRevLett.93.254301, 2004.
Xie, T., Xu, T., Ai, Y., Zeng, Q., Zhang, W., and Zheng, F.: Imaging the shallow crustal velocity structure of the Qingchengzi ore field on the Liaodong Peninsula, China, with a short-period dense array using ambient noise tomography, Tectonophysics, 813, 228913, https://doi.org/10.1016/j.tecto.2021.228913, 2021.
Zhang, X., Curtis, A., Galetti, E., and de Ridder, S.: 3D Monte Carlo Surface Wave Tomography, Geophys. J. Int., 215, 1644–1658, https://doi.org/10.1093/gji/ggy362, 2018.