Near-surface structure of the Sodankylä area in Finland, obtained by passive seismic interferometry

Controlled-source seismic exploration surveys are not always possible in nature-protected areas. As an alternative, the application of passive seismic techniques in such areas can be proposed. In our study, we show results of passive seismic interferometry application for mapping the uppermost crust in the area of active mineral exploration in northern Finland. We utilize continuous seismic data acquired by the Sercel Unite wireless multichannel recording system along several profiles during XSoDEx (eXperiment of SOdankylä Deep Exploration) multidisciplinary geophysical project. The objective of XSoDEx was to obtain a structural image of the upper crust in the Sodankylä area of northern Finland in order to achieve a better understanding of the mineral system at depth. The key experiment of the project was a high-resolution seismic reflection experiment. In addition, continuous passive seismic data were acquired in parallel with reflection seismic data acquisition. Due to this, the length of passive data suitable for noise cross-correlation was limited from several hours to a couple of days. Analysis of the passive data demonstrated that dominating sources of ambient noise are non-stationary and have different origins across the XSoDEx study area. As the long data registration period and isotropic azimuthal distribution of noise sources are two major conditions for empirical Green function (EGF) extraction under the diffuse field approximation assumption, it was not possible to apply the conventional techniques of passive seismic interferometry. To find the way to obtain EGFs, we used numerical modelling in order to investigate properties of seismic noise originating from sources with different characteristics and propagating inside synthetic heterogeneous Earth models representing real geological conditions in the XSoDEx study area. The modelling demonstrated that scattering of ballistic waves on irregular shape heterogeneities, such as massive sulfides or mafic intrusions, could produce a diffused wavefield composed mainly of scattered surface waves. In our study, we show that this scattered wavefield can be used to retrieve reliable EGFs from short-term and non-stationary data using special techniques. One of the possible solutions is application of “signal-to-noise ratio stacking” (SNRS). The EGFs calculated for the XSoDEx profiles were inverted, in order to obtain S-wave velocity models down to the depth of 300 m. The obtained velocity models agree well with geological data and complement the results of reflection seismic data interpretation.

1 Introduction

Exploration of new mineral deposits is an important task because the modern world needs many types of minerals for functioning (Reid, 2011). Most of the shallow mineral deposits around the world nowadays are well known and exploration of new deep mineral deposits becomes more difficult (Vasara, 2018). That is why cost-effective exploration techniques are required. Moreover, there is the problem that application of controlled-source seismic exploration is not always possible in nature-protected areas. Particularly in the Arctic areas, non-invasive, environmentally friendly exploration is relevant. As an alternative, application of passive seismic techniques in such areas has been proposed. The main advantage of passive seismic methods is the possibility to study the subsurface in remote areas with minimum impact on to environment (Polychronopoulou et al., 2020).

Passive seismic interferometry is a cost-efficient methodology with a relatively simple setup of field experiments. This methodology allows for retrieving impulse response of a medium, called an empirical Green function (EGF), from ambient seismic noise recorded at two receivers, assuming that the noise field is diffuse (Lobkis and Weaver, 2001; Campillo and Paul, 2003). If this condition is satisfied, it is possible to retrieve both surface and body waves from seismic noise using either cross-correlation, autocorrelation, deconvolution or cross-coherence of seismic records (Wapenaar et al., 2011). As shown in Rickett and Claerbout (1999), the condition of diffuse noise field is satisfied if the noise sources are distributed isotropically around seismic recorders and the noise registration time is long enough. The methodology to retrieve EGFs, using cross-correlation or autocorrelation of ambient seismic noise, has been successfully applied in numerous studies (e.g. Shapiro and Campillo, 2004; Roux et al., 2005; Ruigrok et al., 2011; Draganov et al., 2009; Poli et al., 2012; Tibuleac and von Seggern, 2012; Wang et al., 2015; Taylor et al., 2016; Afonin et al., 2017; Oren and Nowack, 2016; Romero and Schimmel, 2018). In addition to ambient seismic noise interferometry, the coda wave interferometry was proposed (Aki and Richards, 2002; Campillo and Paul, 2003; Snieder et al., 2002; Snieder, 2006). This methodology is also based on the diffuse filed approximation and it is widely used for different purposes, such as estimating nonlinear behaviour in seismic velocity (Snieder et al., 2002), monitoring of stress changes inside the studied medium (Grêt et al., 2005, 2006) and determination of the third-order elastic constants in a complex solid (Payan et al., 2009). Numerous studies describe results of successful application of passive seismic interferometry for exploration and other applied geophysical purposes (e.g. Cheraghi et al., 2017; Roots et al., 2017; Dantas et al., 2018; Abraham and Alile, 2019; Polychronopoulou et al., 2020; Planès et al., 2020). In spite of the numerous studies, there is the problem that the conditions of the diffuse noise field from sources outside an observation area are difficult to satisfy in many practical situations. One important condition is isotropic (that is, when the waves arrive from all azimuths) and homogeneous (when the waves arrive from all azimuths and have nearly the same energy) azimuthal distribution of noise sources. In order to satisfy this condition, a long registration time is required. However, the diffuse wavefield for EGFs evaluation can exist also under other conditions. Wapenaar (2004) demonstrated that EGFs could be retrieved from cross-correlation of two recordings of a wavefield at different receiver locations at the free surface in the case when diffuse wavefield is produced by many uncorrelated sources inside the subsurface. Wapenaar and Thorbecke (2013) also considered conditions for EGF retrieval from ambient noise originating from a directional scatterer in a homogeneous embedding medium that is illuminated by a directional noise field.

In our paper, we examine application of passive seismic interferometry for the case when noise is strongly directional and the receiver array is semi-linear. For this, we use continuous seismic data recorded during the XSoDEx (eXperiment of SOdankylä Deep Exploration) project in northern Finland (Buske et al., 2019). We analyse the ambient seismic noise recorded by stations located in different subregions of the XSoDEx study area in order to understand spatial and temporal distribution of noise sources. We perform numerical modelling of propagation of seismic wavefield corresponding to identified noise sources through synthetic models that represent certain types of heterogeneities in our study area (mafic intrusions, faults, massive sulfides). We show by numerical modelling that direct waves generated by various sources are scattered on these heterogeneities and produce a diffused wavefield. Therefore, EGFs can be retrieved by cross-correlation of this wavefield recorded at different locations (Wapenaar, 2004; Wapenaar and Thorbecke, 2013; van Manen et al., 2006). For evaluation of EGFs, we apply the method of passive seismic interferometry with the so-called signal-to-noise ratio stacking (SNRS) algorithm (Afonin et al., 2019). We show results of application of passive seismic interferometry for mapping the uppermost crust in the XSoDEx study area.

2 Experiment description

The XSoDEx seismic survey was conducted in the Central Lapland Greenstone Belt in northern Finland, around the Sodankylä region (Fig. 1). The area is famous for its mineral deposits, including the operating Kittilä gold mine west from the survey area and the Kevitsa Ni–Cu mine that also has significant amounts of platinum, palladium, gold and cobalt. The seismic survey lines cross varying geology including outcrops of the Archean basement and layered mafic intrusions (Buske et al., 2019).

Figure 1 The XSoDEx seismic experiment: (a) geographical location of the study area (red rectangle); (b) the XSoDEx survey lines (coordinate system EUREF FIN TM35FIN); (c) geological map of the XSoDEx study area (Buske et al., 2019). The XSoDEx survey lines are shown with black lines. Previous seismic reflection survey lines are plotted with blue lines.

Within the XSoDEx project, Geological Survey of Finland, TU Bergakademie Freiberg and the University of Oulu acquired seismic reflection and refraction data using the Vibroseis © truck of TU Bergakademie Freiberg and partly explosive sources during July and August 2017, resulting in four seismic profiles of total length of approximately 80 km recorded along existing roads (Fig. 1): Pomokairantie line (about 37.5 km); Alaliesintie line (about 14 km); Sakatti line (about 20 km); Kuusivaarantie line (about 16 km). The seismic reflection data were recorded in a roll-along scheme by a maximum 3.6 km long spread of cabled vertical-component geophones with 10 m spacing. The seismic reflection layout was designed to map crustal structures down to a minimum of 3 km depth in detail. The seismic refraction data were recorded separately along all reflection profiles by the other set of instruments that included 60 vertical component 5 Hz geophones and 40 three-component micro electro mechanical system (MEMS) accelerometers with Sercel Unite wireless autonomous data acquisition units by the Sercel Ltd. receiver spacing for Sercel Unite wireless system was approximately 160 m, and the whole spread for recording refraction data was intended to be about 12 km long. The sampling rate was typically 250 sps, except for one high-resolution 1 km long profile. For this profile, the distance between recording units was 10 m and the sampling rate was 500 sps. The plan was to keep 75 RAUs (remote acquisition units) recording continuous data at the spread while charging the remaining 25 RAUs at the accommodation of the field team. When the Vibroseis © truck moved to a large-enough distance from the 25 receivers at the beginning of the spread, those RAUs were demobilized to be charged during the nighttime. Correspondingly, the charged 25 RAUs were deployed at the end of the line. In the field conditions, the deployment and demobilization time schedule was adopted to take into account the progress of active seismic experiment. Refraction seismic data were extracted from the continuous data during periods when active seismic source was in operation, while the passive seismic data were accumulated during periods when active seismic source was not working (nighttime, spare days). Generally, the length of time intervals for continuous passive data recording was about 8–9 h. Thus, the XSoDEx experiment provided a good opportunity to verify results of passive seismic interferometry with controlled-source seismic data, to identify limitations of this technique in areas of generally low level of high-frequency anthropogenic noise and to propose possible improvements of known techniques.

3 Ambient seismic noise in the XSoDEx study area

As the population density in northern Finland is low, this results in a low level of anthropogenic high-frequency noise in the XSoDEx study area. In addition to microseismic noise, there are several local industrial noise sources: Kevitsa mine, traffic noise from the roads and noise from waterpower plants of the Kitinen river.

For estimation and comparison of noise level for different XSoDEx profiles, we used vertical components of ambient seismic noise recorded by wireless seismic receivers with 5 Hz geophones. The data were not acquired during nighttime on Sundays, when the anthropogenic activity is minimal. The noise spectra were estimated using the whole periods of continuous data (usually about 8–9 h) in order to find averaged characteristics of ambient noise. The sources of anthropogenic noise in our study area may be considered in an approximation of quasi-stationary noise sources during considered time intervals. For example, mining activity in Kevitsa mine is not changing during long time intervals, and mining machinery, excavation, transportation of ore and blasting in the production area are producing seismic signals of similar amplitudes and frequencies for different time periods. The dams can be also considered as sources of quasi-stationary signals. Generally, the roads may be used by transport of different types, and, as a result, the traffic noise may have some temporal differences. Nevertheless, the roads in our study area are characterized by generally low traffic. For example, the traffic along all the roads, along which the XSoDEx seismic data were acquired, included several cars per day. Highway no. 4 (Fig. 1) is characterized by slightly larger traffic during weekdays, but it is still very quiet during nighttime and weekends.

As one can see in Fig. 2, the noise amplitude and its frequency spectra differ significantly at stations located at different sites in the XSoDEx area. Station V1 was installed in the area characterized by the highest noise level for all frequencies analysed. It can be explained by location of this station close to the Kevitsa mine and the dam of a waterpower plant. Seismic noise recorded by station P2 has a lower amplitude than the noise recorded by the V1 station. It can be suggested that the main noise sources at that station were the dam and highway no. 4. A narrow peak at frequencies of about 26–30 Hz is seen in the seismic noise spectrum of station V2. We suggest that the main noise source at station V2 was the dam. Unexpected results of noise level estimations were obtained for stations P2 and A. As one can see in Fig. 2, station P2 is characterized by a relatively high level of seismic noise. At the same time, seismic noise at station A is relatively low, despite similar distance from these stations to the dam and to the road. We need to remember, however, that data acquisition along different profiles was made during different time periods, so probably some additional high-frequency noise source was acting at P2 during the data acquisition period.

Figure 2 Illustration of noise level at six stations located in different profiles of the XSoDEx experiment. The positions of strong noise sources (Kevitsa mine (yellow star) and dams of waterpower stations (black triangles)) and seismic stations selected for analysis (blue circles) are indicated. Brown and white lines show roads. The inset figure shows a comparison of ambient noise power spectral density estimated at selected stations.

It is clear that dominating noise sources are different across the area of our study, and the general condition for passive seismic interferometry (the sources need to be isotopically and homogeneously distributed around the study area) is not satisfied. However, local sources of high intensity can be used for evaluation of EGF for selected profiles. According to the results of spectral analysis and a priori knowledge about locations of potential noise sources, the following possible candidates for sources of signals for passive seismic interferometry can be proposed:

1. Kevitsa mine, because all the profiles are located at distances of about 6–42 km from the mine;

2. Kitinen river and waterpower plants located on the river, because three of four profiles are located along the river;

3. the waves that are scattered on heterogeneities and can produce diffused wavefield, as proposed by van Manen et al. (2006), Wapenaar (2004), Wapenaar and Thorbecke (2013); and

4. we can also use the signals from Vibroseis © and explosions recorded in the XSoDEx refraction experiment to utilize propagation of surface waves to long offsets. Such analysis is not possible with the data acquired in typical near-vertical reflection experiments, because only short offsets and limited recording times are used. In addition, active sources have relatively high frequencies, and they can be used only for shallow subsurface investigation.

In the next chapter, we investigate the wavefield produced by these possible sources using numerical modelling.

4 Numerical modelling of seismic wavefield from different sources

There are a few previous theoretical and numerical studies of the wavefield from various sources scattered on heterogeneities of elastic properties (Aki, 1969; Wu and Aki, 1985; Frankel and Clayton, 1986; Gritto et al., 1995; Bohlen et al., 2003). They showed that the scattered wavefield can be quite complicated, depending on the shape of heterogeneity and its elastic properties, location of the receiver in the far field or near field and other factors. For simulation of seismic wavefield propagation in the bedrock typical for the XSoDEx area, we used SOFI3D software, which solves a wave equation with the finite-difference method (https://git.scc.kit.edu/GPIAG-Software/SOFI3D/tree/Release, last access: 10 December 2019). For simulation of the wavefield scattered on heterogeneities we developed synthetic model based on a priori knowledge about geological structure of the study area. These main features are subvertically oriented mafic and ultra-mafic intrusions of irregular shape inside generally felsic bedrock composed of granites, gneisses and quartzites. In some places, the bedrock is overlaid by quaternary sediments with a thickness of up to several dozen metres (Leväniemi et al., 2018; Karjalainen, 2019).

We used background velocities of V_p = 5600 m/s, V_s = 3500 m/s and density = 2650 kg/m³ corresponding to felsic rocks. The embedded vertical high-velocity bodies were representing mafic rocks intruded into the felsic rocks. The bodies were 30–150 m wide, with depths varying randomly from 60 to 600 m with the following physical properties (Fig. 3): V_p = 6500 m/s, V_s = 3700 m/s and density of 2800 kg/m³. We also assumed an uppermost 60 m thick layers representing quaternary sediments with V_p = 2000 m/s, V_s = 1200 m/s and density of 1600 kg/m³. The following elastic properties of air were used as boundary conditions of the model: V_p = 330 m/s, V_s = 0 and density of 1.25 kg/m³. As sources, we used (1) plane waves from sources located out of line in far-field area, with frequencies of 50 and 2.5 Hz and with an incidence angle of 45^∘; (2) blast with dominant frequency of 30 Hz, located in the beginning of the profile; (3) waterpower plant (stationary noise with frequency of 2.5 Hz), located in line with the profile. The assumption about major sources was made based on analysis of spectra, presented in Sect. 3, and our knowledge about locations of objects of human activities (Fig. 2). In the modelling, we used the grid size of 30 m.

Figure 3 The synthetic model used for investigation of wave propagation. Position of the seismic profile marked by the red line. Low velocity corresponding to the air is indicated at the top in blue.

Figure 4 An example of synthetic seismogram (vertical component) of a plane wave that arrived with an incidence angle of 45^∘ and propagated through the synthetic model (0–4 s). The seismogram shows the first arrivals and numerous reflections. From about 1.5 s, we can see scattered arrivals of different directions with apparent velocities of 2100–2500 m/s. Several arrivals of scattered waves are indicated by red arrows.

Figure 5 Result of analysis of scattered arrivals in synthetic seismograms with an apparent velocity of 2100–2500 m/s for frequencies 2–20 Hz (Fig. 4): (a) 3C seismogram for synthetic receiver located at distance of 2000 m from the source (Fig. 4); (b) spectrogram of the scattered arrival, recorded by the vertical component of the synthetic receiver (dashed white line illustrates dispersion properties of the scattered arrival); (c) particle motion diagrams, calculated using part of seismogram indicated by red lines in panel (a), which corresponds to scattered arrival.

The first synthetic signal was a plane wave originating from a source in the far-field area. The wave front propagated from the depth of 6000 m with the angle of 40^∘ with respect to the profile direction and arrived at the surface at the incidence angle of 45^∘. As one can see in synthetic seismogram in Fig. 4, the recorded wavefield consists of the first arrival, several reflected waves and numerous scattered waves with apparent velocities of 2100–2500 m/s corresponding to surface waves. Figure 5 shows an example of particle motion diagrams (Fig. 5c) and results of spectral analysis of these arrivals (Fig. 5b). Due to elliptical polarization and dependence of phase velocity on frequency, one can conclude that these scattered waves are surface Rayleigh waves. As these surface waves have stochastic directivity effects, superposition of them may be considered as diffused wavefield that can be used, in principle, to estimate EGFs.

In this synthetic example, we demonstrate that the diffuse wavefield consisting of low-frequency (5–20 Hz) surface waves (Rayleigh) can be produced by scattering of a high-frequency (50 Hz in our case) plane wave at velocity heterogeneities. We considered monochromatic plane wave, but in real ambient noise many frequencies are usually present; thus, the scattering would be more pronounced.

The second example simulates propagation of the signal originating from a production blast in the Kevitsa mine recorded by sensors of the Sakatti profile (Fig. 1). Figure 6 shows results of modelling of the wavefield produced by the blast and propagating in the model with stochastically distributed heterogeneities. The source function of the blast in this case is delta function and the source is located in line with the profile of seismic sensors. As seen, the wavefield consists of only direct P-wave arrival, reflected P waves, multiples of P waves, reflected from vertical heterogeneities, and surface Love waves. In that case, surface Rayleigh waves are absent from the wavefield. The scattered wavefield can be seen after approximately 2 s at offsets of 2000–3600 m.

Figure 6 Synthetic seismograms of the blast. The source was located at an offset of −300 m: (a) vertical channel Z; (b) horizontal channel X; (c) horizontal channel Y.

The third synthetic example corresponds to the direct wave continuously produced by the waterpower plant, which could also be scattered on heterogeneities and produce diffused wavefield. As positions of all the dams on Kitinen river are well known and all of them are located in line with the Sakatti profile; we used the in-line position of the source in our simulation. As an input signal, we used a real seismic signal recorded by station V1 that was located at the shortest distance from the waterpower plant (Fig. 2). The spectral-time diagram of the signal is presented in Fig. 7. As one can see, there are several ranges of frequencies with some increasing of amplitudes (about 5, 12.5 and 20–50 Hz). According to Antonovskaya et al. (2017, 2019), seismic noise generated by waterpower plants may correspond to a set of spectral peaks between 3.6 Hz and about 50 Hz. The other relatively high amplitudes could be due to production activities (transportation, excavation or others) at the Kevitsa mine. Therefore, in this case, we have a complex contribution of all these sources to the noise wavefield.

Figure 7 Spectrogram of the signal recorded by station V1 (location is shown in Fig. 2) and used as a wavefield of the source in the simulation.

Figure 8a shows synthetic seismograms of the stationary wavefield corresponding to the spectrogram presented in Fig. 7. Analysis of particle motion (Fig. 8c) shows that this stationary field consists of Rayleigh waves with apparent velocities of about 2100–2500 m/s. Figure 8b shows cross-correlations of the first trace with all other traces in Fig. 8a. These cross-correlations show also P- and S-wave arrivals.

Figure 8 Synthetic seismograms of the signal with the spectrogram presented in Fig. 7 and propagating in the synthetic model shown in Fig. 3. The virtual source is located at an offset of 0 m: (a) seismograms of vertical components; (b) cross-correlation functions of vertical components, calculated between the first trace and all other traces in panel (a); (c) particle motion diagram of Rayleigh wave part in panel (b).

Figure 9 Location of passive seismic profiles of the XSoDEx Sakatti line: (a) positions of high-resolution line (green colour) and profile of lower resolution (blue line); (b) the high-resolution line in which red triangles mark positions of virtual sources; (c) profile of lower resolution (red triangles indicate positions of virtual sources).

Results of our synthetic modelling demonstrate that the plane wave scattered at heterogeneities satisfies condition of diffuse wavefield and hence can be used to extract EGFs. The wavefield, produced by scattering of stationary signal from the waterpower plant, can also be used in the cases when receivers are deployed within the first Fresnel volume area. Usage of the diffuse wavefield produced by scattering on local heterogeneities or of the stationary wavefield from a single source will have an advantage that in both cases the long registration time necessary for obtaining isotropic azimuthal coverage of ambient noise sources is not required. However, special analysis of the continuous data would be necessary in order to extract the diffuse wavefield from the data. For this purpose, the SNRS algorithm described earlier in Afonin et al. (2019) can be used. The technique is based on the global optimization algorithm, in which the optimized objective function is a signal-to-noise ratio of an EGF, retrieved at each iteration. Maximizing the signal-to-noise ratio of the retrieved EGF is ensured by stacking only cross-correlation functions coherent with each other and corresponding to the stationary phase area. The details of the algorithm are provided in the Supplement.

5 Verifying passive seismic interferometry with the scattered wavefield using passive seismic data recorded during the XSoDEx experiment

To demonstrate application of passive seismic interferometry with the scattered diffuse wavefield, we used passive seismic data acquired in the XSoDEx experiment. We applied the SNRS algorithm for diffuse field extraction and for the EGF evaluation. From the XSoDEx lines, one is particularly suitable for such demonstration. In this short high-resolution profile (green line in Fig. 9) of total length of 1000 m, both 1C and 3C Sercel wireless units were installed at distances of 10 m. The 3C sensors were installed between 1C sensors at distances varying from 20 to 30 m. The results of passive seismic interferometry along this line can be also verified using the active source seismic data acquired along the same line and results of previous geophysical experiments and drilling in this area.

For retrieving of EGFs and further dispersion curve calculation, we used continuous passive seismic data recorded during the period of 21–23 August 2017 (about 48 h). The workflow for data analysis used in our paper is shown in Fig. 10.

Figure 10 The workflow for data analysis.

Table 1 Starting model of the medium used for inversion of dispersion curves.

Figure 11 Example of EGFs with the correspondent dispersion curve, obtained from passive seismic data for the high-resolution profile shown in Fig. 9: (a) EGFs; (b) dispersion curve, extracted by the multichannel analysis of surface waves (MASW) technique; (c) particle motion diagrams for part of EGF indicated in panel (a) by red lines.

We made no a priori assumptions about the nature and spatial distribution of noise sources. For calculation of EGFs, we applied such pre-processing procedures as removing mean and trend, spectral whitening and prefiltering by the bandpass filter of 1–100 Hz. After this, we calculated cross-correlation functions in such a way that virtual sources of impulse signal were placed in the beginning, in the middle and at the end of the profile. The choice of virtual source positions is based on the dominant wavelength of the analysed signal that is about 400 m. An example of EGFs calculated with the SNRS algorithm and the correspondent dispersion curve are presented in Fig. 11.

As seen from particle motion diagram calculated for one of the 3C sensors, the main arrival seen in EGFs corresponds to the Rayleigh wave. A good coherence of waveforms of dispersed surface wave is also seen.

For calculation of velocity models from these three dispersion curves, we used Geopsy software (http://www.geopsy.org, last access: 9 April 2020). We applied a global optimization algorithm with 500 iterations to obtain the best solution. The starting model consisted of three major layers with properties presented in Table 1. The range of each property was selected using information about physical properties of rocks in the study area available from literature (Dortman, 1992; Leväniemi, 2018; Schön, 2015).

Figure 12 Results of inversion of dispersion curves obtained from passive seismic data of the high-resolution XSoDEx profile: (a) 2-D velocity model obtained by interpolation of 1-D velocity models. Velocity colour scale is shown on the left. (b) 2-D velocity model in which possible rock types are indicated by different colours. The ranges for S-wave velocities are defined according to Dortman (1992); see also Fig. 19. Velocity colour scale is shown on the left. (c) An example of 1-D velocity model (black line) obtained by inversion of the dispersion curve in panel (d), corresponding to the distance of 200 m in panel (a); (d) dispersion curve (black line), used for inversion of the velocity model in panel (c). The colour scale corresponds to values of a misfit function during different iteration steps in the global optimization algorithm.

We calculated 1-D velocity models for every 100 m of the profile using the position of virtual sources at 0, 450 and 900 m (Fig. 11a). After that, using triangular and linear interpolation of 1-D models, we obtained a 2-D model presented in Fig. 12a. Panels (c) and (d) present results of global optimization of one selected dispersion curve. The 2-D velocity model in which possible rock types are indicated by different colours is presented in panel (b). The ranges for S-wave velocities are defined according to Dortman (1992). One interesting feature of this model is a layer with S-wave velocities of about 200–400 m/s and thickness of 20–38 m that may correspond to sediments. The thickness of sediments agrees well with the result of Åberg et al. (2017), who interpolated results of GPR (ground-penetrating radar) measurements and drilling information to obtain sedimentary thickness in this area. According to them, the sedimentary thickness in the studied area is about 25–30 m. The velocities of S waves beneath the sedimentary cover will be discussed in Sect. 6.

We applied the same technique of EGF calculation to the data with lower spatial resolution recorded along the part of Sakatti profile shown by blue in Fig. 9. Particle motion diagram (Fig. 13c) shows that evaluated EGFs contain mainly surface waves. Dispersion curves were calculated for each 500 m of the profile for obtaining 2-D velocity model. For this, the virtual sources were placed at each 500 m and cross-correlation functions were calculated between the virtual source receiver and all receivers located at distances no larger than 1000 m from the virtual source. For calculation of dispersion curves, the multichannel analysis of surface waves (MASW) technique was used.

Figure 13 Example of EGFs with correspondent dispersion curve, obtained from passive seismic data: (a) EGFs on frequency band 3–8 Hz; (b) dispersion curve, calculated from EGFs by the MASW technique; (c) particle motion diagrams of a surface wave part of an EGF indicated by red lines.

As one can see in Fig. 13, the passive seismic interferometry with the SNRS allowed us to evaluate dispersion curve for frequencies of about 3.5–7 Hz. Inversion of dispersion curves was used to obtain 1-D velocity models that were combined into a 2-D model (Fig. 14a).

Figure 14 Velocity models calculated by inversion of dispersion curves along the part of Sakatti profile marked by blue in Fig. 9. (a) S-wave velocity model obtained using passive seismic data; (b) S-wave velocity model obtained using seismic data, which contain signals produced by a controlled source; (c) differences between S-wave velocity models in panels (a) and (b).

For verification of the modelling results, we compared the velocity model in Fig. 14a with the model obtained by inversion of dispersion curves estimated from surface waves produced by scattering of signal from the controlled source for the same part of Sakatti line (Fig. 14b). We compared velocity models, because they are obtained at the last step of data processing, where all the errors from previous steps are accumulated. In addition, we noticed that the differences between EGFs and dispersion curves for these two data sets are insignificant. As seen, the velocity models reveal the same details, and the velocities are generally in good agreement. From Fig. 13b one can see that the width of error bars of dispersion curves is about 500 m/s. Differences in velocities between two 2-D models (Fig. 14) are within these limits. The differences in velocities are of the order of 100 m/s in the central part of the profile. The largest difference up to 600 m/s can be seen in the beginning of the profile (from 15 000 to 15 500 m), and it can be explained either by uncertainty in dispersion curve extraction or by inversion errors.

6 Shear-wave velocity models obtained using Vibroseis © signal scattered at heterogeneities

Surface waves recorded in active source experiments (ground roll) usually considered as unwanted signal and removed from the data during processing. However, S-wave velocity models can be obtained from Vibroseis © surface waves using MASW method (Al-Husseini et al., 1981; Mari, 1984; Gabriels et al., 1987; Park et al., 1999). In this case, the depth resolution for S-wave velocity models is limited to several metres due to short offsets and small registration time in near-vertical reflection data. As the data in the XSoDEx experiment were recorded at long offsets with the wireless equipment, such recordings can be used to obtain the S-wave velocities at larger depths.

Figure 15 Cross-correlation functions of signals produced by Vibroseis © with the data recorded by wireless equipment at large offsets in the frequency band of (a) 1–10 Hz, amplitudes normalized at maximum of each trace; (b) 20–100 Hz. The signal used for cross-correlation is shown in Fig. 16.

Figure 16 An example of Vibroseis © sweep, recorded by a wireless sensor placed near the vibrator in the XSoDEx experiment: (a) seismogram; (b) spectrogram showing the frequency content of a vibrator signal.

Figure 17 An example of EGFs with correspondent dispersion curve of a Rayleigh wave on the frequency band 5–10 Hz, obtained by passive seismic interferometry: (a) EGFs, obtained by conventional method; (b) EGFs, obtained by passive seismic interferometry with SNRS algorithm; (d) dispersion curve, extracted by MASW technique; (c) particle motion diagrams for surface wave part of EGF marked by red in panel (a).

Examples of raw shot gathers of Vibroseis © signals recorded in reflection experiment are presented in Buske et al. (2019). They used bandpass filters of 30–40–100–120 Hz to eliminate surface waves from raw reflection data. The sweep frequencies were from 10 to 170 Hz. In raw reflection data, the surface wave arrivals can be followed up to 2–3 s at rather short offsets of about 350–400 m. An example of the record section compiled from wireless recorder data deployed at 160 m inter-station distances is presented in Fig. 15. In this section, a vibrator signal, presented in Fig. 16, was correlated with the traces recorded at long offsets. In Fig. 15, one can see the first arrival of the P wave with velocity of about 5400 m/s and direct Rayleigh wave with velocity of about 880 m/s in frequency band of 20–100 Hz (Fig. 15b) that can be followed to offsets of about 1 km. In the frequency band of 1–10 Hz (Fig. 15a), the surface waves cannot be followed to long offsets. The surface wave arrivals that can be correlated appear on several traces at short offsets. As seen in Fig. 16b, the frequencies of the vibrator signal start from about 12 Hz and no lower frequencies are present.

We used the SNRS technique and continuous seismic recordings of the XSoDEx experiment to obtain EGFs for all the XSoDEx profiles. As seen in Fig. 17a, the EGFs cannot be extracted from our data using conventional passive seismic interferometry, where simple stacking of cross-correlation functions is used (Shapiro and Campillo, 2004; Bensen et al., 2007; Wapenaar et al., 2011; etc.). In all EGFs, the main phase seen is the Rayleigh wave (Fig. 17b). The EGFs were used to obtain dispersion curves and invert them using Geopsy inversion software and model parameters, presented in Table 1. Figure 18 shows the 2-D velocity models for all XSoDEx profiles to the depth of 300 m obtained by interpolation of 1-D velocity models. The same velocity models, in which different colours indicate possible rock types, are presented in Fig. 19. The ranges of S-wave velocities correspond to major rock types of the Fennoscandian Shield (Dortman, 1992). The boundaries of major lithological units from Fig. 1 are also indicated in Fig. 19. All the velocity models are shown in the same colour scale.

Figure 18 Velocity models calculated by inversion of dispersion curves, obtained by passive seismic interferometry for XSoDEx profiles shown in Fig. 1. (a) Pomokairantie; (b) Alaliesintie; (c) Sakatti; (d) Kuusivaarantie.

Figure 19 Velocity models presented in Fig. 18, in which possible rock types are indicated by different colours. The ranges of S-wave velocities correspond to major rock types of the Fennoscandian Shield (Dortman, 1992): (a) Pomokairantie; (b) Alaliesintie; (c) Sakatti; (d) Kuusivaarantie.

The S-wave velocities along the XSoDEx profiles generally vary from very low values of 200–400 m/s detected in some places up to 3200 m/s. The uppermost layer with velocities of 200–400 m/s and with thickness up to 50 m corresponds to quaternary sediments, and the boundary between this layer and the lowermost part of velocity models is indicated also by the velocity contrast in 1-D models. Independent information about thickness of sediments in our study area obtained by direct drilling (summarized in Karjalainen, 2019) shows that the thickness of quaternary sediments there is not greater than 40–50 m. That is why it can be concluded that the S-wave velocities in the range of 800–3200 m/s correspond to different types of basement rocks. As these values are generally much lower than the values of S-wave velocities for the rocks of the Fennoscandian Shield, obtained by laboratory measurements on rock samples (Kern et al., 1993; Dortman, 1992), they cannot be interpreted directly in terms of rock composition.

As known, non-zero azimuths to noise sources may only increase apparent velocities (e.g. Sadeghisorkhani et al., 2016). Therefore, it is necessary to find another explanation of generally low S-wave velocities of the basement rock in our study.

The S-wave velocities in the uppermost crust down to several kilometres were previously evaluated using surface waves by Pedersen and Campillo (1991), Grad and Luosto (1992) and Grad et al. (1998). They reveal low values of S-wave velocity in the shallow crust and low values of a quality factor that is rapidly increasing at the depth of about 1 km. Moreover, Grad et al. (1998) found that the quality factor (Q factor) in the Archean shallow crust of the Fennoscandian Shield is lower than that in the Proterozoic crust. Grad and Luosto (1992) explained the low Q factor in the uppermost 1 km of the crust by increased crack density. Therefore, increased crack density can explain also generally low S-wave velocities of the basement rocks revealed by our study. Consider that the wave is propagating through fractured medium consisting of granites and the fractures are filled with some clastic rocks. If the S-wave velocities in the granitic rock are about 3300 m/s and those in the clastic rocks are about 400 m/s (like the velocity in quaternary deposits), then the averaged velocity would be about 1700–2000 m/s, depending on crack density. Taking the general effect of fracturing into account, we conclude that the values of S-wave velocities lower than 2000 m/s correspond to felsic rock, while the higher values of velocity correspond to mafic and ultramafic rocks (Fig. 18).

7 Conclusions

In our study, we used the data of ambient noise recorded during a short time period in a generally quiet area with a low level of anthropogenic noise. We showed that the noise was non-stationary and that the azimuthal distribution of noise sources was neither isotropic nor heterogeneous during the whole data acquisition period. In spite of that, we obtained good-quality EGFs, dispersion curves and the S-wave velocity models showing presence and thickness of quaternary sedimentary cover and velocity heterogeneities in the bedrock that agree well with the geological data. We explain this by the fact that the ambient noise recorded during the XSoDEx experiment contained a large proportion of diffused wavefield produced by scattering of plane waves from distant sources and ballistic waves produced by certain types of non-stationary sources on stochastically distributed heterogeneities in the uppermost crust. We demonstrated by numerical modelling that certain types of geological structures, particularly those composed of rocks with contrasting elastic properties, could scatter plane waves and ballistic waves from non-stationary sources and produce scattered Rayleigh waves. This scattered wavefield together with the ballistic wavefield from non-stationary sources can be used for EGF evaluation if the special algorithms are applied. One opportunity is to use the SNRS algorithm. Our result is a practical illustration of the conclusions about retrieval of EGFs from the scattered wavefield revealed previously by Wapenaar (2004) and Wapenaar and Thorbecke (2013). In certain geological areas, extraction of EGFs from the wavefield scattered at heterogeneities provides an opportunity to reduce the time for short-term passive seismic experiments.